Improving Scalable K-Means++

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Improving Scalable K-Means++ © 2020 by the authors. Licensee MDPI, Basel, Switzerland Published version Hämäläinen, Joonas; Kärkkäinen, Tommi; Rossi, Tuomo Hämäläinen, J., Kärkkäinen, T., & Rossi, T. (2021). Improving Scalable K-Means++. Algorithms, 14(1), Article 6. https://doi.org/10.3390/a14010006 2021 algorithms Article Improving Scalable K-Means++ Joonas Hämäläinen * , Tommi Kärkkäinen and Tuomo Rossi   Citation: Hämäläinen, J.; Kärkkäinen, T.; Rossi, T. Improving Scalable K-Means++. Algorithms 2021,14, 6. https://doi.org/10.3390/a14010006 Received: 25 November 2020 Accepted: 21 December 2020 Published: 27 December 2020 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This articleisan openaccess article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license(https://creativecommons.org/ licenses/by/4.0/). Faculty of Information Technology, University of Jyväskylä, 40014 Jyväskylä, Finland; [email protected] (T.K.); [email protected] (T.R.) *Correspondence: [email protected] Abstract: Two new initialization methods for K-means clustering are proposed. Both proposals are based on applying a divide-and-conquer approach for the K-means k type of an initialization strategy. The second proposal also uses multiple lower-dimensional subspaces produced by the random projection method for the initialization. The proposed methods are scalable and can be run in parallel, which make them suitable for initializing large-scale problems. In the experiments, comparison of the proposed methods to the K-means++ and K-means k methods is conducted using an extensive set of reference and synthetic large-scale datasets. Concerning the latter, a novel highdimensional clustering data generation algorithm is given. The experiments show that the proposed methods compare favorably to the state-of-the-art by improving clustering accuracy and the speed of convergence. We also observe that the currently most popular K-means++ initialization behaves like the random one in the very high-dimensional cases. Keywords: clustering initialization; K-meansk; K-means++; random projection 1. Introduction Clustering is one of the core techniques in data mining. Its purpose is to form groups from data in a way that the observations within one group, the cluster, are similar to each other and dissimilar to observations in other groups. Prototype-based clustering algorithms, such as the popular K-means [ 1 ], are known to be sensitive to initialization [ 2 , 3 ], i.e., the selection of initial prototypes. A proper set of initial prototypes can improve the clustering result and decrease the number of iterations needed for the convergence of an algorithm [ 3 , 4 ]. The initialization of K-means was remarkably improved by the work of Arthur and Vassilvitskii [ 5 ], where they proposed the K-means++ method. There, the initial prototypes are determined by favoring distinct prototypes, which in high probability are not similar to the already selected ones. A drawback of K-means++ is that the initialization phase requires K inherently sequential passes over the data, since the selection of a new initial prototype depends on the previously selected prototypes. Bahmani et al. [ 6 ] proposed a parallel initialization method called K-means k (Scalable K-means++). The K-means k speeds up initialization by sampling each point independently and by updating sampling probabilities less frequently. Independent sampling of the points enables parallelization of the initialization, thus providing a speedup over K-means++. However, for example MapReduce-based implementation of K-means k needs multiple MapReduce jobs for the initialization. The MapReduce K-means++ method [ 7 ] tries to address this issue, as it uses one MapReduce job to select K initial prototypes, which speeds up the initialization compared to K-means k . Suggestions of parallelizing the second, search phase of K-means have been given in several papers (see, e.g., [ 8 , 9 ]). On a single machine, distance pruning approaches can be used to speed up K-means without affecting the clustering results [ 10 – 14 ]. Besides parallelization and distance pruning, data summarization is also a viable option for speeding up the K-means clustering [15,16]. Algorithms 2021,14, 6. https://doi.org/10.3390/a14010006 https://www.mdpi.com/journal/algorithms Algorithms 2021,14, 6 2 of 20 Dimension reduction has had an important role in making clustering algorithms more efficient. Over the years, various dimension reduction methods have been applied to decrease the dimension of data in order to speed up clustering algorithms [ 17 – 20 ]. The key idea for improved efficiency is to solve an approximate solution to the clustering problem in a lower-dimensional space. Dimension reduction methods are usually divided into two categories: feature selection methods and feature extraction methods [ 21 ]. Feature selection methods aim to select a subset of the most relevant variables from the original variables. Correspondingly, feature extraction methods aim to transform the original dataset into a lower-dimensional space while trying to preserve the characteristics (especially distances between the observations and the overall variability) of the original data. A particular dimensional reduction approach for processing large datasets is the random projection (RP) method [ 22 ]. Projecting data from the original space to a lowerdimensional space while preserving the distances is the main characteristic of the RP method. This makes RP very appealing in clustering, whose core concept is dissimilarity. Moreover, classical dimension reduction methods such as the principal component analysis (PCA) [ 23 ] become expensive to compute for high-dimensional spaces whereas RP remains computationally efficient [24]. Fern and Brodley [ 18 ] proposed an ensemble clustering method based on RP. They showed empirically that aggregation of clustering results from multiple lower-dimensional spaces produced by RP leads to better clustering results compared to a single clustering in lower-dimensional space produced by PCA or RP. Other combinations of K-means and RP have been studied in several papers [ 17 , 25 – 27 ]. RP for K-means++ was analyzed in [ 28 ]. Generally, the main idea is to create a lower-dimensional dataset with RP and to solve the ensuing K-means clustering problem with less computational effort. On the other hand, one can also optimize clustering method’s proximity measure for small datasets [29]. In general, K-means clustering procedure typically uses a non-deterministic initialization, such as K-means++, followed by the Lloyd’s iterations [ 1 ] —with multiple restarts. Prototypes corresponding to the smallest sum-of-squares clustering error are selected as the final clustering result. In [ 30 ], such a multistart strategy was carried out during the initialization phase, thus reducing the need to repeat the whole clustering algorithm. More precisely, a parallel method based on K-means++ clustering of subsets produced by the distribution optimally balanced stratified cross-validation (DOB-SCV) algorithm [ 31 ] was proposed and tested. Here, such an approach is developed further with the help of K-means k and RP. More precisely, we run K-means k method in a low-dimensional subset created by RP. In contrast to the previous work [ 30 ], the new methods also restrict the number of Lloyd’s iterations in the subsets. Vattani [ 32 ] showed by construction that the number of iterations, and thus the running time, of the randomly initialized K-means algorithm can grow exponentially already in small-dimensional spaces. As stated in the original papers [ 5 , 6 ], the K-means++ and Kmeans k readily provide improvements to this both in theory and in practice. Concerning our work, we have provided time complexity analysis for SK-means k in Section 3.1 and for SRPK-means k in Section 3.2. In terms of time complexity, SRPK-means || reduces the time complexity of the initialization for large-scale high-dimensional data (this was also confirmed by our experimental results) and provides better clustering results; thus, it reduces need for restarts compared to the baseline methods. Reduced need for restarts also improves the overall time complexity of K-means algorithm. In terms of clustering accuracy, SK-means does this same effect for large-scale lower-dimensional datasets. Moreover, we showed for the synthetic datasets (M-spheres) that the random projection variant of the initialization (SRPK-means k ) can provide clear advantage in very high-dimensional cases, where the distances can become meaningless for other distance-based initialization methods. The main purpose of this article is to propose two new algorithms for clustering initialization and compare them experimentally to the initializations of K-means++ and K-means k using several large-scale datasets. To summarize the main justification of the proposed methods: they provide better results compared to baseline methods with bet- Algorithms 2021,14, 6 3 of 20 ter or equal running time. The proposed initialization method reduces data processing with sampling, subsampling, and dimensional reduction solving the K-means clustering problem in a coarse fashion. Moreover, from the perspective of parallel computing, using a parallelizable clustering method in the subset clustering allows fixing the number of subsets and treating each subset locally in parallel, hence improving the scalability. For quantified testing and comparison of the methods, we introduce a novel clustering problem generator for high dimension spaces (see Section 3.4). Currently, challenging simulated datasets for high-dimensional clustering problems are difficult to find. For instance, the experiments with DIM datasets of tens or hundreds of dimensions in [ 4 ] were inconclusive: all clustering results and cluster validation index comparisons behaved perfectly without any errors. Therefore, better experimental datasets are needed and can be produced with the proposed algorithm. 2. Existing Algorithms In this section, we introduce the basic composition of the existing algorithms. 2.1. K-Means Clustering Problem and the Basic Algorithms Let X={x1 , x2 ,..., xN} be a dataset such that xi∈RM∀ 1 ≤i≤N , where M denotes the dimension, and let C={c1 , c2 ,..., cK} be a set of prototypes, where each prototype also belongs to RM . The goal of the K-means clustering algorithm is to find a partition of X into Kdisjoint subsets, by minimizing the sum-of-squares error (SSE) defined as SSE(C) = ∑ x∈X min c∈Ckc−xk2. (1) An approximate solution to the minimization problem with (1) is typically computed by using the Lloyd’s K-means algorithm [ 1 ]. Its popularity is based on simplicity and scalability. Even if the cost function in (1) is mathematically nondifferentiable because of the min -operator, it is easy to show that after the initialization, the K-means type of iterative relocation algorithm converges in finite many steps [4]. Prototype-based clustering algorithms, such as K-means, are initialized before the prototype relocation (search) phase. The classical initialization algorithm, readily proposed in [ 33 ], is to randomly generate the initial set of prototypes. A slight refinement of this strategy is to select, instead of random points (from appropriate value ranges), random indices and use the corresponding observations in data as initialization [ 34 ]. Because of this choice, there cannot be empty clusters in the first iteration. Bradley and Fayyad [ 35 ] proposed an initialization method where J randomly selected subsets of the data are first clustered with K-means. Next, it forms a superset of the J×K prototypes obtained from the subset clustering. Finally, the initial prototypes are achieved as the result of K-means clustering of the superset. Arthur and Vassilvitskii [ 5 ] introduced the K-means++ algorithm, which improves the initialization of K-means clustering. The algorithm selects first prototype at random, and then the remaining K− 1 prototypes are sampled using probabilities based on the squared distances to the already selected set, thus favoring distant prototypes. The generalized form of such an algorithm with different lp -distance functions and the corresponding cluster location estimates was depicted in [4]. The parallelized K-means++ method, called K-means k , was proposed by Bahmani et al. [ 6 ] (see Algorithm 1). In Algorithm 1, and from here onwards, symbol “#” denotes ’number of’. In the algorithm, sampling from X is conducted in a slightly different fashion compared to K-means++. More precisely, the sampling probabilities are multiplied with the over-sampling factor l and the sampling is done independently for each data point. The initial SSE for the first sampled point ψ determines the number of sampling iterations. K-means k runs O(log(ψ)) sampling iterations. For each iteration, the expected number of points is l . Hence, after O(log(ψ)) iterations, the expected number of points added to C is O(l log(ψ)) . Finally, weights representing the accumula- Algorithms 2021,14, 6 4 of 20 tion of data around the sampled points are set and the result of the weighted clustering then provides the K initial prototypes. K-means++ can be used to cluster the weighted data (see Algorithm 1 in [ 36 ]). Selecting r= 5 instead of O(log(ψ)) rounds and setting the over-sampling factor to 2 K were demonstrated to be sufficient in [ 6 ]. Recently, Bachem et al. [ 36 ] proved theoretically that small r instead of O(log(ψ)) iterations is sufficient in K-means k . A modification of K-means k for initializing robust clustering was described and tested in [37]. Algorithm 1: K-meansk Input: Dataset X, #clusters K, and over-sampling factor l. Output: Set of prototypes C={c1,c2, ..., cK}. 1: C←select point c1uniformly random from X. 2: ψ←compute SSE(C). 3: for O(log(ψ)) times do 4: C0←sample each point x∈Xindependently with probability l·d(x)2/SSE(C). 5: C←C∪C0 6: For each xin Cattach a weight defined as the number of points in Xcloser to xthan any other point in C. 7: Do a weighted clustering of Cinto Kclusters. 2.2. Random Projection The background for RP [ 22 ] comes from the Johnson–Lindenstrauss lemma [ 38 ]. The lemma states that points in a high-dimensional space can be projected to a lower dimension space while approximately preserving the distances of the points, when the projection is done with a matrix whose elements are randomly generated. Hence, for an N×M dataset X , let R∈M×P be a random matrix. Then, the random projected data matrix e X is given by e X=1 √PXR . The random matrix R consists of independent random elements (rij) which can be drawn from one of the following probability distributions [ 22 ]: rij = + 1 with probability 1 / 2, or − 1 with probability 1 / 2; or rij = + 1 with probability 1 / 6, 0 with probability 2/3, or −1 with probability 1/6. 3. New Algorithms Next we introduce the novel initialization algorithms for K-means, their parallel implementations, and the novel dataset generator algorithm. 3.1. SK-Meansk The first new initialization method for K-means clustering, Subset K-means k (SKmeans k ), is described in Algorithm 2. The method is based on S randomly sampled nondisjoint subsets {X1 , X2 ,..., XS} from X of approximately equal size, such as X=∪S i=1Xi . First, K-means k is applied in each subset, which gives the corresponding set of initial prototypes Ci . Next, each initial prototype set Ci in Xi is refined with Tinit Lloyd’s iterations. Tinit is assumed to be significantly smaller than the number of Lloyd’s iterations needed for convergence. Then, SSE is computed locally for each Ci in Xi . Differently from the earlier work [ 30 ], this locally computed SSE is now used as the selection criteria for the initial prototypes instead of the global SSE. Computation of SSE for Xi in Step 3 is obviously much faster than to compute it for the whole X . However, a drawback is that if the subsets are too small to characterize the whole data, the selection of the initial prototypes might fail. Algorithms 2021,14, 6 5 of 20 Therefore, S should be selected such that the subsets are sufficiently large. For example, if S is close to the number of samples in the smallest cluster, then this cluster will appear as an anomaly for most of the subsets. On the other hand, this property can also be beneficial to exclude anomalous clusters already in the initialization phase. Currently, there are no systematic comparisons in the literature on the size of the subsets in sampling-based clustering approaches [ 39 ]. In [ 35 ], the number of subsets was set to 10. Based on [ 30 , 35 ], this selection appears reasonable for the sampling-based clustering initialization approaches. Algorithm 2: SK-meansk Input: Subsets {X1,X2, ..., XS}, #clusters K, and #Lloyd’s iterations Tinit. Output: Set of prototypes C={c1,c2, ..., cK}. 1: Ci←for each subset Xirun K-meansk. 2: Ci←for each subset Xirun Tinit Lloyd’s iterations initialized with Ci. 3: Compute local SSE for each Ciin Xi. 4: C←select prototypes corresponding to smallest local SSE. The convergence rate of K-means is fast and the most significant improvements in the clustering error are achieved during the first few iterations [ 40 , 41 ]. Therefore, for the initialization purposes, Tinit can restricted, e.g., to 5 iterations. Moreover, since the number of Lloyd’s iterations needed for convergence might vary significantly (e.g., [ 4 ]), a restriction on the number of Lloyd’s iterations helps in synchronization, when a parallel implementation of the SK-meanskmethod is used. The computational complexity of the K-means k method is of the order O(rlNM) , where r is the number of initialization rounds. Therefore, SK-means k also has the complexity of the order O(rlNM) in Step 1. In addition, SK-means k runs Tinit Lloyd’s iterations with the complexity of O(TinitKNM) , and computes local SSE with the complexity of O(KNM) . Hence, the total complexity of SK-means k is of the order O(rlNM +TinitKNM) . 3.2. SRPK-Meansk The second novel proposal, Subset Random Projection K-means k (SRPK-means k ), adds RPs to SK-means k . Since SK-means k mainly uses time in computing distances in Steps 1 and 2, it is reasonable to speed up the distance computation with RP. The RP-based method is presented in Algorithm 3. Generally, SRPK-means k computes a set of candidate initial prototypes in a lower-dimensional space and then evaluates these in the original space. As with Algorithm 2, the best set of prototypes based on the local SSE are selected. Algorithm 3: SRPK-meansk Input: Subsets {X1,X2, ..., XS}, #clusters K, #Lloyd’s iterations Tinit, and random projection dimension P. Output: Set of prototypes C={c1,c2, ..., cK}. 1: Ri←for each subset Xigenerate M×Prandom matrix. 2: e Xi←for each subset Xicompute 1 √PXiRi 3: e Ci←for each e Xirun K-meansk. 4: Ii←for each e Xirun Tinit Lloyd’s iterations initialized with e Ci. 5: For each partitioning Iicompute prototypes Ciin original space Xi. 6: Compute local SSE for each Ciin Xi. 7: C←select prototypes corresponding to smallest local SSE. Algorithms 2021,14, 6 6 of 20 The proposal first computes a unique random matrix for each subset Xi . Then, the P -dimensional random projected subset e Xi is computed in each subset Xi . Steps 3–4 are otherwise the same as the Steps 1–2 in Algorithm 2, but these steps are applied for the lower-dimensional subsets { e X1 , e X2 ,..., e XS} . Next, the labels Ii for partitioning each subset are used to compute Ci in the original space Xi . Finally, the local SSEs are computed, and the best set of prototypes are returned as the initial prototypes. Please note that the last two steps in Algorithm 3are the same as Steps 3–4 in Algorithm 2. SRPK-means k computes projected data matrices, which require a complexity of O(PNM) (naive multiplication) [ 17 ]. Execution of K-means k in the lower-dimensional space requires O(rlNP) , and Tinit Lloyd’s iterations requires O(TinitKNP) operations. Step 6 requires O(KNM) operations, since it computes the local SSEs in the original space, so that the total computational complexity of the SRPK-means k method is O(PNM +rlNP +TinitKNP+KNM) . Typically, applications of RP are based on the assumption P<< M . Thus, when the dimension of data M is increased, the contribution of the second and the third term of the total computational complexity start to diminish. Moreover, when both M and K are large compared to P , the last term dominates the overall computational complexity. Therefore, in terms of running time, SRPK-means k is especially suited for clustering large-scale data with very high dimensionality into a large number of clusters. Fern and Brodley [ 18 ] noted that clustering with RP produces highly unstable and diverse clustering results. However, this can be exploited in clustering to find different candidate structures of data, which then can be combined into a single result [ 18 ]. The proposed initialization method in this paper uses a similar idea as it tries to find structures from multiple lower-dimensional spaces that minimize the local SSE. In addition, selecting a result that gives the smallest local SSE excludes the bad structures, which could be caused by inappropriate Rior Ci. 3.3. Parallel Implementation of the Proposed K-Means Initialization Algorithms Bahmani et al. [ 6 ] implemented K-means k with the MapReduce programming model. It can also be implemented by the Single Program Multiple Data (SPMD) programming model with message passing. Then all the steps of the parallelized Algorithms 1–3are executed inside an SPMD block. Next, a parallel implementation of K-means k as depicted in Algorithm 1is briefly described, by using Matlab Parallel Computing Toolbox (PCT), SPMD blocks, and message passing functions (see [ 42 ] for a detailed description about PCT). First, data X is split into Q subsets of approximately equal size and then the subsets are distributed to Q workers. Step 1 picks a random point from a random worker and broadcasts this point to all other workers. In Step 2, each worker calculates distances and SSE for its local data. Next, points are aggregated by calling gplus-function, after which the aggregation distributes this sum to other workers. In Steps 4 and 5, each worker samples points from its local data, the next points are aggregated to C0 by calling gop-function, and then C0 is broadcasted to all workers. Again, distances and SSE are calculated similarly as in Step 2. Each worker in Step 6 assigns weights based on its local data, after which the weights are aggregated with gop-function. Finally, Step 7 is computed sequentially. As with the parallel K-means k implementation, a parallel implementation of Algorithm 2with SMPD and message passing is described next. First, each subset Xi from S subsets is split into J approximately equal size subsets and then these subsets are distributed to J×S workers, e.g., subset Xi is distributed to workers (i− 1 )J+ 1,..., (i− 1 )J+J . In Steps 1–3, each subset of workers runs steps for subset Xi in parallel similarly as described in the previous paragraph. For parallel Lloyd’s iterations, a similar strategy as proposed in [ 8 ] can be used in Step 2. Steps 1–3 require calling modified gop-function and gplus-function for the subset of workers; these functions were modified to support this requirement. Finally, prototypes corresponding to the smallest local SSE from the subset i0 allocated workers (i0−1)J+1, ..., (i0−1)J+Jare returned as the initialization. The parallel SRPK-means k in Algorithm 3can be implemented in a highly similar fashion to the parallel SK-means k . More precisely, in Step 1, each worker (i− 1 )J+ 1, Algorithms 2021,14, 6 7 of 20 where i∈ { 1,2,..., S} , generates the random matrix Ri and broadcasts it to workers (i− 1 )J+ 1,..., (i− 1 )J+J . In Step 2, each worker computes random projected data for its local data. Steps 3–4 are otherwise computed similarly to the parallel SK-means k Steps 1 and 2, except these steps are executed for the projected subsets. In Step 5, the prototypes are computed in the original space in parallel. Finally, Steps 6 and 7 are the same as Steps 3 and 4 in Algorithm 2. The parallel implementations of the proposed methods and K-means k are available in (https://github.com/jookriha/Scalable-K-means). 3.4. M-Spheres Dataset Generator The novel dataset generator is given in Algorithm 4. Based on the method proposed in ([ 43 ], p. 586), Algorithm 5generates a random point that is on the M -dimensional sphere centered on c with radius of d . The core principle is to draw M independent values from the standard normal distribution and transform these with corresponding M -direction cosines. The obtained M -dimensional vector is then scaled with the radius d and relocated with the center c . The generated points are uniformly distributed on the surface of the sphere, because of the known properties of the standard normal distribution (see [ 43 ] (p. 587) and articles therein). Algorithm 4: M-spheres dataset generator Input: #clusters K, #dimensions M, #points per cluster NK, nearest center distance dc, radius of M-sphere dr. Output: Dataset X={x1,x2, ..., xN}. 1: C←{(0,...,0)}. 2: if K>1then 3: c2←randsurfpoint(c1,dc). 4: C←C∪ {c2}. 5: k←2. 6: if K>2then 7: while k<Kdo 8: i←rand({1, 2, ..., k}). 9: ccand ←randsurfpoint(ci,dc).. 10: i∗←argminjkcj−ccandk. 11: if i∗== ithen 12: C←C∪ {ccand}. 13: k←k+1. 14: k←1. 15: X← {}. 16: while k≤Kdo 17: n←1 18: while n≤NKdo 19: d∗ r←rand((0, dr]) 20: xnew ←randsurfpoint(ck,d∗ r). 21: X←X∪ {xnew}. 22: n←n+1 23: k←k+1 The generator uses Algorithm 5for generating K cluster centers so that kci−cjk=dc , where i6=j , dc is the given distance between the centers, and both ci , cj then belong to Algorithms 2021,14, 6 8 of 20 the set of centers C . Finally, NK data points for each cluster are generated by applying Algorithm 5with a uniformly random radius from the interval ( 0, dr] . This means, in particular, that points in a cluster are nonuniformly distributed and approximately 100a dr % percentages of the points are within a M -dimensional sphere with radius a for 0 ≤a≤dr . In Algorithm 5,N(0,1)denotes the standard normal distribution. Algorithm 5: randsurfpoint(c,d) Input: Sphere center c, sphere radius d. Output: New point x∗. 1: m←1. 2: S←0. 3: while m≤Mdo 4: xm←N(0,1). 5: S←S+x2 m. 6: m←m+1. 7: x∗←c+dS−1/2x For simplicity, generation of the cluster centers starts from the origin. When the new centers are randomly located with the fixed distance and then expanded as clusterwise data in RM , the generator algorithm does not restrict the actual values of the generated data and centers. Hence, depending on the input, data range can be large. However, this can be alleviated with the min-max scaling as part of the clustering process. 4. Empirical Evaluation of Proposed Algorithms In this section, empirical comparison between K-means++, K-means k , SK-means k , and SRPK-means k is presented by using 21 datasets. In Section 4.1, the results are given for 15 reference datasets. The performance of the methods was evaluated by analyzing SSE, the number of iterations needed for convergence, and the running time. Finally, In Section 4.2, we analyze the final clustering accuracy for six novel synthetic datasets that highlight the effects of the curse of dimensionality in the K-means++ type initialization strategies. The simulated clustering problems have been formed with the novel generator described in Algorithm 4. The MATLAB implementation of the algorithm is available in (https://github.com/jookriha/M_Spheres_Dataset_Generator). 4.1. Experiments with Reference Datasets In this section, the results are shown and analyzed for 15 publicly available reference datasets by considering separately the accuracy (Section 4.1.2), efficiency (Section 4.1.3), and scalability (Section 4.1.4) of the algorithms. 4.1.1. Experimental Setup Basic information about the datasets is shown in Table 1. The parallel implementations of the proposed methods and K-means k (omitting K-means++ readily tested in [ 6 ]) were applied to the seven largest datasets and serial implementations were used otherwise. For the serial experiments, we used the following eight datasets: Human Activity Recognition Using Smartphones (http://archive.ics.uci.edu/ml/index.php) (HAR), ISOLET (ISO), Letter Recognition (LET), Grammatical Facial Expressions (GFE), MNIST (http: //yann.lecun.com/exdb/mnist/) (MNI), Birch3 (http://cs.joensuu.fi/sipu/datasets/) (BIR), Buzz in Social Media (BSM), and Covertype (COV). For the parallel experiments, the following seven large high-dimensional datasets were used: KDD Cup 1999 Data (KDD), US Census Data 1990 (USC), Oxford Buildings (http://www.robots.ox.ac.uk/~vgg/data/ Algorithms 2021,14, 6 15 of 20 where I refers to the cluster labels, L to the ground truth labels, K∗ to the number of unique ground truth labels, and n denotes labels’ frequency counts. Then, NMI can be defined as NMI(I,L) = 2MI(I,L) H(I) + H(L), (3) where H(I) = log2(K) and H(L) = log2(K∗) denote the entropies of the cluster labels and ground truth labels. Details of the six generated datasets with Algorithm 4are summarized in Table 5. The generated datasets are referred as M-spheres. For each dataset, we set NK= 10,000, K= 10, and dr= 1. To demonstrate interesting effects in the clustering initialization, we varied the nearest cluster center distance as dc={ 0.05,0.1,0.2 } and the data dimension as M={ 1000,10,000 } . We set dc values to much smaller than dr in order to increase the difficulty of the clustering problems. In Figure 2, PCA projections on the three largest principal components show that the clusters are more separated for M= 10,000 than for M= 1000. For the M-spheres datasets, we used the serial implementations of the initialization methods with the same settings as before. Table 5. Characteristics of the synthetic datasets. Dataset #Observations (N) #Features (M) #Clusters (K) Center Distance (dc) Radius (dr) M-spheres-M1k-dc0.05 100,000 1000 10 0.05 1.0 M-spheres-M1k-dc0.1 100,000 1000 10 0.1 1.0 M-spheres-M1k-dc0.2 100,000 1000 10 0.2 1.0 M-spheres-M10k-dc0.05 100,000 10,000 10 0.05 1.0 M-spheres-M10k-dc0.1 100,000 10,000 10 0.1 1.0 M-spheres-M10k-dc0.2 100,000 10, 000 10 0.2 1.0 Results for the final clustering accuracy using NMI with 100 repeats are shown in Figure 3. Clearly, SRPK-means k outperforms other methods in terms final clustering accuracy for all the synthetic datasets. Moreover, if we compare the results between the datasets with M= 1000 and M= 10,000, we observe that the clustering accuracy for SRPKmeans k is improved when the dimensionality increases. The most significant difference is obtained for the M-spheres-M10kdc 0.05 dataset, where K-means++ has a total breakdown of the accuracy while SRPK-means k is able find the near optimal clustering result out of 100 repeats. Moreover, the accuracy of K-means++ is clearly worse compared to K-means k and SK-means k for very high-dimensional datasets. We tested K-means with random initialization for this dataset and observed from the statistical testing that K-means++, K-means k and SK-means k are no better than the random initialization in terms of NMI. These results demonstrate that the use of distances in the K-means++ type of initialization strategies can become meaningless in very high-dimensional spaces. In Figure 4, scatter plots of SSE and NMI values show that the relative SSE difference between worst possible clustering result ( NMI = 0) and the optimal clustering result ( NMI = 1) can be surprisingly small for very high-dimensional data. Therefore, the improvements for the final clustering accuracy in Table 2can be much more significant than the impression given by SSE in terms of how spherical clusters are found for highdimensional datasets. Figures 3and 4illustrate the deteriorating behavior of the currently most popular K-means++ initialization method in high dimensions. We especially observe that the K-means++ initialization behaves like (i.e., is not better than) the random one in the very high-dimensional cases. Such finding also suggests further experiments, where as a function of the data dimension, emergence of such a behavior is being studied to identify most appropriate random project dimensions to restore the quality of initialization and the whole clustering algorithm. Algorithms 2021,14, 6 16 of 20 (a) M-spheres-M1k-dc0.05 (b) M-spheres-M1k-dc0.1 (c) M-spheres-M1k-dc0.2 (d) M-spheres-M10k-dc0.05 (e) M-spheres-M10k-dc0.1 (f) M-spheres-M10k-dc0.2 Figure 2. Synthetic dataset projections on the three largest principal components. K-means++ K-means|| SK-means|| SRPK-means|| P = 5 SRPK-means|| P = 10 SRPK-means|| P = 20 SRPK-means|| P = 40 0 0.2 0.4 0.6 0.8 1 NMI (a) M-spheres-M1k-dc0.05 K-means++ K-means|| SK-means|| SRPK-means|| P = 5 SRPK-means|| P = 10 SRPK-means|| P = 20 SRPK-means|| P = 40 0 0.2 0.4 0.6 0.8 1 NMI (b) M-spheres-M10k-dc0.05 Figure 3. Cont. Algorithms 2021,14, 6 17 of 20 K-means++ K-means|| SK-means|| SRPK-means|| P = 5 SRPK-means|| P = 10 SRPK-means|| P = 20 SRPK-means|| P = 40 0 0.2 0.4 0.6 0.8 1 NMI (c) M-spheres-M1k-dc0.1 K-means++ K-means|| SK-means|| SRPK-means|| P = 5 SRPK-means|| P = 10 SRPK-means|| P = 20 SRPK-means|| P = 40 0 0.2 0.4 0.6 0.8 1 NMI (d) M-spheres-M10k-dc0.1 K-means++ K-means|| SK-means|| SRPK-means|| P = 5 SRPK-means|| P = 10 SRPK-means|| P = 20 SRPK-means|| P = 40 0 0.2 0.4 0.6 0.8 1 NMI (e) M-spheres-M1k-dc0.2 K-means++ K-means|| SK-means|| SRPK-means|| P = 5 SRPK-means|| P = 10 SRPK-means|| P = 20 SRPK-means|| P = 40 0 0.2 0.4 0.6 0.8 1 NMI (f) M-spheres-M10k-dc0.2 Figure 3. NMI for the synthetic datasets. 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (a) Random 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (b) K-means++ 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (c) K-meansk 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (d) SK-meansk 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (e) SRPK-meanskp=5 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (f) SRPK-meanskP=10 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (g) SRPK-meanskp=20 2.42 2.43 2.44 SSE 107 0 0.5 1 NMI (h) SRPK-meanskp=40 Figure 4. Scatter plot of SSE and NMI results for the M-spheres-M10k-dc0.05 dataset. Algorithms 2021,14, 6 18 of 20 5. Conclusions In this paper, we proposed two parallel initialization methods for large-scale K-means clustering and a new high-dimensional clustering data generation algorithm to support their empirical evaluation. The methods are based on divide-and-conquer type of K-means k approach and random projections. The proposed initialization methods are scalable and fairly easy to implement with different parallel programming models. The experimental results for an extensive set of benchmark and novel synthetic datasets showed that the proposed methods improve clustering accuracy and the speed of convergence compared to state-of-the-art approaches. Moreover, the deteriorating behavior of the K-means++ and K-means k initialization methods in high dimensions can be recovered with the proposed RP-based approach to provide accurate initialization also for high-dimensional data. Our experiments also confirmed the finding (e.g., [ 52 ]) that the difference between the errors (SSE) of good and bad clustering results in high-dimensional spaces can be surprisingly small also challenge cluster validation and cluster validation indices (see [4] and references therein) in such cases. Experiments with SRPK-means k method demonstrate that use of RP and K-means k is beneficial for clustering large-scale high-dimensional datasets. In particular, SRPK-means k is an appealing approach as a standalone algorithm for clustering very high-dimensional large-scale datasets. In future work, it would be interesting to test a RP-based local SSE selection for SRPK-means k , which uses the same RP matrix in each subset for the initial prototype selection. In this case, use of sparse RP variants [ 54 ] or the mailman algorithm [ 17 ] for the matrix multiplication could be beneficial, particularly in applications where K is close to P . Furthermore, integrating the proposed methods into the robust Kmeans k [ 37 ] would be beneficial for clustering noisy data, because the clustering problems in these cases are especially challenging. Author Contributions: Conceptualization, J.H., T.K. and T.R.; Data curation, J.H.; Formal analysis, J.H.; Funding acquisition, T.K.; Investigation, J.H.; Methodology, J.H.; Project administration, T.K.; Resources, J.H.; Software, J.H.; Supervision, T.K. and T.R.; Validation, J.H., T.K. and T.R.; Visualization, J.H.; Writing—original draft, J.H.; Writing—review & editing, J.H. and T.K. All authors have read and agreed to the published version of the manuscript. Funding: The work has been supported by the Academy of Finland from the projects 311877 (Demo) and 315550 (HNP-AI). Conflicts of Interest: The authors declare no conflict of interest. References 1. Lloyd, S. Least squares quantization in PCM. IEEE Trans. Inf. Theory 1982,28, 129–137. [CrossRef] 2. Emre Celebi, M.; Kingravi, H.A.; Vela, P.A. A comparative study of efficient initialization methods for the k-means clustering algorithm. Expert Syst. Appl. 2012,40, 200–210. [CrossRef] 3. Fränti, P.; Sieranoja, S. How much can k-means be improved by using better initialization and repeats? Pattern Recognit. 2019 , 93, 95–112. [CrossRef] 4. Hämäläinen, J.; Jauhiainen, S.; Kärkkäinen, T. Comparison of Internal Clustering Validation Indices for Prototype-Based Clustering. Algorithms 2017,10, 105. [CrossRef] 5. Arthur, D.; Vassilvitskii, S. k-means++: The advantages of careful seeding. In Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, Society for Industrial and Applied Mathematics, New Orleans, LA, USA, 7–9 January 2007; pp. 1027–1035. 6. Bahmani, B.; Moseley, B.; Vattani, A.; Kumar, R.; Vassilvitskii, S. Scalable K-means++. Proc. VLDB Endow. 2012 ,5, 622–633. [CrossRef] 7. Xu, Y.; Qu, W.; Li, Z.; Min, G.; Li, K.; Liu, Z. Efficient k -Means++ Approximation with MapReduce. IEEE Trans. Paral. Distrib. Syst. 2014,25, 3135–3144. [CrossRef] 8. Dhillon, I.S.; Modha, D.S. A data-clustering algorithm on distributed memory multiprocessors. In Large-Scale Parallel Data Mining; Springer: New York, NY, USA, 2002; pp. 245–260. 9. Zhao, W.; Ma, H.; He, Q. Parallel K-Means Clustering Based on MapReduce. In Proceedings of the 1st International Conference on Cloud Computing, CloudCom ’09, Munich, Germany, 19–21 October 2009; pp. 674–679. Algorithms 2021,14, 6 19 of 20 10. Elkan, C. Using the triangle inequality to accelerate k-means. In Proceedings of the 20th International Conference on Machine Learning (ICML-03), Washington, DC, USA, 21–24 August 2003; pp. 147–153. 11. Hamerly, G. Making k-means even faster. In Proceedings of the 2010 SIAM International Conference on Data Mining, Columbus, OH, USA, 29 April–1 May 2010; pp. 130–140. 12. Drake, J.; Hamerly, G. Accelerated k-means with adaptive distance bounds. In Proceedings of the 5th NIPS Workshop On Optimization for Machine Learning, Lake Tahoe, NV, USA, 7–8 December 2012; pp. 42–53. 13. Ding, Y.; Zhao, Y.; Shen, X.; Musuvathi, M.; Mytkowicz, T. Yinyang k-means: A drop-in replacement of the classic k-means with consistent speedup. Int. Conf. Mach. Learn. 2015,37, 579–587. 14. Bottesch, T.; Bühler, T.; Kächele, M. Speeding up k-means by approximating Euclidean distances via block vectors. Int. Conf. Mach. Learn. 2016, 2578–2586. Available online: http://proceedings.mlr.press/v48/bottesch16.pdf (accessed on 24 December 2020) 15. Bachem, O.; Lucic, M.; Lattanzi, S. One-shot coresets: The case of k-clustering. In International Conference On Artificial Intelligence And Statistics; PMLR: Canary Islands, Spain, 2018; pp. 784–792. 16. Capó, M.; Pérez, A.; Lozano, J.A. An efficient K-means clustering algorithm for massive data. arXiv 2018, arXiv:1801.02949. 17. Boutsidis, C.; Zouzias, A.; Drineas, P. Random projections for k -means clustering. Adv. Neural Inf. Process. Syst. 2010 ,23, 298–306. 18. Fern, X.; Brodley, C. Random projection for high dimensional data clustering: A cluster ensemble approach. In Proceedings of the 20th International Conference on Machine Learning, Washington, DC, USA, 21–24 August 2003; pp. 186–193. 19. Alzate, C.; Suykens, J.A. Multiway spectral clustering with out-of-sample extensions through weighted kernel PCA. IEEE Trans. Pattern Anal. Mach. Intell. 2010,32, 335–347. [CrossRef] 20. Napoleon, D.; Pavalakodi, S. A new method for dimensionality reduction using K-means clustering algorithm for high dimensional data set. Int. J. Comput. Appl. 2011,13, 41–46. 21. Liu, H.; Motoda, H. Feature Selection for Knowledge Discovery and Data Mining; Kluwer Academic Publishers: Norwell, MA, USA, 1998. 22. Achlioptas, D. Database-friendly random projections: Johnson-Lindenstrauss with binary coins. J. Comput. Syst. Sci. 2003 , 66, 671–687. Special Issue on {PODS} 2001. [CrossRef] 23. Jolliffe, I.T. Principal Component Analysis; Springer: New York, NY, USA, 2002. 24. Bingham, E.; Mannila, H. Random Projection in Dimensionality Reduction: Applications to Image and Text Data. In Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, 26–29 August 2001; pp. 245–250. 25. Cohen, M.B.; Elder, S.; Musco, C.; Musco, C.; Persu, M. Dimensionality reduction for k-means clustering and low rank approximation. In Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing, Portland, OR, USA, 15–17 June 2015; pp. 163–172. 26. Boutsidis, C.; Zouzias, A.; Mahoney, M.W.; Drineas, P. Randomized dimensionality reduction for k -means clustering. IEEE Trans. Inf. Theory 2014,61, 1045–1062. [CrossRef] 27. Cardoso, Â.; Wichert, A. Iterative random projections for high-dimensional data clustering. Pattern Recognit. Lett. 2012 , 33, 1749–1755. [CrossRef] 28. Chan, J.Y.; Leung, A.P. Efficient k-means++ with random projection. In Proceedings of the 2017 International Joint Conference on Neural Networks (IJCNN), IEEE, Anchorage, AK, USA, 14–19 May 2017; pp. 94–100. 29. Sarkar, S.; Ghosh, A.K. On perfect clustering of high dimension, low sample size data. IEEE Trans. Pattern Anal. Mach. Intell. 2019. Available online: https://www.isical.ac.in/~statmath/report/74969-cluster.pdf (accessed on 24 December 2020). 30. Hämäläinen, J.; Kärkkäinen, T. Initialization of Big Data Clustering using Distributionally Balanced Folding. In Proceedings of the European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning-ESANN 2016, Bruges, Belgium, 27–29 April 2016; pp. 587–592. 31. Moreno-Torres, J.G.; Sáez, J.A.; Herrera, F. Study on the Impact of Partition-Induced Dataset Shift on k-fold Cross-Validation. IEEE Trans. Neural Netw. Learn. Syst. 2012,23, 1304–1312. [CrossRef] [PubMed] 32. Vattani, A. K-means requires exponentially many iterations even in the plane. Discr. Comput. Geom. 2011 ,45, 596–616. [CrossRef] 33. MacQueen, J. Some methods for classification and analysis of multivariate observations. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability; University of California Press: Berkeley, CA, USA, 1967; Volume 1, pp. 281–297. Available online: https://stacks.stanford.edu/file/druid:xb208zr6261/xb208zr6261.pdf (accessed on 24 December 2020). 34. Forgy, E.W. Cluster analysis of multivariate data: Efficiency vs. interpretability of classifications. Biometrics 1965,21, 768–769. 35. Bradley, P.S.; Fayyad, U.M. Refining Initial Points for K-Means Clustering. ICML 1998,98, 91–99. 36. Bachem, O.; Lucic, M.; Krause, A. Distributed and provably good seedings for k-means in constant rounds. Int. Conf. Mach. Learn. 2017,70, 292–300. 37. Hämäläinen, J.; Kärkkäinen, T.; Rossi, T. Scalable robust clustering method for large and sparse data. In Proceedings of the European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning-ESANN 2018, Bruges, Belgium, 25–27 April 2018; pp. 449–454. 38. Johnson, W.B.; Lindenstrauss, J. Extensions of Lipschitz mappings into a Hilbert space. Contemp. Math. 1984,26, 1. 39. Rezaei, M.; Fränti, P. Can the Number of Clusters Be Determined by External Indices? IEEE Access 2020 ,8, 89239–89257. [CrossRef] Algorithms 2021,14, 6 20 of 20 40. Bottou, L.; Bengio, Y. Convergence properties of the k-means algorithms. Adv. Neural Inf. Process. Syst. 1995 , 585–592. Available online: https://papers.nips.cc/paper/1994/file/a1140a3d0df1c81e24ae954d935e8926-Paper.pdf (accessed on 24 December 2020). 41. Broder, A.; Garcia-Pueyo, L.; Josifovski, V.; Vassilvitskii, S.; Venkatesan, S. Scalable k-means by ranked retrieval. In Proceedings of the 7th ACM International Conference on Web Search and Data Mining, New York City, NY, USA, 24–28 February 2014; pp. 233–242. 42. Sharma, G.; Martin, J. MATLAB®: A Language for Parallel Computing. Int. J. Parallel Progr. 2009,37, 3–36. [CrossRef] 43. Muller, M.E. Some continuous Monte Carlo methods for the Dirichlet problem. Ann. Math. Stat. 1956,27, 569–589. [CrossRef] 44. Fränti, P.; Sieranoja, S. K-means properties on six clustering benchmark datasets. Appl. Intell. 2018,48, 1–17. [CrossRef] 45. De Vries, C.M.; Geva, S. K-tree: Large scale document clustering. In Proceedings of the 32nd International ACM SIGIR Conference on Research and Development in Information Retrieval, Boston, MA, USA, 19–23 July 2009; pp. 718–719. 46. Kriegel, H.P.; Schubert, E.; Zimek, A. The (black) art of runtime evaluation: Are we comparing algorithms or implementations? Knowl. Inf. Syst. 2017,52, 341–378. [CrossRef] 47. Gallego, A.J.; Calvo-Zaragoza, J.; Valero-Mas, J.J.; Rico-Juan, J.R. Clustering-based k-nearest neighbor classification for large-scale data with neural codes representation. Pattern Recognit. 2018,74, 531–543. [CrossRef] 48. Kruskal, W.H.; Wallis, W.A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 1952,47, 583–621. [CrossRef] 49. Saarela, M.; Hämäläinen, J.; Kärkkäinen, T. Feature Ranking of Large, Robust, and Weighted Clustering Result. In Proceedings of 21st Pacific Asia Conference on Knowledge Discovery and Data Mining-PAKDD 2017, Jeju, Korea, 23–26 May 2017; pp. 96–109. 50. Äyrämö, S. Knowledge Mining Using Robust Clustering; Vol. 63 of Jyväskylä Studies in Computing; University of Jyväskylä: Jyväskylä, Finland, 2006. 51. Amdahl, G.M. Validity of the single processor approach to achieving large scale computing capabilities. AFIPS Conf. Proc. 1967 , 30, 483–485. 52. Aggarwal, C.C.; Hinneburg, A.; Keim, D.A. On the surprising behavior of distance metrics in high dimensional space. In International Conference On Database Theory; Springer: New York, NY, USA, 2001; pp. 420–434. 53. Strehl, A. Relationship-Based Clustering and Cluster Ensembles for High-Dimensional Data Mining. Ph.D. Thesis, The University of Texas at Austin, Austin, TX, USA, 2002. 54. Li, P.; Hastie, T.J.; Church, K.W. Very sparse random projections. In Proceedings of the 12th ACM SIGKDD International Conference on Knowledge diScovery and Data Mining, Philadelphia, PA, USA, 20–23 August 2006; pp. 287–296.