Systematically and efficiently improving $k$-means initialization by pairwise-nearest-neighbor smoothing
Authors: Carlo Baldassi
Venue: Accepted by TMLR ()

Abstract: We present a meta-method for initializing (seeding) the $k$-means
clustering algorithm called PNN-smoothing. It consists in splitting
a given dataset into $J$ random subsets, clustering each of them
individually, and merging the resulting clusterings with the pairwise-nearest-neighbor
(PNN) method. It is a meta-method in the sense that when clustering
the individual subsets any seeding algorithm can be used. If the computational
complexity of that seeding algorithm is linear in the size of the
data $N$ and the number of clusters $k$, PNN-smoothing is also almost
linear with an appropriate choice of $J$, and quite competitive in
practice. We show empirically, using several existing seeding methods
and testing on several synthetic and real datasets, that this procedure
results in systematically better costs. In particular, our method
of enhancing $k$-means++ seeding proves superior in both effectiveness
and speed compared to the popular ``greedy'' $k$-means++ variant. Our implementation
is publicly available at \href{https://github.com/carlobaldassi/KMeansPNNSmoothing.jl}{https://github.com/carlobaldassi/KMeansPNNSmoothing.jl}.

Published in Transactions on Machine Learning Research (12/2022)
Systematically and efficiently improving k-means initializa- tion by pairwis--nerest-neighbo smoothing
Carlo Baldassi
cato.baldassi@hnibocoo....
Deportment of Computing Sciences, ArtLab, BIDSA
Bocconi Universily, Milan
ELLIS Scholar
Reviewed on OpenReview: https: // openreview. net/ formm? id= id=FFtFAg3peR
Abstract
We present a meta-method for initializing (scoding) the k-means clustering algorithm called PNN-smoothing. It consists in splitting a given dataset into random subsets, clustering each of them individually, and merging the resulting clusterings with the pairwise-neares neighbor (PNN) method. It is a meta- in the sense that when clustering the indi- vidual subsets any seeding algorithm can be used. If the computational complexity of that seeding algorithm is linear in the size of the data N and the number of clusters k. PNN- smoothing is also almost linear with an appropriate choice of J, and quite competitive in practice. We show empirically, using several existing seeding methods and testing on several synthetic and real datasets, that this procedure results in systematically better costs. In particular, our method of enhancing k-means++ seeding proves superior in both effective- ness and speed compared to the popular "greedy" k-means++ variant. Our implementation is publicly available at https:/ httph//w.cooncccaooaaaaasssss/es Satpaaldasss//eeeaaaPssnnoooooio
1 Introduction
The classical k-means algorithm is one of the most well-known and widely adopted clustering algorithms (Berkhin, 2006; Wu et al., 2008). Given N data points X' = (zi) where each point is D-dimensional, 24 E RD, and given an integer & 2 2, the algorithm aims at minimizing the sum-of-squarce-errross (SSE) cost, defined function of k centroids P (p:)=aI.N E (1,...,k) k)",, such:
(Ca)a=t.k (RR)*, and of a partition of the data points
N Σ 1=1
SSE (C,P:X)
||:
12
For fixed C, the optimal partition is obtained by associating each point to its nearest centroid, and conversely, for fixed P, the optimal centroids are given by the barycenter of each cluster. The k-means algorithm starts from an initial guess for the configuration and alternates optimizing C and P until a fixed point is reached. This alternating procedure, due to Lloyd (Lloyd, 1982), greedy and converges to a local minimum. Its computational cost is O (kND) per iteration; both optimization steps can be straightforwardly parallelized over N. It is notoriously sensitive to the choice of the initial configuration, i.e. the seeding, both in terms of the final value of the SSE and of the number of Lloyd's iterations required to converge. Several schemes have been proposed, with various degrees of complexity. An extremely basic and cheap option is to sample the initial centroids uniformly at random from X' (MacQueen, 1967). Other popular methods typically produce considerable improvement and also scale like O (kND), e.g. k-means++ (Arthur & Vassilvitskii, 2007) and maxmin (Gonzalez, 1985; Katsavounidis et al., 1994), among several others. Yet other methods have worse scalings and thus tend to dominate the computational time, e.g. the pariviis--narre--niihhe (PNN) method (Equitz, 1989) which scales like Ω (N2).

---

Published in Transactions on Machine Learning Research (12/2022)
In this paper, we propose novel scheme, PNN-smoothing, based on randomly splitting the dataset X' into J subsets, clustering them individually with k-means, and then merging the resulting Jk clusters following the PNN procedure until only k clusters remain: these constitute the new seed for Lloyd's algorithm. It is a meta-method, in the sense that it can use any seeding procedure for the subsets. We denote this with (NNT where INIT is any seeding algorithm. If INPT is O (kNV, then PNNS(INIT) is also almost linear, as long as we set J = O (VN/k). Our empirical tests indicate that PNNS(INIT) gives systematically better SSEs compared to INIT, at very small computational cost, even in parallel implementation. Its results are surprisingly good even when INIT is one of the worst seeding methods, random uniform initialization.
Throughout the paper, we focus exclusively on the effect of seeding on basic k-means, which can be regarded as a basic tool in the optimization of the SSE cost. Evaluating and contrasting the effectiveness of different methods is not straightforward. Our main metric will be the wall-clock computational time, measured using state-of-the-art implementation of all methods under uniform conditions. It is also desirable, however, to provide some hardware- and langnage-independent metric, beyond the asymptotic scaling analysis. It is often the case that computing distances wet vectors takes up the majority of the computational effort and thus that the number of distance computations correlates well with the running time, especially in high-dimensional cases (although of course various circumstances, espocially caching, can significantly affect the computational cost of distance computation in practice) (Newling & Fleuret, 2016). Indeed, several methods employed to accelerate the Lloyd's iteration procedure are explicitly aimed at reducing the number of distance computations (in the P-from-- step) at the cost of some additional bookkeeping (Kaukoranta et al., 1999; Elkan, 2003; Hamerly, 2010; Ding et al., 2015; Newling & Fleuret, 2016; Xia et al., 2022).
We thus define the number of normalized distance computations (NDC) as the total number of distance (or squared-distance) computations between D-dimensional vectors, wherever they may appear in a procedure, divided by Nk. The normalization makes this measure comparable across datasets, so that a value of 1 corresponds to one computation of the partition P from C, performed from scratch. Due to the variability in the optimization procedure, and the use of accelerators, it generally quite difficult to estimate the NDC required by a given method. However, for some seeding algorithms the nunber of NDC can be computed exactly, and for others it is rather straightffrward to at least provide some useful lower bounds. Note that our implementation of Lloyd's amgorithm starts by optimizing C from P, and thus we regard sceding algorithms as producing an initial partition P. This choice is justified by the fact that several seeding algorithms, even if they are based on picking some initial centroids, include the computation of P as a byproduct, and those that don't would need to perform this step at least once in any case"; therefore, this definition allows to compare all seeding algorithms consistently. As a consequence, all the seeding algorithms that we will consider in this work require at least NDC.
The rest of the paper is organized as follows. In sec. 2 we review prior literature and describe a few seeding methods that will be considered in the tests. In sec. 3 we describe and discuss in detail the PNNS(INIT) scheme. In sec. 4 we present and analyze detailed numerical results on several challenging synthetic and non-synthetic datascts. Sec. 5 has a final discussion.
2 Relation to prior works
As mentioned in the introduction, a large number of seeding schemes for k-means have been proposed. Extensive reviews and benchmarks can be found in refs. (Celebi et al., 2013; Fränti & Sieranoja, 2019). Here, we only cover a few selected ones, chosen on criteria of simplicity, popularity, similarity with our scheme, and effectiveness. Their scaling and NDCs characteristics (logether with those of PNNS(INIT) described in the next section) are also summarized in table 1.
UNIF. Uniformly sampling (prefcrably without replacement) i centroids from the datasct (MacQueen, 1967) is arguably the most popular method. We'll 'll call this seeding method UNIF. It's extremely simple, but its performance is generally very poor, even in moderately hard circumstances, as it often leads to poor local minima and long convergence times.
I We only consider the exact Lloyd's iterations here, in which case sub-linear (in N) seeding methods such as AF-KMC2 (Bachem al., 2016) don't provide an advantage.
2

---

Published in Transactions on Machine Learning Research (12/2022)
For the reasons explained in the introduction, although the cost of sampling the centroids is O (k), the overall cost of UNIF is still O (kND), entirely due to performing NDC.
MAXMIN. Another simple method (Gonzalez, 1985; Katsavounidis et al., 1994) goes under the name of "furthest point", or "maxmin", ", or "maximin". We will refer to the same variant that was used in Celebi et al. (2013); Fränti & Sieranoja (2019), and call it MAXMIN. It. consists in selecting the first. centroid at random from the dataset, after which the process is iterative and deterministic: at each step, each successive centroid is chosen as the furthest point from the centroids selected so far. More preciscly, it's the point that maximizes the distance from its nearest centroid: co argmax,ex minsce || coll?.
This algorithm scales as O (KND), D): it also computes the optimal partition (with respect to the chosen centroids) as a byproduct of the selection procedure: it requires precisely NDC. The results of Fränti & Sieranoja (2019), on synthetic datasets, report this method as being among the best of those that were tested. On the other hand, in Celcbi et al. (2013), in which an array of real datasots was also tested, the authors advise against this method. They suggest instead, among the algorit in the same complexity class, to use grredy-k-meanstae or Bradley and Fayyad's "refine" (both described below).
[G]KM++. The "k-means++ sceding method (Arthur & Vassilvitskii, 2007) can be regarded as a more stochastic version of MAXMIN. While the first centroid is also selected at random from the dataset, the remaining centroids are sampled from the dataset with a probability proportional to the squared distance from the closest centroids, More precisely, the probability of selecting a point x as the next centroid e, when a 2 is PP (c. = x (co)sea) x mins<c ||= a This procedure thus also computes the optimal partition as byproduct, like MAXMIN. It can be further extended in greedy manner: at. each step a 2 2, candidates are sampled according to the previous probability distribution, the new SSE (with clusters) is computed, and the candidate with the lowest SSE is chosen. The number of candidates per step is usually logarithmic in the number of clusters k; in all our tests, we have used s = [2 + logk] kj. We refer to the original variant as KM++ and to the greedy one as G.M++.
The computational complexity of KM++ is O (kND) and it requires NDC, like MAXMIN, while + scales like O Dogg k) and the NDC required are slightly less than the number of candidates s. In Celebi et al. (2013) GKM was reported as superior to KM++ and overall as one among the best lincar (in N) stochastic methods; in Fränti & Sieranoja (2019) the results of KM++ were considered comparable to slightly worse than MAXMIN; however, GKM++ was not tested.
REF(T Bradley and Fayyad's "refine" seeding algorithm (Bradley & Fayyad, 1998) is the one that most resembles our proposed scheme. Indeed, the initial step of the two methods is basically the same, i.. it consists in splitting the dataset into random subsets and clustering them individually with k-means, thus obtaining groups of & centroids. The crucial difference relies in the way in which these J solutions are merged, which in Bradley & Fayyad (1998) is referred to "smoothing" procedure: in the refine method, the whole pool of Jk centroids is used as a new dataset and clustered for J times. Each time, one of the previous centroid configurations is used as seed for the Lloyd algorithm. Out of the resulting J configurations, the one with the smallest cost (computed on the pooled dataset) is finally chosen. Originally, the authors used UNIF as the seeding method for clustering the subsets, but this can be trivially generalized to other methods. We thus consider a meta-method like PNN-smoo,hing and denote it with REF (INIT) where INIT is the sceding method for the initial step.
If the computational cost of INIT is linear, O OKND)), then REF(INTT) scales as O 0(&N + JKP) D). The algorithm requires J < N/k in order to perform the first step, and thus REF (INII) is always at most O D). In terms of NDC, the initial J clusterings overall require at least the same amount as INIT (e.g. if INFT is UNIF or MAXMIN) plus some additional ones for the Lloyd's iteration which are hard to estimate; then at least another Jk/N NDC are needed for the merging; finally, 1 NDC is required at the end; overall, the lower bound on the NDC is 1 + Jk/N more than INIT. In the original publication, Bradley & Fayyad (1998), only REF (UNIF) with J = 10 was tested. This was also the setting used in Celebi et al. (2013); Fränti & Sieranoja (2019). We also used the same value of J in our tests (the only exception being dataset for which N/k 10, in which case we used J 5), but we tested more INIT algorithms. As mentioned above,
Thhis is the default value used by the seikit-learn library (Pedregosa et al., 2011) and seems to work well; Celebi et al. (2013) used 8 log (k) (it's unclear if it was truncated or rounded).
3

---

Published in Transactions on Machine Learning Research (12/2022)
in Celebi et al. (2013) REF(UNIF) was found to be among the best methods, whereas in Fränti & Sieranoja (2019) it was shown to perform rather poorly on synthetic datasets.
PNN. Our method starts out identically to REF(INIT), but it uses the PNN procedure to merge the resulting J clusterings. This procedure was originally introduced in Equitz (1989) as a seeding algorithm. The original algorithm, which we call PNN, is hierarchical. Tt starts with N clusters, one cluster per data point, and then merges pairs of clusters iteratively until only & clusters remain. Merging two clusters means that the partition P is updated by substituting the two clusters with their union. The centroids set C. is also updated, by substituting the two starting centroids Ca and a with the centroid of the new cluster, Cnew (zoCa + >bc)) / (=a + 2b); where &a and 2n are the number of clements in each of the two original clusters. The algorithm is deterministic and greedy: the two clusters to be merged at each step are those whose merging will result in the smallest increase in the SSE cost. It is easy to see from eq. 1 that the cost increment of merging two clusters of sizes ze and 2, and with centroids Ca and a Aab = =0 ||ca - coll" / (za + 26). Thus both the pairwise merging costs and the new centroid can be computed using only the centroids and the cluster sizes.
The initial computation of all the merging costs Aab requires N (N 1) distance computations, thus O (DN?) operations. The computational complexity of merging step, assuming that we are going from ++1 clusters to & clusters, would be O (DR²) if performed straighuforwaardly, due lo the need to updale the Aub after each merge. However, in Franti & Kaukoranta (1998) it was shown that significant speedup can be obtained by considering that most cluster pairs are unaffected by individual merges, and the complexity can be reduced to O (DRTe) where Tg [1. is essentially the number of (potentially) affected clusters
and is generally much smaller than h. The number of distance computations of a merging step can easily be limited to exactly à by memoizing and keeping up-to-date all the distances; however, since T is usually very small, in practice this technique ouly helps the running time to a limited extent, and only in very high-dimennsioal cases, and thus for simplicity our implementation does not use it. The update operations can be straightforwardly parallelized over k which proves advantageous above a certain threshold. In the original PNN algorithm the merging step must be performed N - k times with k ranging from N - 1 to k. which amounts at roughly /2/2 distance computations (assuming N > k). Thus, the overall iomputational complexity is Ω (DN2) and the NDC required are about N/k. This is quite expensive for large datasets. On the other hand, the results in terms of the SSE objective are generally very good.
The PNN scheme was also employed the genetic algorithm of Fränti (2000), where however it was used as a crossover step rather than a sceding procodure. In that algorithm, two given configurations with clusters each are first merged into a single configuration with 2 clusters, which is then used as the starting point for the PNN iterative merging, until k clusters remain. The resulting cost is O (DR?T)). which (crucially) does not involve a factor of N thanks to the fact that only centroid computations are involved in the merge. Our sceding scheme is similar to this, in that we also use the itcrative PNN merging procedure with an initial number of clusters much smaller than N.
3 The PNN-smoothing scheme
In this section we describe in detail the (INIS seeding scheme and discuss its properties.
The inputs of the procedure are the same as for any other algorithm (the dataset X and the number of clusters k), plus the subset-seeding algorithu INIT, and one extra parameter p, used to determine the
Significantly more distance computations could be skipped by koeping lower and upper bounds based on the triangle inequality, similarly to the approach used to accelerate Lloyd's iterations in Elkan (2003) and others. This appears be promising optimization for high-dimensional data, which left for future work.
4

---

Published in Transactions on Machine Learning Research (12/2022)
Table 1: Summary of seeding methods characteristics: computational complexity (only N and k dependence, all methods are linear in D) and normalized distance computations. For REF(INIT) and PNNS(INIT) only lower bounds on NDCs are available since the expressions don't account for the internal Lloyd's iterations. The тx terms for (NNSS and PNN represent hard-to-estimant terms, upper-bounded by X but much smaller in practice. For G we set s [2 + log k.. For REF(INIT) we normally set J = 10. For (INS() we used the value throughout the main text.
method
complexity
NDCs
UNIF
kN kN
1 1 1
MAXMIN
KM++
KN
G
N
1/k + (k- 11 1)/K /k
(IF(INIT)
KN
NDC (INIT) + +1+ Jk/N
(NNS(INIT)
NNT.EN
2 NDC((INTT) + p/2 N/k
PNN
Nº
number of subsets J. The output is configuration to be used as starting point for local optimization. The high-level summary (see also the illustration of each step in fig. 1) is as follows:
Algorithm 1: PNNS seeding
Input: X, k, INIT, P
Set J [VpN/(2))
Cap the result [N/k].
2. Split X' into J random subsets (Aa) a=1.."
3. Cluster each Xa independently, using INIT for seeding followed by Lloyd's algorithm; obtain J configurations (C., Pa), with & clusters each.
4. Collect the J configurations (Ca, into single configuration (Co, Po) for the entire X, with kJJ centroids and clusters.
Merge the clusters of (Co, Po), two at time, using the PNN procodure, until & clusters remain; obtain set of k centroids C.
Compute the optimal partition P associated to C.
Output: PC
Next, we describe and discuss in more detail each step.
1. The parameterization in terms of rather than J ensures that the asymptotic behavior of the algorithm almost linear in N and k as discussed in the introduction (cf. point 5 below). The upper bound J necessary for the following steps and 3. The effective range of p, assuming for simplicity that N is divisible by k, is [2k/N, 2N/k). The lower bound leads to J PNNS(INIT==NIIT The upper bound corresponds to assumption that INIT invoked on & points will assign each one to its own cluster). In other words, by changing P we can interpolate between any sceding algorithm INIT and the pairwise-ncarest- neighbor algorithm. This explains why, as a rule of thumb, increasing improves quality at the cost of performance, although this is not strictly true in all cases.
and N/k and PNNS (INT(DN))=PNX (under the
which i is much smaller than The allowed range for p is quite wide under normal circumstances, N. Throughout the main text we use the value p = 1, which from our preliminary analysis seems result. in good trade-off in all cases and for all INIT algorithms. This is valid (Le. non-degenerate)
5

---

Published in Transactions on Machine Learning Research (12/2022)
Lloyd's alg,
Figure 1: Example of PNNS(INIT) in action. Here D = 2, N 72, & == and INIT=UNIF. The numbering of the steps follows algorithm (step simply yields J= = 3). The crosses represent centroids. After step 3 the subsets clustorings contain some clear mistakes, but after the merge (4) and the PNN procedure (5) the remaining centroids are close to the correct positions. After step 6 the seeding procedure completed; the final local optimization solves the problem.
choice whenever k < N/2, which is arguably always the case in realistic scenarios, thus making it a reasonable default value. Additional results with p p= 10 are reported in Appendix A.
We split the data as evenly as possible, i.c. we create N mod J subsets of size [N/J] 1 and J- N mod J subsets of size [N/J]. This is easy to implement efficiently by just constructing sorted list of indices, each index being repeated for the appropriate number of times, and shuffling it. Our preliminary testing showed that the algorithm is not sensitive to the details of the splitting procedure.
The computational cost of each individual clustering of one the subsets depends on the choice of INIT. Assuming that be linear, like in the examples mentioned in the previous section, this scales like N/J rather than N. Just like for refine, the NDC required for this step are least the same for INIT (thus at least 1) plus those for Lloyd's iterations, which are hard to estimate. This step is also trivially parallelizable.
In order to obtain the new configuration, we just take the union of the centroids, i.c. Co = U, č.; the partition Po would also be simply the union of the (Pa). with remapped indices, but since it is not even needed for the algorithm (only the cluster sizes are used) it can be skipped. It is also interesting to note that the new configuration will, in general, be nowhere near optimal for problem with kJ clusters, since points near each other will likely be assigned to clusters coming from different subscts. Nevertheless, al least under some favorable scenarios, themselves be approximately clustered into & groups (see fig. 1). This is the same intuition at the root of the refine method. The question then becomes how to best find a consensus configuration among the different results.
can expect that the centroids in Co may
The PNN procedure will start, the name implies, by merging the closest centroids (accounting for their associated cluster size). Tf the centroids in Co mostly clustered already, the procedure will

---

Published in Transactions on Machine Learning Research (12/2022)
likely pick first the centroids that appeared in multiple subset clusterings. Each time two centroids are merged, their associated size (and thus weight) increases, such that even if in the last stages some very sub-optimal partitions are merged with a large cluster, the centroid will be heavily skewed toward the latter.
The computational cost of the merging scales like O (D)(J))' analysis of the PNN procedure of the previous section. This is indeed consequence of our choice for the scaling of J. The additional factor TJR is hard to estimate; it is bounded by Jk = pNN/22 but in practice it appears to be quite small. This is the step that dominates the computational complexity of the whole algorithm; however, it is seldom the step that takes up the majority of computing time in practice when using p = 1. Indeed, the number of NDC can be estimated as just p/2 (where p/4 come from the initial step and p/4 from the merging process, analogously to the analysis presented for PNN).
T.K) O (DNNTJH), as per the
6.
At the end of the procedure the partition is recomputed (1 NDC). This is done to compute the SSE and provide a consistent starting point for Lloyd's algorithm (as discussed in the introduction), but it also means that outlier points in the partitions that could be (virtually) produced during the PNN merge are eliminated (like the blue points in the bottom-right cluster and the green points in the top-left cluster in fig. 1).
Overall, the number of NDC in the PNNS(INIT) procedure is larger by at least + p/2 compared to INIT (roughly, see steps 3, 5 and 6). This estimate appears to capture the biggest time penalty of PNNS(INIT) in practice, even though it does not account for the additional Lloyd's iterations performed in step 3. This is because those usually require much less than 1 NDC per iteration thanks to the use of accelerator algorithms; furthermore they are compensated by obtaining a seed which is closer to a local optimum, so that the overall number of Lloyd's iterations at the end of the process is comparable (the data showing this is reported in Appendix A).
4 Numerical experiments
4.1 Experimental setup
We performed series of tests comparing PNNS with all the algorithms mentioned in sec. 2, namely: UNIF, MAXMIN, KM++, GKM++, PNN, and REF(INIT) and PNNS(INIT) with INITE {UNIF, MAXMIN, KM++, GKM+ We used the same data structures and programming language (Julia v1.7.3) for all of them. For the local optimization part, ie. Lloyd's algorithm, we implemented a number of techniques that can accelerate the computation, while keeping exact, by skipping some updates at the cost of some bookkeeping: the "reduced computation" method (RC) by Kaukoranta al. (1999); the Elkan method (ELK) from Elkan (2003); the Hamerly method (HAM) from Hamerly (2010); the Yinyang method (YY) from Ding et al. (2015); the exponion method (EXP) from Newling & Fleuret (2016); the ball-kmeans method (BALL) from Xia et al. (2022). For ELK and YY we implemented the simplified versions described in Newling & Fleuret (2016). For each of the datasets that we tested we chose the accelerator technique that resulted in the fastest average convergence time, when using KM++ for scoding, across at Icast 30 random repetitions, and used that accelerator for all other tests with that dataset. The winner for each dataset reported in Table 2. This choice puts PNNS the maximum disadvantage, since, as a general tendency, the better the accclerator method, the larger the relative penalty of applying the PNNS scheme. 4 run until convergence to fixed point. Our code is available In all our tests Lloyd's algorithm saaarraalllass/MeeeeanPSSSnnnnnn
All the timings that we report refer to tests performed on the same hardware (Intel Core 17-9750H 2.60GHz CPU with 6 physical cores, 64Gb DDR4 2666MHz RAM, running Ubuntu Linux 20.04 with 5.15.0 kernel)
*This is because, as mentioned above, the number of overall Lloyd's iteration ends up being comparable between INIT and PNNS(INIT), but accelerators become more effective the iterations progress, and PXNS Lloyd's scheme is effectively started twice (this effect indeed captured by the NDC count discussed in sec. 3, point 6).
7

---

Published in Transactions on Machine Learning Research (12/2022)
Table 2: Characteristics of the datasets used in the tests and best accelerator for each dataset.
dataset
D
NN
A:
accel.
A8
2
7500
50
EXP
synthetic
Birch1
2
100000
100
EXP
Birch2
2
100000
100
EXP
Unbalance
6500
EXP
Dim1024
1024
1024
16
EXP
Bridge
16
4096
256
RC
real world
House
3
34112
256
EXP
Miss America
16
6480
256
RC
UrbanGB
2
360177
469
EXP
Olivelti
4096
400
40
ELK
Isolet
617
7792
26
ELK
USCensus
68
2458285
100
YY
with no other computationally intensive processes rumning while testing. All codes were carefully optimized and can run in parallel with iterations, and over the & centroids when k 2 500 during the PNN procedure, see sec. (2)); most of our results are shown for the single-threaded case except where otherwise noted. We also report the average total number of normalized distance computations (NDC), summing up those performed during seeding (discussed in the previous sections) and in the final optimization phasc.
(the parallelization is over the N data points during Lloyd's
We tested a number of synthetic and real-world datasets whose characteristics are shown in table 2. The synthetic datasets are mainly intended to measure the ability of the algorithms find the solution when one can be clearly identified, and for direct comparison with the results of ref. (Fränti & Sieranoja, 2019). In these kind of datasets, when k-means gets stuck in a sub-optimal minimum, it is usually due to having made one or more clearly identifiable mistakes (see e.g. the supplementary materials of ref. (Baldassi, 2022)). As we shall show in sec. 4.2, PNNS(INIT) is the only family of algorithms (among the linear or quasi-linear ones) that is capable finding the solution in 100% of these cases; this gives a degree confidence that, when using PNNS, k-means will not get stuck into clear and easily-fixed mistakes, and that any variation in the optimization result more likely to emerge from the features of the datasets not being aligned with the SSE objective.
The real-world datasets on the other hand do not generally have simple structure or a simple "solution", and their optimum SSE is unknown. Our tests in sec. 4.3 will show that the PNNS scheme can provide systematic improvements in the SSE with a bounded increase in computing time, generally offering a better trade-off than alternative methods.
4.2 Synthetic datasets
In Fränti & Sieranoja (2019), several synthetic datasets with different characteristics were chosen and tested in order to probe the strengths and weaknesses of several sceding algorithms with respect to propertics of the data. We picked of the most challenging ones, all obtained from the UEF repository (Fränti & Sieranoja, 2018): A3, Birch1, Birch2, Unbalance and Dim1024 (sce table 2). We scaled each datasot uniformly order to make them span the range [0, 1) The difficulty for A3, Birchl and Birch2 is in their relatively large size and abundance of local minima; for Unbalance, it's the fact that some clusters are small and far from the bigger ones; for 10m1 it's the large dimensionality. All algorithms tested in Fränti Sieranoja (2019) showed poor results in at least some of these datasets; in particular, the authors report 0% success rate (as defined below) on Birchl for all algorithms.
OOur implementations of ELK, HAM, and EXP generally roughly comparable the very optimized ones (written in C+++)) that accompany Newling & Flcuret (2016), available https:// /gggtuh.hh.oom.m www.cooooiiioaaeeeeeeee htt//diaa/cakmmassss ccakmcans; our implementation of BALL gen- erally faster than the C++ one provided by the original authors of Xia et al. (2022), available htttps//gitthbb.c. httg//gttthhgomm//g /syxiaa/yaaaaaa /ball- l-means.
8

---

Published in Transactions on Machine Learning Research (12/2022)
Table 3: Results on synthetic datasets
success rate
A3
Birchl
Birch2
Unbalance
Dim1024
UNIF
0
0.001
MAXMIN
0.004
0
0.223
KM
O
0.541
0.996
G
0.051
0.02
0.06
0.946
REF(UNIF)
0.055
ARF(GKM++)
0.239
0.03
0.44
PNNS(UNIT)
0.709
0.26
0.696
0.919
(MAXMIN)
1
PNNS(KM++)
0.98
PNNS (GKM
1
1
PNN
convergence time (meanttstdvv)
A3
Birchl
Birch2
Unbalance
Dim1024
(in 10 10-
(in 107
(in 10 10"1x)
1s)
(in 10 10)
5.1 5.111. 1.5
2.3 0.6
0.63 1 0.10
3.1 1 1.8
5.5 1: 1.3
MAXMIN
4.0 1 0.9
2.0 0.4
0.75 + 0.10
0.7 = 0.8
6.4 + 0.8
KM++
4.1 + 1.2
1.78 + 0.38
0.66 + 0.10
1.1 + 1.3
5.4 + 1.5
+
6.5 1.4
2.44 L026
1.64 0.11
0.9 1.0
12.8 1 1.7
REF(UNIF)
14.7 2.5
2.97 + 0.33
1.39 + 0.09
4.7
20.5 +
REF(GKM
13.8 2.2
3.06 0.23
2.06 1 0.12
19.1 +
PNNS(UNIP)
11.2
1.78 0.11
3.7
14.2
(NNS(MAXMIN)
8.7
2.51 0.10
1.60 + 0.09
1.8 1.4
12.9 + 2.5
P)NS(KM++)
9.1
2.55 4 0.07
1.59 + 0.07
2,0
13.1 + 2.8
PNNS(GKM)
12.1 2.0
3.32 1 0.08
2.49 0.09
2.5 1 14
19.9 1 3.3
PNN
396 = 16
660 + 6
709.6 + 3.1
292 = 11
1678
+ 34
normalized distance computations (mean
Ay
Birch1
Birch2
Unbalance
Dim
UNIF
0.18
2.52
1.146 0.027
0.6
1.89 0.32
L.33 0.09
2.08 0.29
1.073 + 0.020
1.26 1.26 0.10
1.07812 0.0
KM+
1.41 + 0.12
2.00 0.25
0.25
1.063 0.001
1.35
1.079 0.009
GRM
5.08 1 0.05
6.46 0.16
5.972 L 0.004
1.111
3.89062 0.0
REF(UNIF)
3.94 0.12
0.22
2.393 + 0.019
4.6
6.6 1 0.4
REF (GGM++++))
6.992 + 0.029
0.12
7.1128 + 0.0030
5.068 +: 0.028
6,9261 +: 0.0014
PNNS(UNIP)
3.84 = 0.06
4.07 0.04
3.191 0.012
4.63 - 0.37
4.27 = 0.16
(MAMMM)X
3.283 1 0.025
3.662 0.025
3.044 0.007
3.257 1 0.026
3.380 0.027
+)NS(KM+++
3.439 + 0.039
4.033
3.084 0.008
3.40 0.05
3,380 0.027
PNNS(GKM)
7.137 1. 0.022
8.416 1 0.020
7.957 0.007
5.903 1 0.034
6.193 0.027
PNN
262.324
1730.39
1738.25
1422.21
531.759

---

Published in Transactions on Machine Learning Research (12/2022)
Table 4: Results on real-world datasets
mean SSE cost (mean-stder 100 repetitions)
Bridge
House
Miss A.
Urb.GB
Olivetti
Isolet
(x10")
(x10²)
(x10")
(x10²)
(x10*)
(x10")
UNIF
1.178 + 0.009
10.11 = 0.13
6.07 = 0.05
6.8 1.1
1.296 + 0.027
1.196 + 0.010
MAXMIN
1.139 E 0.005
10.18 + 0.07
5.794 +: 0.028
2.99 + 0.07
1.252 0.015
1.233 + 0.014
KM++
1.154 0.007
9.60 0.04
5.687 0.037
2.72 0.06
1.276 0.021
1.196 = 0.010
1.124 0.004
9.529 i 0.027
5.307 0.016
2.430 0.020
1.227 0.013
1.190 1 0.007
REP(UNIF)
1.162 0.006
9.95 0.10
5.850 0.035
5.3 0.4
1.298 0.023
1.189 + 0.006
RRF (GKM+++
1.159 0.005
9.562 = 0.027
5.826 0.031
2.420 + 0.022
1.255 0.020
1.1837 0.0036
(UNS(UNIF)
1.124 1 0.005
9.553 1 0.034
5.607 1 0.040
2.70 1 0.06
1.215 1 0.014
1,1800 1 0.0019
()MS(KM++)
1.1076 + 0.0033
9.486 + 0.020
5.403 + 0.013
2.323 + 0.010
1.209 + 0.013
1.1795 + 0.0016
PNNS(GKM 1)
1.0947 + 0.0029
9.476 + 0.020
5.342 + 0.010
2.297 1.000
1.189 0.007
1.1790 0.0015
PNN
1.08279
9.49701
5.31588
2.3153
1.16238
1.17692
minimum SSE cost (over 100 repetitions)
Bridge
House
Miss A.
Urb. GB
Olivetti
Isolet
(x10")
(x10²)
(x10º)
(x10?)
(x10*)
(x10*)
UNIP
1,15723
9.87713
5.95536
5.42031
1,23747
1.1795
MAXMIN
1.12869
10.0159
5.72134
2.90322
1.21735
1.20639
KM++
1.14219
9.52638
5.59949
2.60788
1.23035
1.18025
G
1.11421
9,45992
5.47583
2.37869
1.19908
1. 1.17709
REP(UNIP)
1.14662
9.77128
5.76604
4.39285
1.24007
1.17688
+G
1.14344
9.49776
5.74112
2.37821
1.2116
1.17737
(UNIS
1.11345
9.47461
5.51748
2.5921
1.17155
1.17685
(NNS(KM++)
1.10041
9.44324
5.37209
2.30731
1.18579
1.17687
+NNS(GKM+++)
1.08729
9.42885
5.31768
2.28822
1.17066
1.17695
PNN
1.08279
9.49701
5.31588
2.3153
1.16238
1.17692
convergence time (meannstdev))
Bridge
House
Miss A.
Urb. GB
Olivetti
Isolet
10
(in 10 10
10 s)
(in s)
10 Ps)
(in 10 10=1s)
UNIP
4.17 = 0.38
0.6
1.13
1.63 = 0.23
2.7 0.6
1.52 0.30
4.6 + 0.5
3.5
1.78 0.33
1.73 = 0.20
4.1 0.6
1.58 : 0.21
KM++
4.22 = 0.36
2.40 0.37
0.11
1.56 + 0.13
4.10 = 0.39
1.75 0.29
GKM
6.45 1 0.35
3.33 1 0.33
5.51 0.16
8.8 0.6
3.02 1. 0.29
REP(UNIF)
22.0 - 0.5
6.5
3.26
3.25
8.0
1.68 - 0.24
RRF(GKM++)
24.6 - 0.6
6.48 0.36
0.99
5.36
12.0
2.19 0.22
PNNS(UNIP)
5.90 = 0.31
4.74
0.09
4.12
6.2
1.63 0.16
PNNS(KM++
4.44 0.31
1.10
3.58 0.16
6.74 0.38
+NNS(GKM+++)
8.13 = 0.39
5.22 0.32
1.89 0.11
5.99 0.14
9.51 0.39
2.07 0.18
PNN
36.2 1.2
93.68 = 0.20
9.356 0.033
819
30.5 = 1.8
297.2 + 2.8
normalized distance computations (mranistdov)
Bridge
House
Miss A.
Urb. GB
Olivetti
Isolet
UNIF
8.4 0.5
5.7 = 0.9
13.5 4.9
1.60 + 0.23
1.79 0.14
3.35 + 0.36
MAXMIN
8.9 + 0.8
0.6
+
1.196 0.087
2.21 0.07
3.14 0.25
KM++
8.2 0.6
3.71 0.39
13.7 1.3
1.229 0.029
2.26 = 0.07
3.36 0.35
GKM++
13.2
0.32 8.3
1.5
8.168 0.033
6.04 = 0.05
6.92 0.33
REF(UNIF)
44.8 1 0.8
1.04
3.15 0.15
5.68 1 0.10
4.62 0.26
REF REF (GKM++)
49.5 = 0.8
12.25 + 0.35
48.3
9.465 0.002
10.34 0.07
8.20 0.22
(NNS(UNIF)
10.81 + 0.37
6.60 + 0.36
16.6
4.03 0.09
4.25 0.10
5.06 + 0.15
PNNS(KM++)
10.6 = 0.4
6.05 + 0.29
16.7
3.31 + 0.06
5.00 = 0.10
5.24 + 0.17
PNNS(GKM ++)
15.9 1 0.4
11.59 1 0.30
22.2 11.
10.191 1 0.014
8.76 1 0.09
8.73 0.16
PNN
51.8177
271.383
87.8464
1351
19.7431
739.388
10

---

Published in Transactions on Machine Learning Research (12/2022)
All of these datasets were generated from isotropic Gaussians centered around ground-truth centroids, and for all of them the global optimum of the SSE is very close to the ground truth. Under these circumstances, it is reasonable to classify the local minima configurations that the algorithms produce by their "centroid index" (CT), as defined in Fränti et al. (2014). The CI is computed by matching each centroid of a clustering with its closest one from the ground truth, and counting the number of unmatched ground truth centroids; in formulas:
Aet Σ11
CI
ta E (l....,k)
.b3 : b arg min d (Ca: (ca:cg)) )
arg
(2)
where of with b E {1, {1...., bay are the ground-truth centroids (in our setup k = h*t), and 1 (-) is an indicator function. The CI can be interpreted as the number of mistakes in the resulting clustering, and therefore we define the success rate of an algorithm as the frequency with which it finds a solution, i.c. a configuration with CI= O.
We present our most representative results in table 3: the complete results are reported in Appendix A.1. All algorithms except for PNN and those in the PNNS family fail badly in at least some datasct. Conversely, all the PNNS(INTT) algorithms solve all the datasets in 100% of the cases (98% for PNNS (KM+++) on A3), with the only notable exception being INIT=UNIP. Even in that case, however, the performance is much better than UNIF and REF(UNIF). and at times even than REF +CF())+++) (which is the variant in the refine family that gives the best overall results). This demonstrates that (at least for this scenario) PNN-smoothing is a considerably better smoothing technique than refine.
The timings of the PNNS(INIT) algorithms are generally comparable with, and often better than, those of the refine family, and within a constant factor (smaller than 3) of the basic methods; on the other hand, the only other algorithm capable of solving all the datasets, PNN, is Ω (N2) and indeed orders of magnitude slower. These conclusions are corroborated by inspecting the NDC, which depend on the accelcrator used but not on the implementation or the hardware, and are roughly correlated with the running time.
Overall, our results show that the PNNS(INIT) scheme is able to improve the success rate of any INIT algorithm at the cost of a modest time penalty, and that it is a considcrably better at achioving this than the alternatives. It is particularly interesting to consider the case of the Birch1 dataset, for which UNIF, MAXMIN and KM++ all have a 0% success rate. Both GKM++ and REF (KM++) are intended to improve KM++ at the cost of additional computations, and REF(GKM++) combines the two approaches; yet, none of them improves the success rate beyond 3%. The +)M+(KM++) method, on the other hand, achieves 100% success rate Birch1 in the same smaller amount of time as those other methods (and fewer NDCs).
4.3 Real-world datasets
We tested 7 challenging real-world datasets (see table 2). The first three, Bridge, House and Miss America, were obtained from the UEF repository (Fränti & Sicranoja, 2018); UrbanGB, Isolet and USCensus are from the UCI repository (Dua & Graff, 2017) and Olivetti from scikit-learn (Pedregosa et al., 2011). All of these are comparatively large and quite challenging for SSE optimization. The Isolet dataset was chosen as the hardest one among those tested in Celebi et al. (2013) (based on the results reported there); the USCensus dataset was chosen for its large size, to provide a test of parallelization efficiency; the other 5 datasets were also tested in Baldassi (2022), where was shown that even sophisticated, stale-of-the-art evolutionary algorithms cannot casily find their global minima (and indood it is not even clear if they can find them at all). We did not scale the datasets, except UrbanGB for which we scaled the longitude by a factor of 1.7 to make distances roughly proportional to geographical distances, and USCensus for which we linearly transformed each dimension individually to make them fit into the range [-1, For all datascts, except UrbanGB and Isolet, the choice of k follows existing literature and is otherwise arbitrary.
We performed 100 tests for each dataset (30 for UrbanGB and USCensus) and computed the average and minimum SSE cost achieved across the runs (the latter metric can provide an indication about what can be achieved with a multiple-restarts strategy by each algorithm), as well as the average running time, NDCs, Lloyd's iterations. As for the previous section, we report our most representative results here; the full results
11

---

Published in Transactions on Machine Learning Research (12/2022)
Baslge
x10
10.3
1.36
10.2
10.3
1.1
9.9
KM
GRO+++
1.11
MAXMIN REF UNNSONIT)
9.7
1.00
9.3
0.00
0.1
0.15
0.2
0.25
9.1
02
03
0.4
0.3
0.5
0.7
tinne (4)
time (6)
X10
X10
3.8
1.2
5.
1.24
5.6
*
1.20
1.18
0.00
0.1
0.15
time
time
Iecbst
UihanGGB (am-peralad))
10
1.54
3.2
2.9
2.8
SSC
1.21
X
2.6
1.29
2.4
1.17
2.2
119Cazam (paralid))
UrlanCB (paraDd)
1.35
3.2
2.9
SSE
2.7
2.5
X
1.15
23
1.10
2.2
5
10
15
20
30
10
45
9.5
2
25
3.3
4
40
Hime (4)
Lnne (4)
Figure 2: SSE cost vs convergence time, averages over 100 samples (only for UrbanGB and USCensus), for the real-world datasets. The meta-methode symbols enclose the symbol for their INIT algorithm. The dashed lines the costs achieved by PNN, whose timings would be off-scale. The top 6 panels non-parallel test: the bottom row are tests performed with 4 threads in parallel.
12

---

Published in Transactions on Machine Learning Research (12/2022)
can be found in Appendices A.2 and A.3. We first present the results for the non-parallel case (on all datasets except USCensus). Some results are presented in table 4, confirming that the PNNS family of algorithms attains results superior to all other linear algorithms, in comparable times (and comparable NDC), both in terms of the average and of the minimum SSE cost. More specifically, (INT(INIT) consistently achieves better costs than both INIT and REF(INTT), for all dat and all the lested INIT (note that REF(INIT) is not consistently better than INIT). This is confirmed by a statistical significance analysis (standard two-sided Wilcoxon rank-sum test, all p-values smaller than 10-11;; in fact, in nearly all cases the average cost for PNNS((NIT) is even smaller than the minimum cost for INIT. In many cases, the average time for PNNN(INIT) is shorter than for (IF(IN PNNS (GKM++) generally achieves the best SSEs (possibly on par with other PNNS methods), and in 2 cases out of 6 it even outperforms PNN, which belongs to a different computational class and is considerably slower.
In fig. 2, we plot the mean SSE vs convergence time for the best linear "plain" methods, KM++, GKM++ and MAXMIN, their REF (INIT) versions", and their (NNS(INIT) versions. They all have comparable timings, but the (NNS(IN family is consistently below the others. The dashed lines denote for reference the costs achieved by PNN (its timings would all be off-scale, ef. table 4), showing that even when it performs better than the algorithms in the PNNS family, the difference is generally relatively small.
In order to test the effect of parallelizatio ion on the performance of the algorithm, we also performed some tests with tulti-threading enabled, using threads. We tested the two largest datasets, USCensus and UrbanGB. The results are shown in the last two rows of fig. 2 (complete data in Appendix A.3); the case of UrbonGB allows a comparison with the non-parallel ones; note that we did not run PNN on USCensus because it would be impractical. The results are qualitatively similar to the non-parallel versions.
As for the synthetic datasets, it is particularly interesting to compare the different existing methods that can be used to improve on KM++ namely CKM++, REF +KM+++) and REF (EF(GKM++), with our method (N Our method achieves better SSE costs than all of the others; in 5 out of 7 datasets it is also faster (for House only GKM++ is faster, for UrbanGB only REF(KM++) is faster). Even in this real-world setting, our PNNS scheme provides better trade-off between cost and computational time than the existing alternatives.
5 Discussion
We have presented a scheme for k-means seeding called PNN-smoothing that can be applied to any existing linear (0 (IND)) algorithm, with a limited impact on the scaling and on the convergence times. Our experiments, performed with an efficient implementation and state-of-the-ar techniques, show clear and consistent improvements on challenging synthetic and real-world datasets, and systematically superior results (both in terms of quality and of speed) with respect to the similar "refine" smoothing scheme (also note that to the best of our knowledge we were the first to report tests for REF(INIT) with /NITUUNIF, finding that it does not systematically improve over INTT). One particularly interesting case that of the very popular k-means++ seeding algorithm: we showed that our scheme outperforms alternative enhancement techniques (i.c. making it & "groedy", using "refine", or both) in terms of both quality and speed." 7
The overall picture is unchanged in a parallel multi-threading context. We also verified that, as one would expect, the results are qualitatively the same regardless of the scheme used to accelerate k-means iteration. The experiments also indicate that PNN-smoothing is not particularly susceptible to the original seeding scheme.
Overall, our results do not highlight any clear best among the PNNS algorithms: PNNS(KM+++) is ususally the fastest, and it produces consistently good results, but PNNS (MAXMIN) is not much different, while the best costs are usually obtained by PNNS(CKM++)) When trying to improve the cost sacrificing performance in difficult cases, the best strategy is probably to use PNNS (ONN(G with a larger From the results reported above, one can obtain an estimate of the potential gain of increasing p by looking at the differonce
For the Olivetti datasct N/k = 10 and thus we used J 5 instcad of the default J 10. 7 At the time of writing, +eeans+- is the default seeding methad the very popular seikit-learn package; according to our results, PNNS(KM+++) would be a superior default.
13

---

Published in Transactions on Machine Learning Research (12/2022)
between the costs achieved with p == 1 and those achieved with PNN. Since this difference is rather small, one would not expect large improvement (in Appendix A we report the results with p 10). Another option would be to apply the PNN-smoothing scheme recursively (e.g. PNNS(PNNS(KM+H) we thoroughly explored this possibility, 9 but concluded that does not bring significant advantages over simply using a larger p. On the other hand, minimizing the costs of paramount importance, there exist better optimization schemes than simple k-means, most notably evolutionary algorithms that use k-means as their starting point (Fränti, 2000; Baldassi, 2022); for those, the seeding process comparatively less important in general, but it still be crucial in some cases (as discussed in ref. (Baldassi, 2022), when trying to access low-cost regions of the configuration space, improving the seeding is usually more efficient than blind exploration via random mutations). Assessing the effect of better seeding on algorithmic schemes that go beyond simple k-means (like the above-mentioned evolutionary schemes, or for scenarios in which full Lloyd's iterations are too costly, or data cannot fit into memory, etc.) is left for future work.
and similar):
References
David Arthur and Sergei Vassilvitskii. k-means+ The advantages of careful seeding. In Proceedings of the eighteenth annual ACM-SIAM symposium on Discrele algorithms, pp. 1027 1035. Society for Industrial ottp://r....aaaoooodd///////////
and Applied Mathematics, 2007. URL http://theory.......// httttteerrssstaffrdoooee wew.eseegei/ppperreeeeee htwrrg///pppessstlllaass-----ooo /kteeassPPPPPso...... pdf.
Olivier Bachem, Mario Lucic, Hamed Hassani, and Andreas Krause. Fast and provably good seedings for k-means. In D. D. Lee, M. Sugiyama, U. V. Luxburg, I. Guyon, and R. Garnett (ods.), Aduances in Newral Information Processing Systems, volume 29, pp. 55 63. Curran Associates, Inc., 2016. URL http::/papeee
Carlo Baldassi. Recombinator---mmeann An evolutionary algorithm that exploits k-means++ for recombina- tion. IEEE Transactions Evolutionary Computation, pp. 1-1, 2022. URL https:///oi.orr/10 1109/
Pavel Berkhin. A survey of clustering data mining techniques. In Grouping multidimensional data, pp. 25 71. Springer, 2006. URL tttps://doi:
Paul S Bradley and Usama M Fayyad. Refining initial points for k-means clustering. In ICML, volume 98, pp. 91-99. Citeseer, 1998. URL B528arep=reppttyee=pdt
M Emre Celebi, Hassan A Kingravi, and Patricio A Vela. A comparative study of efficient initialization methods for the k-means clustering algorithm. Expert systems with applications, 40(1):200 210, 2013. eswa. 2012.07 07. 021.
URL https: //wi. http://t....o...../////
Yufei Ding, Yue Zhao, Xipeng Shen, Madanlal Musuvathi, and Todd Mytkowicz. Yinyang k-means: A drop-in replacement the classic k-means with consistent specdup. In Francis Bach and David Blei (eds.), Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 579-587, Lille, France, 07-09 Jul 2015. PMLR, PMLR. URL https: / //proceedings. mlr. www.ppprsss///eeee /37ei// html.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http::/archiivee ics. www.iii.......... m1.
Charles Elkan. Using the triangle inequality accelerate k-means. In Proceedings of the 20th international conference on Machine Learning (ICMI-03), pp. 147-153, 2003. URL http://www. aaai aaai.aaa.aaaaaa.aaa....aa.aaa org/Papers/ ICML/2003/ICLL03-0222 pdf.
Indeed, results also be summarized by stating that PNNS can provide solutions of comparable quality to PNN but with a quasi-lincar scaling rather than quadratic one.
We also considered fully recursive version, which we call PNNSR, that calls rNNs recursively until the subsets have size N < 2k, at which point uses PNN. It. achieves good costs, but it scales at least as O (Nlog log N). See the associated code for implementation details.
14

---

Published in Transactions on Machine Learning Research (12/2022)
William H Equitz. A new vector quantization clustering algorithm. IEEE transactions acoustics, speech, and signal processing, 37(10):1568 1575, 1989. doi: 00.1109/29.3556 URL https://doi. http::://ooo..o//111110000
29.35395.
Pasi Fränti. Genetic algorithm with deterministic crossover for vector quantization. Pattern Recognition Letters, 21(1):61-68, 2000. URL 44t (99))60333))))
Pasi Franti and Timo Kaukoranta. Fast. implementation of the optimal pmm method. In Proceedings 1998 International Conference on Image Processing. ICIP98 (Cat. No. 98CB33266), pp. 104 108 vol.3. IEEE,
1998. URL https://doi.org:: 11w.:11119IIIII.1111 999001.
Pasi Fränti and Sami Sicranoja. K-means properties six clustering benchmark datascts. Applied Intelli- gence, 48(12):4743 4759, 2018. URL http://cs. uef wtspp//aaaaatsss.///
Pasi Fränti and Sami Sicranoja. How much k-means can be improved by using better initialization and repeats? Pallern Recognilion, 2019. URL https://doi org/ 1016/ wwt.oat.c....... 2014
2019. 04. 014.
Pasi Fränti, Mohammad Rezaei, and Qinpei Zhao. Centroid index: cluster level similarity measure. Pattern Recognition, 47(9):3034 3045, 2014. URL
patcoa............. 20.4.03.00 017.
Teofilo F Gonzalez. Clustering to minimize the maximum intercluster distance. Theoretical Computer Science, 38:293 2009 1985. URL https://doi. http:////ooooorg/1...
http::t:://o...or//10..0000 2ttp:/00.1006/0004-3337770 (85) 90224-5.
Greg Hamerly. Making k-means even faster. In Proceedings of the 2010 STAM international conference on data mining, pp. 130- 140. SIAM, 2010. URL
11 97881171.97111117781111111 12.
Ioannis Katsavounidis, C-C Jay Kuo, and Zhen Zhang. A new initialization technique for generalized lloyd iteration. IEFE Signal Processing Letters, 1(10):14-146, 1994. URT. http::::/doo.oo//////.sss http://dd...rr/0.111////9/999/// 329844.
Timo Kaukoranta, P Franti, and Olli Nevalainen. Reduced comparison search for the exact gla. In Proceedings DCC'99 Data Compression Conference (Cat. No. PR00096), pp. 33 41. IEEE, 1999. URL https://doi. org/10 0ts/10.11000 1009/10CD 1999. 755651.
Stuart Lloyd. Least squares quantization in pem. JEEE transactions on information theory, 28(2):129-137, 1982. URL https:///ol org. 1109/T11.. 1982. 1056489.
James MacQueen. Some methods for classification and analysis of multivariate observations. In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 1: statistics, volume 1, pp. 281 297. Oakland, CA, USA, University of California Press, 1967. URL org/ org/euclid. www.p/ 120051/2005
James Newling and François Fleuret. Fast k-means with accurate bounds. In Maria Florina Balcan and Kilian Q. Weinberger (eds.), Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pp. 936-944, New York, New York, USA, 20-22 Jun 2016. PMLR, PMLR. URL https://prrceeeding: .mw.g.ssssrr....m... wwt.:/pesss////eeeeeeee http////////eee/n///em html.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Bruchor, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825 2830, 2011. URL http://www.jml wwtgppaers///eeeooooot.....seeee aaapprss/v2/ppdregeaaaaa html.
Xindong Wu, Vipin Kumar, J Ross Quinlan, Joydeep Ghosh, Qiang Yang, Hiroshi Motoda, Geoffrey J McLachlan, Angus Ng, Bing Liu, S Yu Philip, et al. Top 10 algorithms data mining. Knowledge and information systems, 14(1):1 37, 2008. URL https: http:://....ooo....//////
org/10. 00000/011100//100100000000-----0
Shuyin Xia, Daowan Peng, Deyu Meng, Changqing Zhang, Guoyin Wang, Elisabeth Giem, Wei Wei, and Zizhong Chen. Ball k-means: Fast adaptive clustoring with no bounds. JEEE Transactions on Pattern Analysis and Machine Intelligence, 44(1):87-99, 2022. URL sttps://doi.org/11 1109/TPAMI 2020. 3008694.
15

---

Published in Transactions on Machine Learning Research (12/2022)
A Appendix: Complete numerical results
In this section we report the full data for all numerical lests of sec. of the main text. This complements tables 3 and 4 the main text, and includes additional measures. When we report the average number of Lloyd's iterations, we include the ones preformed during soeding by REFINE(INIT) and PNNS(INIT), but we scale them by a factor of J keep them comparable to the ones performed during the final optimization.
Whenever report means, also report the standard deviations.
A.1
Synthetic datasets
We tested 1000 samples for A3, Undalance and Dim Pim1024, and 100 samples for Birch1 and Birch2.
A.1.1 Convergence time (in s)
seeder
A3 (x10-9)
Bircht (x10-1)
Birch2 (x10-1)
Unbolance (x10-3)
Dim Rkm/024 (x10-3)
UNIF
5.1 1 1.5
2.3 0.6
0.63 1 0.10
3.1 1.118
5.5 1 1.3
MAXMIN
4.0 0.9
2.0 0.4
0.75 + 0.10
0.7
+
KM++
4.1 1.2
1.78 0.38
0.66 + 0.10
1.1
5.4
GKM I-+
6.5 1.4
2.44 1 0.26
1.64 1.64 0.11
1.0
12.8 1 1.7
REF (UNIF)
2.5
2.97 + 0.33
1.39 +.30
1.5
20.5
REF (MAXMIN)
1.5
2.40 0.21
1.25 0.111
1.111
13.8 2.8
REF (KM++)
12.1 + 2.7
2.49 + 0.25
1.28 0.10
18
13.6 = 2.9
REF (GKM)))+-
13.8 1 2.2
3.06 0.23
2.06 0.12
2.
19.1 + 3.0
PNNS(UNIP)
11.2 2.0
2.79 0.09
1.78
3.7
14.2
PXNS(MAXMIX)
8.4 1.7
2.51 0.10
1.60 0.09
18
2.5
(KNS(KM++))
9.1 1.8
2.30 0.07
1.39 0.07
1.2
13.1 2.8
PNNS PNNS(GRM +
12.1 2.0
0.08
2
2.5 2.4 1.4
19.9 = 3.3
PNNS UNIF; 10)
30.0 1.9
6.60 0.16
5.59 1.111
5.1120
50.9 1 2.5
PNNS MAXMIN; p 10)
26.1 3.5
0.12
0.12
3.7 3.7 0.8
54.9 2.0
PNNS KM+-; 10)
1.00
6.26 0.06
5.31 0.07
0.7
55.3 -
PKNS(CKKM++: p= 10)
31.8 1.2 1.2
7.07 0.21
6.24
4.8 0.8
60.9 = 2.0
1 16
16
1
292
34
A.1.2 SSE cost
seeder
A3
BirchJ
Birch2
Unholance
Dim1024
UNIF
11.5
110
1.56 + 0.22
6.3
(1.4 + 0.4) 104
MAXMIN
8.5
106.5 + 3.2
0.84 + 0.10
6 4
5.3938 + 0.0
KM+
9.4 0.9
104.7 1 3.0
0.85 1 0.10
1.1 1.1::00 0.6
20 : 220
GKM+++
7.7 0.5
0.528 0.035
0.68 0.16
5.3938
REF(UNIF)
=
104.0 2.8
1.22 + 0.11
4.5 1.3
(5.3 2.4) 100
REF(MAXMIN)
7.17 1: 0.28
100.4 1: 1.9
0.56 1 0.05
0.646935 0.0
5.3938 1 0.0
RBF(KM++)
8.4 + 0.5
101.0 + 2.5
0.71 + 0.05
0.648 0.030
5.3938 + 0.0
REF(GKM
7.16 1 0.26
96.7 1. 1.7
0.481 1 0.021
0.646935 1. 0.0
5.3938 1 0.0
PNNS(UNIP)
6.90 + 0.27
92.77290 + 0.00006
0.497 + 0.030
1.0 0.6
(3 9) 102
PNNS(MAXMIN)
6.73784 0.00012
92.77290 + 0.00005
0.456724 0.0
0.646935 0.0
5.3938 + 0.0
PNNS(KM 1)
6.75 0.07
92.77290 + 0.00006
0.456724 0.0
0.646935 0.0
5.3938
+)
6.73782 0.00011
92.77290 1 0.00006
0.456724
0.0
0.646935
0.0
5.3938
PNNS(UNIF; p 10)
6.75 0.07
92.77290 0.00006
0.468
0.021
0.82 - 0.39
5.3938 0.0
PNNS(MAXMIN; p=
6.73786 0.00013
92.77290 + 0.00006
0.456724 0.0
0.646935 0.0
5.3938 + 0.0
PNNS(KM) +: 10)
6.73785 1 0.00012
92.77290 1 0.00006
0.456724 1. 0.0
0.646935 0.0
5.3938 1 0.0
PNNS/CKM++: 10)
6.73785 + 0.00012
92.77290 + 0.00006
0.456724 0.0
0.646935 + 0.0
5.3938 + 0.0
PNN
6.73774 0.0
92.7729
0.456724 0.0
0.646995
5.3938
16

---

Published in Transactions on Machine Learning Research (12/2022)
A.1.3 CI (average)
seeder
A3
Birch1
Bircha
Unbalance
Dim1024
6.6 = 1.6
6.7 + 1.7
2.5
3.93 + 0.35
3.7 3.7
MAXMIN
2.9 1 1.0
5.6 1 1.3
7.4 d 1.8
0.9 i 0.6
KM-++
4.1 = 1.2
4.8 + 1.2
7.1 = 1.6
0.5 = 0.6
0.004 + 0.063
1.7 1 0.8
2.7 1.1
1.6 - 0.8
0.05 0.23
REP(UNIF)
4.7
4.5
124418 1.8
3.4
1.4
REP(MANMIN)
0.8
2.4
RRP(KM++++
2.8
1.0
4.9
0.00
REP(CROC++)
0.8 0.3
1.6 0.7
0.6 0.3
0
PNNS(UND))
0.3 0.5
0.9 0.7
0.4 + 0.6
0.08 + 0.27
(SASSMAANMIN))
0
0
0
PNNS(KM+-)
0.02 1 0.14
0
O
P
PANS(UNII; p= 10)
PINS(MAXMIN; p 10;
PANS(KM+-= 10)
PNNS(GKM+++; 10)
0
0
PNN
0
0
A.1.4 Success rate
seeder
A8
Birch1
Birch2
Unbalance
Dima 1024
UNIP
0
0
0.001
MAXMIN
0.004
0
0
0.223
KM++
0.541
0.996
GKM+H-
0.051
0.02
0.06
0.946
(UNIP)
0.055
REF (MANMIN)
0.254
0.01
1
REF (KM++++
0.002
0.998
REF (GKM 1)
0.239
0.03
0.44
1
PNNS(UNIP)
0.709
0.26
0.696
0.919
(MAMMI)
1
PNNS(KM 1)
0.98
1
PNNS(CKM++++)
PNNS(UNIF; P 10)
0.982
0.76
0.806
PNNS/ MAXMIN; p 10)
PNNS(KM 10)
1
1
PNNS(GKM++H p= 10)
1
1
1
PNN
1
1
17

---

Published in Transactions on Machine Learning Research (12/2022)
A.1.5 Normalized distance computations
seoder
A3
Birch1
Birck2
Unbalance
Dim1024
UNIF
1.80 + 0.18
2.52 0.37
1.146 + 0.027
2.7 0.6
1.89 0.32
MAXMIN
1.33 + 0.09
2.08
L073 +
1.26 + 0.10
1.07812 0.0
KM++
1.41 0.12
2.00 = 0.25
1.063 + 0.014
1.35 + 0.31
1.079 = 0.009
GKM+++
5.08 0.05
6.16
5.972 + 0.004
3.78
3.89062 0.0
REF(UNIP)
3.94 + 0.12
3.75
2.393 + 0.019
4.5
6.6 = 0.4
REF (MANMIX)
3.139 1 0.038
3,08 10.11
2.207 1 0.008
2,468 0.001
4.1136 1 0.0015
REF (KM++)
3.39 # 0.07
0.14
2,240 # 0.012
2.61 + 0.08
4.1140 + 0.0023
REF (GKG)+)++
6.992 + 0.029
0.11
7.1128 = 0.0030
5.068 0.028
6.9261 0.0014
(NNS(UNIP)
3.84 = 0.06
0.04
3.191 + 0.012
4.63 0.37
4.27 0.16
PNNS(MAXMIN)
0.025
3.662 0.025
1 3.084 0.008
0.026
3.380 1 0.027
P)NS(KM+++)
3.439 + 0.039
3.791 + 0.031
3.40 0.05
3,380 = 0.027
PNNS(CKM+++)
7.137 + 0.022
8.416 + 0.020
7.957 + 0.007
5.903 0.034
6.193 = 0.027
PNNS(UNIE; p 10)
12.37 1 0.10
12.027 1 0.037
11.314 - 0.035
12.38 1 0.35
20.6 1 0.5
PNNS(MAXMIN; p= 10)
11.78 + 0.09
11.567 + 0.030
11.089 = 0.035
1.13
22.4 12.4:000 0.5
PNNS(KM 1: p= 10)
11.88 + 0.08
11.716 + 0.037
11.154 = 0.034
1.14
0.5
PNNS(CKM+++: p= 10)
15.63 + 0.08
16.385 = 0.037
16.015 = 0.032
13.98 0.13
25.2 0.5
262.324
1730.39
1738.25
1422.21
531.759
A.1.6 Lloyd's iterations
seader
A3
Birch1
Birch2
Unbalance
Dim1024
UNIF
26 1 10
120 1: 40
45 1 12
33 i 15
2.7 1 1.3
MAXMIN
10
95 =
41
3.2 1.1
1
KM
99
41
L
1.00 1.0.7
GKM++
13
60 23
23
3
1
REF(UNIP)
26 13
100 100+:00 40
56 - 18
32
4.6
REF (MANMIN)
141 10
68 1 34
35 19
2.3
2.001 1: 0.031
REF (KM++)
21
73 37
4
5#4
2.001 + 0.031
REF REF(GKMM +G
1
58 32
17
2.
2.000 + 0.019
(NNS(UNIE)
18
36
38 19
18 16
3.1 1.0
PNNS(MAXMIN)
10.3
29
1214
2.1
210
+)
13:4
31
1544
4.4
210
P)NS(CKM+++)
9.4
23
9.1
2.9
2
PNNS(UNIE; P 10)
3.2
2
18
13=14
3.0 0.4
PNNS(MAXMIN; p= 10)
2.2
20
2.8
2.2
2.00 - 0.04
PNNS(KM++H ++; p= 10)
8.8
22 22-4
2.5
4.0
2.01 + 0.08
PNNS(CKM++:: p= 10)
2.1
18.2 + 3.9
5.9
3.3
2.000 0.015
PNN
18

---

Published in Transactions on Machine Learning Research (12/2022)
A.2 Real-world datasets (non-parally
We tested 100 samples for each dataset (30 for Urban GB).
A.2.1 Convergence time (in s)
sooder
Bridge (x10-3)
House (x10-1)
Miss A. (x10-1)
Urb.GB
Olivetti (x10-2)
Isolct (x10-1)
UNIP
4.17 1: 0.38
0.6
1.13 0.09
1.63 0.23
2.7 0.6
1.52 1 0.30
MAXMIN
4.6 0.5
3.5
1.78 + 0.33
1.73 0.20
4.1 0.6
1.58 + 0.21
KM
4.22 1 0.36
2.40 1 0.37
1.15 0.11
1.56 0.13
4.10 1 0.39
1.75 1 0.29
GKM++
6.45 + 0.35
3.33 = 0.33
1.58 = 0.13
5.51
8.8 0.6
3.02 + 0.29
REF(UNIF)
22.0 +.0 0.5
6.4
3.26 0.11
3.65 3.65 0.28
8.0
1.68 0.24
REF (MAXMIN)
21.3 0.6
0.3
3.49
1.244.224
8.5
1.56
REF(KM+++
21.5
5.76 0.38
3.35
1.310.20
9.1
REF(CKK+++)
0.6
6.38
0.19
5.36
12.0
2.19 0.22
(NNS(UNIF)
5.90 0.31
4.74
1.52 0.09
4.12
6.2 1.0
1.63 + 0.16
(NNS(MAX
5.88 0.29
4.94 1 0.36
1.68 0.14
4.33 0.18
6.6 1.44
1.79 1 0.25
PNNS(KM+++
5.9
4.44 0.31
1.100
3.58 0.16
6.74 + 0.38
1.66 + 0.17
PNNS(GKM +)
5.22 1: 0.32
1.80
5.11
9.51 0.39
2.07 0.18
PNNS(UNIP; p= 10)
15.8 0.9
8.88
3.31 1 0.08
10.76
21.1
3.85 0.21
PNNS(MAXMIN; p= 10)
0.8
0.12
10.84
21.5
3.71
PNNS(KM+++3 p= 10)
1.
8.33 0.30
0.10
9.97 0.18
218400
3.87 4 0.20
PNNS(GKM) 1:p= 10)
0.9
9.4 0.4
3.47 0.11
12.04 0.15
25.0 1 0.7
4.20 1 0.17
PNN
36.2
93.68 + 0.20
9.356 = 0.033
819 +
30.5 1.8
297.2 = 2.8
A.2.2 SSE cost (average)
seeder
Bridge (x10")
House (x10%)
Miss A. (x10P)
Urb.GB (x103)
Olivetti (x10
Isolet (x109)
UNIP
1.178 + 0.009
10.11 0.13
6.07 + 0.05
6.8 1.1
1,296 + 0.027
1.196 = 0.010
MAXMIN
1.139 + 0.005
10.18 = 0.07
5.794 = 0.028
2.99 0.07
1.252 + 0.015
1.233 0.014
KM+++
1.154 + 0.007
2.60::0000
5.687 1 0.037
2.72 1 0.06
1.276 0.021
1.196 = 0.010
GKM
1.124 0.004
9.529 0.027
5.507 + 0.016
2.430 + 0.020
1,227 0.013
1.190 = 0.007
REF(UNIF)
1.162 0.006
9.95
5.850 0.035
5.3
0.023
1.189 0.006
REF(MAXMIN)
1.160 0.005
9.96 1 0.07
5.845 1 0.033
2.94 1 0.05
1.258 0.020
1.234 0.013
REF(KM++)
1.158 0.005
9.592 + 0.038
5.792 + 0.033
2.63 0.04
1.283 0.022
1.191 0.008
REF(GKM
1.159 1 0.005
9.562 1 0.027
5.826 1 0.031
2.420 1 0.022
1.255 0.020
1.1837 T 0.0036
PNNS(UNIF)
1.124 + 0.005
9.553 + 0.034
5.607 + 0.040
2.70 + 0.06
1.215 00014
1.1800 + 0.0019
PNNS(MAXMIN)
1.0947 + 0.0027
9.508 + 0.025
5.339 + 0.008
2.388 + 0.013
1.192 0.008
1.1844 + 0.0029
PNNS(KM
1.1076 + 0.0033
9.486 1 0.020
5.403 0.013
2.323 + 0.010
1.209 1.013
1.1795 1.
+)
1.0947 1 0.0029
9.476 1 0.020
5.342
2.297 1
1.1790 1 0.0015
PNNS(UNIE; p= 10)
1.1020 + 0.0032
9.526 + 0.030
5.447 0.022
2.604 + 0.032
1.181 0.008
1.1807 + 0.0020
PNNS(MAXMIN; p 10)
1.0861 + 0.0022
9.482 + 0.021
5.327 + 0.007
2.349 + 0.010
1.1609 0.0034
1.1849 + 0.0025
PNNS(KM) +:; 10)
1.0884 T 0.0026
9.486 1 0.023
5.336 1 0.008
2,311 T 0.008
1.174 1 0.006
1.1806 1 0.0020
PNNS(CKM++: 10)
1.0850 + 0.0020
9.474 + 0.022
5.325 + 0.008
2,2992 - 0.0033
1.1671 + 0.0032
1.1795 + 0.0018
PNN
1.08279
9.49701
5.31588
2,3153
1.16238
1.17602
19

---

Published in Transactions on Machine Learning Research (12/2022)
A.2.3
SSE cost (minimum)
seeder
Bridge (x10?)
House (x109)
Miss A. (x10%)
(x102)
Olivetti (x10%)
Isolet (x109)
UNIP
1.15723
9.87713
5.95536
5.42031
1.23747
1.1795
MAXMIN
1.12869
10.0159
5. 5.72134
2.90322
1.21735
1.20639
KM++
1.14219
9.52638
5.59949
2.60788
1.23035
1.18025
GKM++
1.11421
9.45992
5.47583
2.37869
1.19908
1.17709
REF(UNIP)
1.14562
9.77128
5.76604
4.39285
1.24007
1.17688
REF(MAXMIN)
1.14719
9.81462
5.7451
2.84408
1.22378
1.20994
+KM((
1.13947
9.50834
5,72392
2.5686
1.22014
1,17933
REF(CKM+++)
1.14344
9.49776
5.74112
2.37821
1.2116
1.17737
PNNS(UNIF)
1.11345
9.47461
5.51748
2.5921
1.17155
1.17685
PNNS(MAXMIN)
1.09926
9.45145
5.32153
2.368
1.17443
1.17977
PNNS(KM 1)
1.10041
9.44324
5.37209
2.30731
1.18579
1.17687
+)
1.08729
9.42885
5.31768
2.28822
1.17066
1.17695
PUNS(UNIF; 10)
1.09472
9.45075
5.39731
2.51072
1.16219
1.17689
PANS(MAXMIX; p 10)
1.07998
9.43216
5.30742
2.33416
1.15955
1.18024
PNNS(KM+ +;; p= 10)
1.08029
9.42166
5.31398
2,29434
1.16141
1.17694
PNNS(CKM+++: 10)
1.08122
9.41836
5.30621
2,29477
1.16036
1.17689
1.08279
9.49701
5.31588
2.3153
1.16238
1.17692
A.2.4 Normalized distance computations
seeder
Bridge
House
Miss A.
Urb.GB
Olivetti
Isole!
UNIP
8.4 + 0.5
5.7 0.9
13.5 + 0.9
1.60 + 0.23
1.79 + 0.14
3.35 = 0.36
MAXMIN
8.9 + 0.8
4.5 0.6
23 1 4
1.196 = 0.037
2.21 + 0.07
3.14 0.25
KM+
8.2
3.71 + 0.39
13.7 + 1.3
1.229 + 0.029
2.26 1 0.07
3.36 + 0.35
GKM+
0.5
0.32
19.5 4
8.168 + 0.033
6.04 0.05
6.92 + 0.33
REF(UNIF)
0.8
8.6
39.4
3.15 0.15
5.68 0.10
4.62 + 0.26
REF(MAXMIN)
0.8
7.4
42.2 1 1.9
2.509 1 0.029
6.40 0.08
4.48 1 0.21
REF(KM++)
0.7
4.38
40.5
2.62 0.08
6.53
4.75 + 0.25
REP(GKM I)
4..8
12.25 1 0.35
48.3
9.465 1 0.021
10.34.0.00
8.20 1 0.22
PNNS(UNIP)
10.81 0.37
6.80 0.36
1.96
4.03 + 0.09
4.25 0.10
5.06 - 0.15
PNNS(MAXMIN)
10.82 + 0.37
6.11 + 0.31
1
3.335 + 0.025
4.82 182:0000
5.27 + 0.24
PNNS(KMI 1)
1..4
6.05 1 0.29
16.7110
3.31 0.06
5.00
5.24 - 0.17
PNNS(GRM+
1.4
1 0.30
22.2
10.191 0.014
0.09
8.73 1 0.16
PNNS(UNIF; p 10)
22.3 0.4
15.20 + 0.34
32.2 0.9
12.21 + 0.09
14.77 0.22
16.26 0.33
PNNS(MAXMIN; p 10)
20.00 + 0.35
14.25 + 0.32
28.9 + 1.4
11.301 0.032
15.48 0.19
15.46 + 0.26
PNNS(KMI 10)
21.11 1. 0.37
14.35 1 0.29
31.1 1 1.2
11.298 1 0.026
15.64 0.22
16.38 L 0.30
PNNS(GKM++: 10)
26.18 + 0.32
19.98 + 0.33
35.4 + 1.2
18.139 + 0.021
19.38 0.20
19.52 + 0.25
PNN
51.8177
271.383
87.8464
1351
19.7431
739.388
20

---

Published in Transactions on Machine Learning Research (12/2022)
A.2.5 Lloyd's iterations
seeder
Bridge
House
Miss A.
Urb. CB
Obiectii
Isolct
UNIF
= 5
124 1: 29
31 7
120 + 28
8.1 2.6
38 + 14
MAXMIN
134
96 + 29
2.6
28 :
KM++
100 26
17
8.1 2.4
38 = 14
GKM+++
90 25
86 25
6.1 1.5
34 = 12
(RF(UNIR)
26 10
116 + 27
40 = 13
130 + 60
10
37 + 17
REF REF (MAXMIN)
27 11
140 40
53 18
100 1 60
914
30 1 15
REF REF (KO++))
25 = 10
104 30
44 15
110 50
9#4
36 = 17
RRF RRF(GRR (GKM++)
27 12
52 20
100 50
8
33 16
(UNTP
24 6
94
38 9
109 109
2.8
33 13
PNNS PNNS(MAXMI)) (MAXMIN)
114
140
2.2
40 1
(KM+++))+++
25
27
93 30
2.6
34 14
PNNS (CKM++++++
24 7
93
42 11
86 28
5.6 2.3
30 + 14
(UNS, p= 10)
18 15 5
24
32 L1 7
103 1 30
4.2 1.5
32 1 13
PNNNS PNNS(MAXMING MAXMIN; p= 10)
17
93
36 = 10
111 + 38
3.5
32 - 14
PNNS(KM p: 10)
17 5
85
34
86 1 34
3.9
29 1 12
PNNS p= 10)
16
9
74 24
3.7
26 - 11
21

---

Published in Transactions on Machine Learning Research (12/2022)
A.3 datasets (parallel, 4 threads)
We tested 30 samples for each dataset.
A.3.1 Convergence time (in s)
seeder
Gr
USCensus
UNIF
0.86 + 0.13
13.2 # 3.3
MAXMIN
1.06 + 0.12
19.8 + 2.1
KM++
0.91 + 0.13
17.8 + 2.2
GKM
4.71 1 0.17
42.5
REF(UNIF)
1.74 0.12
17.1
RRF (RF(MIXMIN)
1.59 0.19
23.7
(EF(KM] 1)
1.51
22.3 2.2
REF(GKM++)
2.34
44.2 1: 2.0
PNNS(UNIF))
2.02 2.02::.111
16.0 = 1.0
(NNN(MAXMIN)
2.11 0.17
20.1 + 1.5
PNNS(KM
0.10
18.9
(KKS(C)M-+))
2.45 0.13
34.8 +
(NNS(UNIF;; P 10)
6.29 0.288
50.2
PNNS(MAXMIN; p= 10)
6.32
53.1
(NNS(KM++; 10)
6.04 + 0.24
51.6
PNNS(GKM 1: p= 10)
6.38 1: 0.18
59.5 1.7
A.3.2 SSE cost (average)
seeder
Urb. GB (x10?)
USCensus (x107)
UNIP
6.8 + 1.1
1.200 + 0.010
MAXMIN
2.99 + 0.07
1.32 - 0.04
KM++
2.72 0.05
1.192 0.011
GKM++
2.430 + 0.020
1.172 0.007
REF(UNIF)
5.26 - 0.35
1.181 +: 0.008
REF(MANMIN)
2.93 0.05
1.300 0.032
(EF(KM] 1)
2.64 1 0.05
1.180 1 0.006
REF RGF(GKM++)
2.413 + 0.017
1.156 + 0.005
PNNS(UNIP)
2.70 + 0.05
1.1423 + 0.0014
(MNS(MAXMIN)
2.382 1 0.013
1.1760 0.0030
PNNS(KM+++)
2.325 0.008
1.1434 0.0017
+)NS(CKM+++)
2.301 + 0.005
1.1405 0.0007
PNNS(UNIF; p 10)
2.595 + 0.031
1.1422 + 0.0014
PAKS(MAXMIN; 10)
2.349 1 0.011
1.1744 1 0.0031
PNKS(KM++; p=
2,311 + 0.006
1.1432 + 0.0016
PNNS(GKM +1: p= 10)
2,297 +: 0.004
1.1402 1 0.0006
22

---

Published in Transactions on Machine Learning Research (12/2022)
A.3.3 SSE cost (minimum)
seeder
Urb. GB (x109)
USCensus (x107)
UNIF
5,42031
1.18059
MAXMIN
2,90322
1.25624
KM+++
2.61092
1.17627
GKM+-
2.37869
1.15802
REF(UNIF)
4.65076
1.16404
(MAXMIX)
2.83593
1.2423
REF (KM++)
2.50803
1.16902
REP REF (GKG))++)
2.37949
1.14828
PNNS(UNIF)
2.61665
1.14017
(NNS(MAXMIN)
2.36183
1.16818
F-NS(KM+++)
2.31092
1.14155
PNNS(CKM++))
2.2858
1.13901
PNNS(UNIF; p 10)
2.52545
1.1392
PANS(MAXMIN; 10)
2.32724
1.16717
(KNS(KM++; p=
2,3022
1.14102
PNNS(CKM+++: P= 10)
2.2881
1.13892
23


Reviews:
[{"review_id": "e5IKlABhbG", "reviewer": "Reviewer_eetN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength\n\nThe proposed method is simple and easy to understand.\n\n\nWeakness\n\n1. This paper propose a strategy or trick (rather than an algorithm) to improve $k$-means initialization. Insides, it employs existing clustering algorithms, such as kmeans++, and existing merging method PNN. So the strategy is more like a combination and the contribution is limitted.\n\n2. The used strategy is easy and straight and can be easily devoloped with a low threshold. \n\n3. The clustering methods are not evaluated properly in experiment. There only one performance metric, i.e. success rate. Also success rates are not reported on real-world data (the results on synthetic data cannot be used in evaluation. They are only a showcase of your motivations etc.). In addition, the definition of success rate is not introduced properly.\n\n4. The proposed strategy shows worse results on multiple metrics, such as convergence time and normalized distance computations.\n\n5. The authors are expected to compare the complexiy of all considered methods with the proposed one. Only comparing the runtime is not enough, since there are so many factors can affect the runtime, such as the code quality, the programming language, running evironment, etc.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a meta-method for initializing the k-means clustering algorithm called PNN-smoothing. It consists in splitting a given dataset into J random subsets, clustering each of them individually, and merging the clustering results with the pairwise-nearest-neighbor (PNN) method.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Please see *Strengths And Weaknesses", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "HmabmU1qOA", "reviewer": "Reviewer_5asV", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: The idea is simple and the authors have performed sufficient numerical experiments.\n\nWeaknesses: \n\n1. According to the numerical results, the advantages of the method proposed in this work over the similar \"refine\" smoothing scheme seem to be minor for most cases (except for the case of convergence time for real-world datasets).\n\n2. There is no theoretical justification for the proposed method and some mathematical notions are vague. For example, \n* What is the precise definition of \"the minimum SSE cost\" (first appears on Page 10, Table 3)? \n\n* For Lloyd's algorithm, the set of centroids $\\mathcal{C}$ and the partition of data points are mutually dependent. I suggest that in (1), the sum-of-squared errors (SSE) should be defined only with respect to $\\mathcal{C}$ (or $\\mathcal{P}$), instead of both $\\mathcal{C}$ and $\\mathcal{P}$. \n\n* The sentence \"We thus define the number of normalized distance computations (NDC) as the total number of distance computations between $D$-dimensional vectors, divided by $Nk$\" (Page 2) is confusing to me.  I hope that the authors can provide a precise mathematical definition. \n\n3. The whole procedure (Algorithm 1) is used as an initialization method for the $k$-means clustering, but in the subsets, Lloyd's iterations are performed (which is the whole clustering procedure, rather than just initialization). This is a bit weird for me, especially when the number of subsets $J$ is pretty small (e.g., according to Table 1, if setting $\\rho =1$, for the Bridge and Olivetti datasets, $J = 3$).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a meta-method for the initialization of the $k$-means clustering. The method first splits the dataset into $J$ random subsets and clusters them separately, then merges the $Jk$ clusters to $k$ clusters by applying the pairwise-nearest-neighbor method to the $Jk$ centroids. The proposed initialization method can be readily combined with any seeding algorithm (followed by Lloyd's algorithm), whose computational complexity is linear in $k$ and the data size $N$, that is performed on the $J$ subsets. Extensive numerical experiments on both synthetic and real datasets have been performed to demonstrate systematically superior results with respect to the similar \"refine\" smoothing scheme. ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "See the Weaknesses above. ", "broader_impact_concerns": "Not applicable. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6AKTwQHWyo", "reviewer": "Reviewer_CQMT", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n---\n\nExtensive experiments with many popular seeding methods are tested, and it is helpful for both practitioners and theorists to see how all of them compare. The proposal of NDC which is a hardware-independent measure for benchmarking the performance of clustering algorithms. This NDC measure is quite natural and indeed captures the “unit” of computation for clustering algorithms, which may be applied more widely and have a lasting impact.\n\nWeakness:\n---\n\nWhile the effort of conducting the experiment and benchmarking is appreciated, I don’t see the significant claimed advantage of the PNN method. To me, it seems PNN performs very much similar to existing methods; while it achieves marginally better performance, it also uses more computations. Hence, it is not convincing that PNN has a clear advantage.\n\nAnother weakness is the lack of discussion/exploration of some potentially interesting experimental findings. For instance, it is claimed that \\rho = 1 gives the best tradeoff - why is this? This is somewhat counterintuitive. Moreover, this paper does not really suggest new/improved methods, and hence lacks technical depth. For instance, regarding this PNN which is a meta method, is it useful to nest it, i.e., running another PNN as the basic clustering algorithm for each of the J parts?\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the performance of PNN-based (PNN = pairwise nearest neighbor) init seeding method for k-means clustering. The PNN method is a meta method, in the sense that it needs to be combined with an existing seeding algorithm, and enhances its performance. In particular, the PNN method first partitions the entire dataset into J random parts (where J is a parameter), and then performs an existing method on each part to obtain kJ centers, and eventually greedily and repeatedly finds and merges the two parts whose merging leads to the smallest overall cost. \n\nThe focus of this paper is to compare the empirical performance of PNN, to many other popular methods, including the basic/non-meta methods such as uniform-sampling (called UNIF), MAXMIN, KM++, and another meta method, the refine method (called REF). In addition to comparing the running time and accuracy, the paper also proposes a hardware independent measure for the running time of clustering algorithms, the normalized distance computations (NDS), denoting the number of distance computations between two data points.\n\nExtensive experiments have been done on both synthetic datasets and real datasets to evaluate the accuracy, convergence time for the follow-up Lloyd steps, and the NDC. For some larger datasets, parallel implementations are also provided and compared. Overall, the PNN methods work well compared with other baselines, and it almost always outperforms the corresponding basic seeding method without PNN, even though it requires a slightly more running time.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Critical to securing my recommendation for acceptance:\n----\n\n1. In the beginning of section 3, you said the experiments have \\rho = 1. Can you justify why this is the best choice? Perhaps you will need additional experiments with different \\rho’s, and also give a (theoretical) argument/explanation why this should happen. Moreover, it would be interesting to explore how this parameter should be set, for different types of datasets and init seeding methods?\n\n2. You pick different accelerator methods for every dataset. Does the choice of accelerator potentially the relative merit of PNN over other methods? If so, shouldn’t one use the same accelerator (or none), for fair comparison? \n\nChanges that strengthen your work:\n----\n\n1. Footnote in page 4: “triangular inequality” -> “triangle inequality”\n2. It is suggested to see how much percentage of time is spent on NDC, and how much is on other parts.\n3. In the test of synthetic data, it is suggested to also list the actual costs that the algorithms achieve, since finding the solution optimally seems a goal that is too restrictive.\n4. How do you implement the parallel algorithms? What part is parallelized?", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a meta-algorithm for improving the initialization of k-means clustering. The idea is to split the data into subsets, cluster each subset, and then merge the resulting clusters using a pairwise-nearest neighbor algorithm.

### Reviewer recommendations

From the reviewers' final recommendations, two were (leaning) positive, while another remained negative. The main critiques raised by the reviewers were as follows.

(1) _Limited novelty_. Reviewers found the core ideas in this paper a modest departure from prior work, with the final solution combining existing ideas. From my reading, there appears to be some merit to this view. However, the TMLR guidelines suggest we do not index on technical novelty, provided the results are sufficiently convincing and interesting.

(2) _Lack of theoretical results_. Reviewers noted that the paper is purely empirical in nature, with no theoretical results. The authors agree with this in the response, noting that most prior works do not have such guarantees. I agree with the critique, and certainly some theory would strengthen the paper, but do not find it disqualifying, as empirical papers could still be interesting to the community.

(3) _Lack of clarity_. Some reviewers raised concerns about clarity of presentation, including of some mathematical notions. I agree that there is some scope for improved clarity, particularly in terms of take-home messages (see below).

(4) _Mixed empirical results_. All reviewers had questions about whether the proposed method shows enough empirical promise over existing methods. The author response noted that their method is the only one that yields 100% accuracy on synthetic tasks, while also yielding consistent improvements over the INIT baseline on real-world tasks. I tend to agree with the authors that the synthetic results are a useful sanity check that the proposed method helps avoid "clear and easily-fixed mistakes". On the real-world datasets, the authors argue that an increase in compute time over INIT might be acceptable if there is a consistent improvement in quality; and that consistent improvements over INIT are non-trivial, e.g., not being achieved by the REF method. I tend to agree with the authors' argument, though the consistent reviewer confusion on this point suggests scope to be much clearer in the discussion of results.

To summarize, while critiques (1)-(4) are reasonable, (1) and (2) do not strongly influence my decision per the TMLR guidelines. While I understand the arguments behind (4), on balance I tend to view the results as being of potential interest.


### AE recommendation

The TMLR guidelines for acceptance are based on the following two questions:

(a) _Are the claims made in the submission supported by accurate, convincing and clear evidence?_

(b) _Would some individuals in TMLR's audience be interested in the findings of this paper?_

For (a), the major potential issue is whether the empirical results indeed convey the potential benefits of the PNNS method. As above, I believe they do reasonably show the potential benefits of the method in settings where consistently improved SSE cost, with bounded increase in compute cost, is desirable.

For (b), two of the reviewers believe that readers may find this paper interesting. From my reading, I agree with this assessment.

Given these, we recommend acceptance with minor revisions as suggested below.


### Suggested revisions

- The title is a bit verbose. Perhaps "Systematically and efficiently improving k-means initialization by pairwise-nearest-neighbor smoothing" conveys the same effect.

- Section 2: it would be useful to have a table summarizing the complexity (+ possibly NDCs) of all methods listed.

- Section 2: for each method, the summary is provided in a single para. This ends up quite long for the last two methods, making it hard to parse. It would be better to split the paras, e.g., so that the discussion of computational complexity appears on a separate para.

- When referring to the refine method, e.g., "Our method starts out identically to refine", it might be useful to instead use "REF(INIT)", or use a different typeface for refine.

- Since all reviewers had questions about the significance of the empirical results, it would be useful to explicitly and succinctly state (e.g., above Section 4.1) what the goals of the experiments are: to highlight that the proposed method helps avoid "clear and easily-fixed mistakes" in synthetic settings, and offers consistently improved SSE cost with bounded increase in compute cost in real-world settings.

- Further to the above, the points on experiments raised in the "Revised version and General response to reviewers" should be incorporated in the discussion of the experimental results.

Decision: Accept (poster)