# A Robustk-Means Clustering Algorithm Based on Observation Point Mechanism

Authors: Xiaoliang Zhang; Yulin He; Yi Jin; Honglian Qin; Muhammad Azhar; Joshua Zhexue Huang
Journal: Complexity
 (2020)

---

## Abstract

Thek-means algorithm is sensitive to the outliers. In this paper, we propose a robust two-stage k-means clustering algorithm based on the observation point mechanism, which can accurately discover the cluster centers without the disturbance of outliers. In the first stage, a small subset of the original data set is selected based on a set of nondegenerate observation points. The subset is a good representation of the original data set because it only contains all those points that have a higher density of the original data set and does not include the outliers. In the second stage, we use the k-means clustering algorithm to cluster the selected subset and find the proper cluster centers as the true cluster centers of the original data set. Based on these cluster centers, the rest data points of the original data set are assigned to the clusters whose centers are the closest to the data points. The theoretical analysis and experimental results show that the proposed clustering algorithm has the lower computational complexity and better robustness in comparison with k-means clustering algorithm, thus demonstrating the feasibility and effectiveness of our proposed clustering algorithm.


---

## Body

## 1. Introduction

Clustering is an important research branch of data mining. Thek-means algorithm is one of the most popular clustering methods [1]. When performing k-means clustering, we usually use a local search to find the solution [2, 3], i.e., selecting k points μ1,μ2,…,μk as the initial cluster centers and then optimizing them by an iterative process to minimize the following objective function (see, for example, [4, 5]):(1)E=∑i=1k∑Xj∈CiXj−μi22,where Xj is the j-th data point belonging to the i-th cluster Ci. It is well known that the solution of equation (1) is affected by the initial values of μii=1,2,…,k.In order to chooseμi properly, the k-means++ algorithm [6] picks out a set of points as the initial center points whose distances between each other are as large as possible. However, this method for choosing the initial center points is sensitive to outliers [7–9]. Some methods use the subsets of the original data set to determine μi. For instance, the CLARA [10] and CLARANS [11] algorithms use PAM [12] to calculate the initial cluster centers from the random subsets of the original data set. The sampling-based methods weaken the sensitivity because the sampling process can discard some outliers in the original data set, but it cannot guarantee all outliers to be ignored in the sampling process. Therefore, the remaining outliers in subsets still affect the clustering results.The automatic clustering algorithms are attracting more and more attention from the academic community, e.g., the density-based spatial clustering of applications with noise (DBSCAN) algorithm [13–15], depth difference-based clustering algorithm [16], and Tanir’s method [17]. Recently, a new automatic clustering algorithm named I-nice was proposed in [18]. Inspired by the observation point mechanism of I-nice algorithm, we propose a two-stage k-means clustering algorithm in this paper to find the cluster centers from a subset of the original data set with all outliers removed. In the first stage, we select a small subset of original data set based on a set of nondegenerate observation points. The subset contains only all the higher density points of the original data set and does not have the outliers. Therefore, it is a good representation of the original data set for finding the proper cluster centers. In the second stage, we perform the k-means algorithm on the subset to obtain a set of cluster centers and then the other points in the original data set can be clustered accordingly.Selecting the subset in the first stage is based on a set ofd+1 nondegenerate observation points that are assigned to the data space ℝd, where d is the dimension of data points. For each observation point, we compute a set of distances between it and all data points in the original data set. The set of distances generates a distance distribution with respect to the observation point. From the distance distribution, we identify the dense areas and extract the subset of data points in the dense areas. Then, we take the intersection of all d+1 subsets of data points in all dense areas from those d+1 distance distributions. After refining this intersection subset of data points, we obtain a subset without outliers of the original data set. Therefore, it can be used to find the proper cluster centers. Finally, we conduct some convictive experiments to validate the effectiveness of our proposed algorithm and the experimental results demonstrate that our proposed algorithm is robust to outliers.The remainder of this paper is organized as follows. We describe the related mathematical principles of our algorithm in Section2. The details of two-stage k-means clustering algorithm and its pseudocode are presented in Section 3. In Section 4, we present a series of experiments to validate the feasibility of our proposed algorithm. Finally, we summarize the conclusions and future work in Section 5.


## 2. Mathematical Principles

Definition 1.
SupposeD=X1,X2,…,XN is a data set including N data points with d dimensions. Given an observation point O∈ℝd, we say D˜=dXi,O:i=1,2,…,N is a generated distance set of D with respect to the observation point O, where dX,Y denotes the Euclidean distance between X and Y.
Given a data setD and an observation point O, we have(2)dX,O−dY,O≤dX,Y,by the triangle inequality for every X,Y∈D. Hence, the distance between two data points in D is larger than the difference of their corresponding two distances in D˜. Therefore, for any positive number r and a point X∈D, the number of points in D with distances to X less than r is not greater than the number of points in D˜ whose distances to dX,O are less than r. In particular, if X is a proper cluster center in D, dX,O will be a data point in D˜ which has more points close to it. That is to say, if X is a dense point in D, it is also corresponding to a dense point in D˜.
Unfortunately, the converse is not true. Because two points inD˜ which have a small distance may correspond to two points in D that have a large distance, a proper cluster center of D˜ may not be corresponding to a proper cluster center of D. Hence, we can deduce that D˜ retains the partial clustering information of D. In order to obtain more clustering information of D, one possible way is to choose more observation points to generate more distance sets and then combines all those different clustering information together. This is the main idea behind our new algorithm. We provide the following two theorems to guarantee the correctness of the abovementioned statements.Definition 2.
Given a set ofd+1 points O0,O1,…,Od where Oi=ai1,ai2,…,aid for i=0,1,2,…,d, define the generating matrix A of O0,O1,…,Od as(3)A=a11−a01a12−a02⋯a1d−a0da21−a01a22−a02⋯a2d−a0d⋮⋮⋮ad1−a01ad2−a02⋯add−a0d.
If the determinant ofA is not equal to zero, we say A is nondegenerate and O0,O1,…,Od is a set of nondegenerate points.Theorem 1.
SupposeO0,O1,…,Od is a set of nondegenerate points and let X1,X2∈ℝd. If(4)dX1,Oi=dX2,Oi,for all i=0,1,…,d, then X1=X2.Proof.
AssumeX1=x11,x12,…,x1d and X2=x21,x22,…,x2d. For each i=1,2,…,d, we have(5)dX1,Oi2−dX1,O02=dX2,Oi2−dX2,O02,i.e.,(6)∑t=1dx1t−ait2−∑t=1dx1t−a0t2=∑t=1dx2t−ait2−∑t=1dx2t−a0t2.
By simplifying equation (6), we can obtain(7)∑t=1d2x1t−ait−a0ta0t−ait=∑t=1d2x2t−ait−a0ta0t−ait.
Thus, we have(8)∑t=1dait−a0tx1t−x2t=0,i=1,2,…,d.
Since the coefficient matrix of (8) is nondegenerate, we can get that(9)x1t−x2t=0,t=1,…,d.
Thus,X1=X2.Remark 1.
For the convenience of calculation, we can chooseO0=0,0,…,0,Oi=0,…,0i−1,1i,0i+1,…,0,i=1,2,…,d as the set of nondegenerate observation points. In fact, if X=x1,x2,…,xd, then for each i, it has(10)dX,Oi2=∑jxj−oij2=∑j≠ixj−02+xi−12=∑jxj2−2∗xi+1=dX,O02−2∗xi+1.
Thus, if we have obtained the distance betweenX and O0, then computing the square of the distance between X and Oii=1,2,…,d will convert to addition operation three times, which will decrease the time complexity in generating those distance sets.Remark 2.
If the number of observation points is less thand+1, Theorem 1 does not hold true. For example, if we choose(11)O0=0,0,O1=1,0,O2=0,1,X1=0.2,0.2,X2=0.6,0.6,then, it has(12)dX1,O1=dX2,O1,dX1,O2=dX2,O2,but X1≠X2. By Theorem 1, we can confirm that all different clustering points can be distinguished by choosing a set of nondegenerate points as the observation points. Thus, d+1 is the minimum number of the observation points to distinguish all the cluster centers of the original data set.Theorem 2.
SupposeO0=0,0,…,0,Oi=0,…,0i−1,1i,0i+1,…,0, i=1,2,…,d. Let X1,X2∈ℝd and set(13)dXj,Oi=dji,for i=0,1,…,d and j=1,2. If, for each i=0,1,…,d,(14)d1i−d2i<r,M=maxj=1,2,i=0,1,…,ddji,then,(15)dX1,X2≤2dMr.Proof.
Forj=1,2, we set Xj=xj1,xj2,…,xjd. We have(16)∑t=1dxjt2=dj02,j=1,2;∑t=1dxjt2−2xji+1=dji2,i=1,2,…,d,j=1,2.
Solving the system of equation (16) results in(17)xji=121+dj02−dji2,i=1,2,…,d,j=1,2.
Then, we can obtain(18)dX1,X22=∑t=1dx1t−x2t2=∑t=1d121+d102−d1t2−121+d202−d2t22=14∑t=1dd102−d202−d1t2−d2t22≤14∑t=1drd10+d20+d1t+d2t2≤4dM2r2,which yields dX1,X2≤2dMr.Remark 3.
If we normalize the original data setD, for example, we perform the min-max normalization on D, and we can deduce that M≤d.Remark 4.
SupposeA is a generated distance set of D with respect to the observation point O. We cannot confirm whether two elements in A which have a small difference are corresponding to two data points in D which also have a small distance. But by Theorem 2, if all the d+1 pairs of generated distances of X and Y have a small difference, then X must have a small distance to Y. This can be used to adjust the dense of the selected subset.Remark 5.
The observation point mechanism aims to transform the original multidimensional data points into one-dimensional distance points, which is different from the landmark point or representative point mechanisms. The landmark point [19] is the core of landmark-based spectral clustering (LSC) algorithm which generates some representative data points as the landmarks and represents the remaining data points as linear combinations of these landmarks. The representative points [20] are the subset of original data set and used in the ultrascalable spectral clustering (U-SPEC) algorithm to alleviate the huge computational burden of spectral clustering. The observation points are designed to enhance the robustness of k-means clustering, while the landmark points or representative points are used to speed up the spectral clustering.


## 3. The Proposed Two-Stagek-Means Clustering Algorithm

Given a data setD with N objects, we want to partition D into k clusters. The main idea of our two-stage k-means clustering algorithm is that we only need to deal with a small subset of D which has a similar clustering structure to D. In order to select a proper subset with the abovementioned property, we need to discard all outliers in D and retain a portion of those points that are close to the cluster centers.


### 3.1. Description of Algorithm

#### 3.1.1. Generatingd+1 Distance Sets in the First Stage

First of all, we conduct the normalization operation on the original data setD. Set D=X1,X2,…,XN, where Xj=xj1,xj2,…,xjd for j=1,2,…,N. Suppose(19)M=maxj=1,2,…,N;i=1,2,…,dxjiandm=minj=1,2,…,N;i=1,2,…,dxji.Then, we transformD into D˜ with Xj=xj1,xj2,…,xjd corresponding to X˜j=xj1−m/M−m,xj2−m/M−m,…,xjd−m/M−m. Obviously, the transformation on D is a composition of a translation transformation and a dilation transformation. The dilation factor 1/M−m is the same for each dimension; hence, the dilation transformation does not change the cluster structure. Because the translation transformation also does not change the cluster structure of a dataset, the cluster structure of D is totally the same as that of D˜. We also note that the value of every component of X˜j is in the interval 0,1.Let(20)O0=0,0,…,0,Oi=0,…,0i−1,1i,0i+1,…,0,i=1,2,…,d,be the set of observation points. Denote D˜j the generated distance set of D˜ with respect to the observation point Ojj=0,1,…,d, and we get d+1 sets D˜0,D˜1,…,D˜d. For each data point X˜, we actually have mapped it to a d+1-dimensional vector dX˜,O0,dX˜,O1,…,dX˜,Od. Theorem 1 shows that we can identify X˜ by the d+1-dimensional vector, and hence, it is reasonable to expect that the clustering structure about D˜ can be deduced by those d+1 distance sets.


#### 3.1.2. Selecting a Representative Subset ofD˜ in the First Stage

For eachD˜i, we can get a set Si consisting of all candidate higher density points of D˜i by using the grid-based clustering methods (e.g., [21]). For example, first, we arrange D˜i in the ascending order. Second, a fixed value δi is selected to be a quantile of diffD˜i. Third, for each s in D˜i, we counter the number of elements of D˜i in the interval s−δi,s+δi. Thus, we obtain a positive integer sequence, where each member indicates the relative size of the density of the corresponding element of D˜i. Finally, we select out those s in D˜i such that the corresponding integer number is either a local maximum or beyond a threshold.In the following experiments, we will setδi as two times the p-th percentile of set diffDi for some p, where Di is the rearrangement of D˜i in the ascending order and diffDi is the sequence of the first-order difference on Di. Denote N as the cardinality of D˜i. If N is small, we usually choose a smaller p; for example, p=75. If N is very large, we choose a bigger p; e.g., p=99. Otherwise, we can choose a proper p between them; e.g., p=90.Now, we have obtainedd+1 sets S0,S1,…,Sd with each one containing all the higher-density points of the corresponding distance set. By the triangle inequality, we have the following property:If there is an
i
∈
0,1
,
…
,
d
such that
d
i
=
d
X
,
O
i
is not in
S
i
, then X cannot be a higher-density point of
D
˜.According to this property, we can select a subsetS of D˜ whose distances to the i-th observation point are in Si for all i∈0,1,…,d.For each pointX∈D˜, we have mapped it to a d+1-dimensional vector. By Remark 4 of Theorem 2, all the d+1 pairs of the corresponding components of two points that belong to the same cluster will have a little difference between them. But it is possible that there are some data points which have some components that have the little difference with that of one cluster center and have the other components that have the little difference with that of another cluster center. In such case, few outliers may be missed by the above selection criterion. To discard those few outliers in S and decrease the number of elements of S, we need to refine S. We have the following criterion according to Remark 4 of Theorem 2.Suppose
X
1
has been selected. Given a data point
X
2
, if, for every
i
∈
0,1
,
…
,
d
,
d
1
i
=
d
X
1
,
O
i
and
d
2
i
=
d
X
2
,
O
i
have a small difference, then we discard
X
2.We denoteS=Y1,Y2,…,Ym and set(21)Sc=∅,Sa=∅.We also need a counter to indicate the density of each data point ofSa. Firstly, we let Y1∈Sa and make the indicative number of Y1 to be 1. We then sequentially choose the data points in S and dynamically construct Sc and Sa according to the following process. Suppose we choose Yi from S, then we compute the distance between Yi and each data point in Sa. If there are some distances less than a threshold value δ, we add 1 to each of the counter of data point that corresponds to these distances and then discard Yi. Meanwhile, if the counter number of a data point in Sa is bigger than another threshold value n, then we remove this data point from Sa and add it into Sc. But if each one of the point in Sa has distance to Yi bigger than δ, we continue to check whether there is a point in Sc that has distance to Yi less than δ; we will discard Yi if there is any and we will add Yi to Sc if not. Finally, we obtain a set Sc that closely represents the original data set and the size of it is smaller than that of the original data set. Furthermore, all outliers of original data set are not included in this selected subset.


#### 3.1.3. ClusteringSc and D˜ in the Second Stage

Since the selected set has discarded all outliers and has a smaller size than the original data set, the running time will decrease significantly when performing thek-means algorithm on the selected subset. Furthermore, because the subset closely represents the original data set, the cluster centers will also be suitable to be chosen as the cluster centers of the original data set. When we have identified the cluster centers, it is easy to cluster the whole data set. The pseudo-code of our proposed algorithm is presented in Algorithm 1.Algorithm 1: Two-stage k-means clustering algorithm.
Input:
The number of clusters:k;
The number of percentile:p;
The original data set:D=Xi=xi1,xi2,…,xid:i=1,2,…,N;
Method:
NormalizeD and generate D˜=X˜i:i=1,2,…,N;
for t=0 to ddo
SetD˜t=dit:i=1,2,…,N where dit=dX˜i,Ot;
GenerateDt by rearranging D˜t in ascending order and set δt be p-th percentile of diffDt;
Setmt=minDt and Mt=maxDt/δt+1∗δt;
Equally divide the intervalmt,Mt into intervals with length δt;
LetSt be the union of those intervals that contain local maximum number of elements in Dt;
end for
Select out a subsetS:Xi∈S if and only if dit∈St for all t=0,1,…,d;
RefineS and obtain the subsets Sc:
Perform thek-means algorithm on Sc;
Assign eachX˜ii=1,2,…,N to the nearest center of the obtained cluster of Sc;
Output:
The result of clustering.


### 3.2. Analysis of Computational Complexity

In this section, we analyze the computational complexity of the proposed algorithm. When running the classicalk-means algorithm, each iteration needs to compute the distances between each data point in the whole data and those new modified cluster centers, which has a time complexity of ONkd. In our algorithm, the time cost in the first stage mainly consists of four parts. The first part is to generate d+1 one-dimensional data sets, which has a time complexity of ONd. The second one is to find those intervals which contain the local maximum of distances, which has a time complexity of ON. The third part is to select S, which has the time complexity ONd. In the fourth part, we refine S and obtain the subset Sc, which has the time complexity ON˜2, where N˜ denotes the cardinality of S. Thus, the time complexity of the first stage is ONd+N˜2. At the second stage, we will perform the k-means algorithm on Sc and each iteration will have a time complexity less than N˜1kd, where N˜1 denotes the cardinality of Sc. Because N˜1<N˜, the total time complexity of the new algorithm is ONd+N˜kd+N˜2.We note that the time complexity of the fourth part in the first stage is usually much less thanON˜2. Since many data points have been discarded when constructing the sets Sc and Sa, we do not have to compute the distances with all N˜ data points in S.


## 3.1. Description of Algorithm

### 3.1.1. Generatingd+1 Distance Sets in the First Stage

First of all, we conduct the normalization operation on the original data setD. Set D=X1,X2,…,XN, where Xj=xj1,xj2,…,xjd for j=1,2,…,N. Suppose(19)M=maxj=1,2,…,N;i=1,2,…,dxjiandm=minj=1,2,…,N;i=1,2,…,dxji.Then, we transformD into D˜ with Xj=xj1,xj2,…,xjd corresponding to X˜j=xj1−m/M−m,xj2−m/M−m,…,xjd−m/M−m. Obviously, the transformation on D is a composition of a translation transformation and a dilation transformation. The dilation factor 1/M−m is the same for each dimension; hence, the dilation transformation does not change the cluster structure. Because the translation transformation also does not change the cluster structure of a dataset, the cluster structure of D is totally the same as that of D˜. We also note that the value of every component of X˜j is in the interval 0,1.Let(20)O0=0,0,…,0,Oi=0,…,0i−1,1i,0i+1,…,0,i=1,2,…,d,be the set of observation points. Denote D˜j the generated distance set of D˜ with respect to the observation point Ojj=0,1,…,d, and we get d+1 sets D˜0,D˜1,…,D˜d. For each data point X˜, we actually have mapped it to a d+1-dimensional vector dX˜,O0,dX˜,O1,…,dX˜,Od. Theorem 1 shows that we can identify X˜ by the d+1-dimensional vector, and hence, it is reasonable to expect that the clustering structure about D˜ can be deduced by those d+1 distance sets.


### 3.1.2. Selecting a Representative Subset ofD˜ in the First Stage

For eachD˜i, we can get a set Si consisting of all candidate higher density points of D˜i by using the grid-based clustering methods (e.g., [21]). For example, first, we arrange D˜i in the ascending order. Second, a fixed value δi is selected to be a quantile of diffD˜i. Third, for each s in D˜i, we counter the number of elements of D˜i in the interval s−δi,s+δi. Thus, we obtain a positive integer sequence, where each member indicates the relative size of the density of the corresponding element of D˜i. Finally, we select out those s in D˜i such that the corresponding integer number is either a local maximum or beyond a threshold.In the following experiments, we will setδi as two times the p-th percentile of set diffDi for some p, where Di is the rearrangement of D˜i in the ascending order and diffDi is the sequence of the first-order difference on Di. Denote N as the cardinality of D˜i. If N is small, we usually choose a smaller p; for example, p=75. If N is very large, we choose a bigger p; e.g., p=99. Otherwise, we can choose a proper p between them; e.g., p=90.Now, we have obtainedd+1 sets S0,S1,…,Sd with each one containing all the higher-density points of the corresponding distance set. By the triangle inequality, we have the following property:If there is an
i
∈
0,1
,
…
,
d
such that
d
i
=
d
X
,
O
i
is not in
S
i
, then X cannot be a higher-density point of
D
˜.According to this property, we can select a subsetS of D˜ whose distances to the i-th observation point are in Si for all i∈0,1,…,d.For each pointX∈D˜, we have mapped it to a d+1-dimensional vector. By Remark 4 of Theorem 2, all the d+1 pairs of the corresponding components of two points that belong to the same cluster will have a little difference between them. But it is possible that there are some data points which have some components that have the little difference with that of one cluster center and have the other components that have the little difference with that of another cluster center. In such case, few outliers may be missed by the above selection criterion. To discard those few outliers in S and decrease the number of elements of S, we need to refine S. We have the following criterion according to Remark 4 of Theorem 2.Suppose
X
1
has been selected. Given a data point
X
2
, if, for every
i
∈
0,1
,
…
,
d
,
d
1
i
=
d
X
1
,
O
i
and
d
2
i
=
d
X
2
,
O
i
have a small difference, then we discard
X
2.We denoteS=Y1,Y2,…,Ym and set(21)Sc=∅,Sa=∅.We also need a counter to indicate the density of each data point ofSa. Firstly, we let Y1∈Sa and make the indicative number of Y1 to be 1. We then sequentially choose the data points in S and dynamically construct Sc and Sa according to the following process. Suppose we choose Yi from S, then we compute the distance between Yi and each data point in Sa. If there are some distances less than a threshold value δ, we add 1 to each of the counter of data point that corresponds to these distances and then discard Yi. Meanwhile, if the counter number of a data point in Sa is bigger than another threshold value n, then we remove this data point from Sa and add it into Sc. But if each one of the point in Sa has distance to Yi bigger than δ, we continue to check whether there is a point in Sc that has distance to Yi less than δ; we will discard Yi if there is any and we will add Yi to Sc if not. Finally, we obtain a set Sc that closely represents the original data set and the size of it is smaller than that of the original data set. Furthermore, all outliers of original data set are not included in this selected subset.


### 3.1.3. ClusteringSc and D˜ in the Second Stage

Since the selected set has discarded all outliers and has a smaller size than the original data set, the running time will decrease significantly when performing thek-means algorithm on the selected subset. Furthermore, because the subset closely represents the original data set, the cluster centers will also be suitable to be chosen as the cluster centers of the original data set. When we have identified the cluster centers, it is easy to cluster the whole data set. The pseudo-code of our proposed algorithm is presented in Algorithm 1.Algorithm 1: Two-stage k-means clustering algorithm.
Input:
The number of clusters:k;
The number of percentile:p;
The original data set:D=Xi=xi1,xi2,…,xid:i=1,2,…,N;
Method:
NormalizeD and generate D˜=X˜i:i=1,2,…,N;
for t=0 to ddo
SetD˜t=dit:i=1,2,…,N where dit=dX˜i,Ot;
GenerateDt by rearranging D˜t in ascending order and set δt be p-th percentile of diffDt;
Setmt=minDt and Mt=maxDt/δt+1∗δt;
Equally divide the intervalmt,Mt into intervals with length δt;
LetSt be the union of those intervals that contain local maximum number of elements in Dt;
end for
Select out a subsetS:Xi∈S if and only if dit∈St for all t=0,1,…,d;
RefineS and obtain the subsets Sc:
Perform thek-means algorithm on Sc;
Assign eachX˜ii=1,2,…,N to the nearest center of the obtained cluster of Sc;
Output:
The result of clustering.


## 3.1.1. Generatingd+1 Distance Sets in the First Stage

First of all, we conduct the normalization operation on the original data setD. Set D=X1,X2,…,XN, where Xj=xj1,xj2,…,xjd for j=1,2,…,N. Suppose(19)M=maxj=1,2,…,N;i=1,2,…,dxjiandm=minj=1,2,…,N;i=1,2,…,dxji.Then, we transformD into D˜ with Xj=xj1,xj2,…,xjd corresponding to X˜j=xj1−m/M−m,xj2−m/M−m,…,xjd−m/M−m. Obviously, the transformation on D is a composition of a translation transformation and a dilation transformation. The dilation factor 1/M−m is the same for each dimension; hence, the dilation transformation does not change the cluster structure. Because the translation transformation also does not change the cluster structure of a dataset, the cluster structure of D is totally the same as that of D˜. We also note that the value of every component of X˜j is in the interval 0,1.Let(20)O0=0,0,…,0,Oi=0,…,0i−1,1i,0i+1,…,0,i=1,2,…,d,be the set of observation points. Denote D˜j the generated distance set of D˜ with respect to the observation point Ojj=0,1,…,d, and we get d+1 sets D˜0,D˜1,…,D˜d. For each data point X˜, we actually have mapped it to a d+1-dimensional vector dX˜,O0,dX˜,O1,…,dX˜,Od. Theorem 1 shows that we can identify X˜ by the d+1-dimensional vector, and hence, it is reasonable to expect that the clustering structure about D˜ can be deduced by those d+1 distance sets.


## 3.1.2. Selecting a Representative Subset ofD˜ in the First Stage

For eachD˜i, we can get a set Si consisting of all candidate higher density points of D˜i by using the grid-based clustering methods (e.g., [21]). For example, first, we arrange D˜i in the ascending order. Second, a fixed value δi is selected to be a quantile of diffD˜i. Third, for each s in D˜i, we counter the number of elements of D˜i in the interval s−δi,s+δi. Thus, we obtain a positive integer sequence, where each member indicates the relative size of the density of the corresponding element of D˜i. Finally, we select out those s in D˜i such that the corresponding integer number is either a local maximum or beyond a threshold.In the following experiments, we will setδi as two times the p-th percentile of set diffDi for some p, where Di is the rearrangement of D˜i in the ascending order and diffDi is the sequence of the first-order difference on Di. Denote N as the cardinality of D˜i. If N is small, we usually choose a smaller p; for example, p=75. If N is very large, we choose a bigger p; e.g., p=99. Otherwise, we can choose a proper p between them; e.g., p=90.Now, we have obtainedd+1 sets S0,S1,…,Sd with each one containing all the higher-density points of the corresponding distance set. By the triangle inequality, we have the following property:If there is an
i
∈
0,1
,
…
,
d
such that
d
i
=
d
X
,
O
i
is not in
S
i
, then X cannot be a higher-density point of
D
˜.According to this property, we can select a subsetS of D˜ whose distances to the i-th observation point are in Si for all i∈0,1,…,d.For each pointX∈D˜, we have mapped it to a d+1-dimensional vector. By Remark 4 of Theorem 2, all the d+1 pairs of the corresponding components of two points that belong to the same cluster will have a little difference between them. But it is possible that there are some data points which have some components that have the little difference with that of one cluster center and have the other components that have the little difference with that of another cluster center. In such case, few outliers may be missed by the above selection criterion. To discard those few outliers in S and decrease the number of elements of S, we need to refine S. We have the following criterion according to Remark 4 of Theorem 2.Suppose
X
1
has been selected. Given a data point
X
2
, if, for every
i
∈
0,1
,
…
,
d
,
d
1
i
=
d
X
1
,
O
i
and
d
2
i
=
d
X
2
,
O
i
have a small difference, then we discard
X
2.We denoteS=Y1,Y2,…,Ym and set(21)Sc=∅,Sa=∅.We also need a counter to indicate the density of each data point ofSa. Firstly, we let Y1∈Sa and make the indicative number of Y1 to be 1. We then sequentially choose the data points in S and dynamically construct Sc and Sa according to the following process. Suppose we choose Yi from S, then we compute the distance between Yi and each data point in Sa. If there are some distances less than a threshold value δ, we add 1 to each of the counter of data point that corresponds to these distances and then discard Yi. Meanwhile, if the counter number of a data point in Sa is bigger than another threshold value n, then we remove this data point from Sa and add it into Sc. But if each one of the point in Sa has distance to Yi bigger than δ, we continue to check whether there is a point in Sc that has distance to Yi less than δ; we will discard Yi if there is any and we will add Yi to Sc if not. Finally, we obtain a set Sc that closely represents the original data set and the size of it is smaller than that of the original data set. Furthermore, all outliers of original data set are not included in this selected subset.


## 3.1.3. ClusteringSc and D˜ in the Second Stage

Since the selected set has discarded all outliers and has a smaller size than the original data set, the running time will decrease significantly when performing thek-means algorithm on the selected subset. Furthermore, because the subset closely represents the original data set, the cluster centers will also be suitable to be chosen as the cluster centers of the original data set. When we have identified the cluster centers, it is easy to cluster the whole data set. The pseudo-code of our proposed algorithm is presented in Algorithm 1.Algorithm 1: Two-stage k-means clustering algorithm.
Input:
The number of clusters:k;
The number of percentile:p;
The original data set:D=Xi=xi1,xi2,…,xid:i=1,2,…,N;
Method:
NormalizeD and generate D˜=X˜i:i=1,2,…,N;
for t=0 to ddo
SetD˜t=dit:i=1,2,…,N where dit=dX˜i,Ot;
GenerateDt by rearranging D˜t in ascending order and set δt be p-th percentile of diffDt;
Setmt=minDt and Mt=maxDt/δt+1∗δt;
Equally divide the intervalmt,Mt into intervals with length δt;
LetSt be the union of those intervals that contain local maximum number of elements in Dt;
end for
Select out a subsetS:Xi∈S if and only if dit∈St for all t=0,1,…,d;
RefineS and obtain the subsets Sc:
Perform thek-means algorithm on Sc;
Assign eachX˜ii=1,2,…,N to the nearest center of the obtained cluster of Sc;
Output:
The result of clustering.


## 3.2. Analysis of Computational Complexity

In this section, we analyze the computational complexity of the proposed algorithm. When running the classicalk-means algorithm, each iteration needs to compute the distances between each data point in the whole data and those new modified cluster centers, which has a time complexity of ONkd. In our algorithm, the time cost in the first stage mainly consists of four parts. The first part is to generate d+1 one-dimensional data sets, which has a time complexity of ONd. The second one is to find those intervals which contain the local maximum of distances, which has a time complexity of ON. The third part is to select S, which has the time complexity ONd. In the fourth part, we refine S and obtain the subset Sc, which has the time complexity ON˜2, where N˜ denotes the cardinality of S. Thus, the time complexity of the first stage is ONd+N˜2. At the second stage, we will perform the k-means algorithm on Sc and each iteration will have a time complexity less than N˜1kd, where N˜1 denotes the cardinality of Sc. Because N˜1<N˜, the total time complexity of the new algorithm is ONd+N˜kd+N˜2.We note that the time complexity of the fourth part in the first stage is usually much less thanON˜2. Since many data points have been discarded when constructing the sets Sc and Sa, we do not have to compute the distances with all N˜ data points in S.


## 4. Experimental Results and Analysis

We have conducted a series of experiments on 6 synthetic data sets and 3 benchmark data sets (UCI [22] and KEEL [23]) to validate the effectiveness of the proposed two-stage k-means clustering algorithm in this section. The synthetic data sets can be downloaded from BaiduPan (https://pan.baidu.com/s/1MfS8JfQdJLHYSlpZdndLUQ) with the extraction code “p3mc.” We first present the clustering results of our proposed algorithm and the k-means algorithm on two synthetic data sets, i.e., the data set #1 and data set #2. The experimental results are shown in Figures 1 and 2. For simplicity, we only use the experimental results on data set #1 to explain the advantage of our proposed algorithm. There are two clusters in data set #1, where each cluster includes 41 data points. The data points obey the 2-dimensional normal distributions with mean vectors (3, 11) and (12, 5) and covariance matrices 3002 and 2003, respectively. There are also two outliers in data set #1.Figure 1
Experimental results on synthetic data set #1 (the red ‘+’s are the cluster centers). (a) Normalized data. (b) Selected subset. (c) Clustering centers ofk-means algorithm. (d) Clustering centers of our proposed algorithm.
(a)
(b)
(c)
(d)Figure 2
Experimental results on synthetic data set #2 (the red ‘+’s are the cluster centers). (a) Normalized data. (b) Selected subset. (c) Clustering centers ofk-means algorithm. (d) Clustering centers of our proposed algorithm.
(a)
(b)
(c)
(d)Figure1(b) gives the selected data points of normalized data points corresponding to the data set #1 as shown in Figure 1(a). In Figure 1(b), we can find that outliers have been removed in the first stage of our proposed method. Figure 1(c) shows the clustering result of the k-means algorithm. We can see that outliers seriously impact the clustering result of the k-means algorithm, although there are only two outlier data points in the data set #1. The clustering result of our proposed method is presented in Figure 1(d), where the cluster center can be found correctly without the disturbance of outliers. The similar results can be found in Figure 2 for the data set #2 which includes 7 clusters and 10 outliers. The experimental results reflect that our proposed two-stage k-means clustering algorithm is not sensitive to outliers and can obtain the better clustering results than that of k-means clustering algorithm.Furthermore, we choose another four synthetic data sets as shown in Figure3 (only 2-dimensional illustration) and three real-world data sets to compare the clustering performances of our proposed algorithm with the k-means algorithm. The details of these data sets and experimental results are summarized in Table 1, where N is the number of the elements of the data set, t is the proportion of the outlier in the data set, k is the number of clusters, d is the dimension of data point, p is the percentile number, nc is the cardinality of selected subset, ARIkmeans and Timekmeans are the adjusted Rand index (ARI) and time consumption of k-means algorithm, and ARIour and Timeour are ARI and time consumption of our proposed algorithm. In Table 1, we can see that our proposed algorithm obtains the larger ARIs with the lower time consumption in comparison with k-means clustering algorithm on these synthetic data sets. For the real data sets without outliers, our algorithm can obtain the ARIs comparable to that of k-means algorithm. Nevertheless, the ARIs of k-means algorithm are severely degraded when the outliers are deliberately arranged in the real data sets, while the experimental results in Table 1 demonstrate that our proposed clustering algorithm is robust to the outliers. Table 2 shows the details of comparison on four large-scale synthetic data sets. The variables in Table 2 have the same meaning as that in Table 1. The comparison of time complexity between our proposed algorithm and k-means algorithm in Table 2 reflects that our algorithm has less time consumption than k-means algorithm. Especially, we can find that the superiority of our proposed method on time consumption is more obvious for data set with the larger size and dimension. Furthermore, the most time-consuming procedure in our algorithm, i.e., the selection of high-density distances for each generated distance set can be ran in the parallel way, which make our algorithm to be easily extended to perform the clustering task for large-scale data set.Figure 3
Two-dimensional illustrations corresponding to four synthetic data sets (the red data points are outliers.). (a) Data set 3 (five clusters). (b) Data set 4 (six clusters). (c) Data set 5 (six clusters). (d) Data set 6 (seven clusters).
(a)
(b)
(c)
(d)Table 1
Comparison betweenk-means clustering algorithm and our proposed clustering algorithm.
Data sets
N
t (%)
k
d
p
n
c
ARI
k
means
ARI
our
Time
k
means
Time
our
Synthetic #1
84
2.3
2
2
75
27
0.084
0.954
0.020
0.000
Synthetic #2
236
4.2
7
2
80
107
0.791
0.860
0.063
0.016
Synthetic #3
2557
2.2
5
3
92
761
0.774
0.972
0.047
0.031
Synthetic #4
3670
1.9
6
4
96
656
0.815
0.977
0.188
0.063
Synthetic #5
3655
1.5
6
5
96
573
0.816
0.982
0.313
0.078
Synthetic #6
2830
1.1
7
6
95
296
0.848
0.988
0.250
0.063
Iris
150
0
3
4
80
27
0.730
0.730
0.031
0.016
Iris∗
152
1.3
3
4
82
26
0.531
0.743
0.031
0.016
Seeds
210
0
3
7
80
37
0.717
0.728
0.047
0.016
Seeds∗
212
0.9
3
7
80
37
0.462
0.694
0.047
0.016
Wine
178
0
3
13
81
6
0.870
0.850
0.031
0.016
Wine∗
180
1.1
3
13
81
10
0.365
0.882
0.031
0.016
∗The real data set which includes two synthetic outliers as shown in our BaiduPan.Table 2
Comparison betweenk-means clustering algorithm and our proposed clustering algorithm on large-scale data sets.
Data sets
N
k
d
p
n
c
ARI
k
means
ARI
our
Time
k
means
Time
our
Synthetic #7
22501
5
4
99
1413
1
1
0.858
0.718
Synthetic #8
50001
5
4
99.5
1148
1
1
1.732
1.248
Synthetic #9
60001
6
4
99.5
1257
1
1
2.262
1.466
Synthetic #10
60001
6
5
99.5
1038
1
1
2.387
1.810In addition, we provide a real application, i.e., the tyre inclusion identification, to validate the clustering performance of our proposed clustering algorithm. Figure4 shows two tyres with different kinds of inclusions, where each picture includes 1027 × 768 pixels. Figures 5 and 6 present the clustering results of our proposed algorithm and k-means clustering algorithm on Tyre #1 and Tyre #2, respectively. In these figures, we can see that our proposed method can accurately identify the cluster centers without the disturbance of outliers. The inclusions can be clearly recognized by our proposed algorithm in the tyres, while the k-means clustering algorithm does not find the inclusions distinctly, e.g., Figures 6(b), 6(d), and 6(f) include not only the inclusions but also the tyre traces. Above all, the experimental results demonstrate the better clustering performance in comparison with the classical k-means clustering algorithm when handling the clustering tasks with the disturbance of outliers.Figure 4
Application of tyre inclusion identification. (a) Tyre #1. (b) Tyre #2.
(a)
(b)Figure 5
Clustering results on Tyre #1 (318×80 pixels, p=99, N=25440, and nc=7696). (a) Our algorithm with k = 2. (b) k-means with k = 2. (c) Our algorithm with k = 3. (d) k-means with k = 3. (e) Our algorithm with k = 4 (f) k-means with k = 4.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Clustering results on Tyre #2 (345× 152 pixels, p=99, N=52440, and nc=16852). (a) Our algorithm with k = 2. (b) k-means with k = 2. (c) Our algorithm with k = 3. (d) k-means with k = 3. (e) Our algorithm with k = 4. (f) k-means with k = 4.
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Conclusions and Future Work

In this paper, we proposed a robust two-stagek-means clustering algorithm which can accurately identify the cluster centers without the disturbance of outliers. As the direct application of the observation point mechanism of I-nice [18], we select a small subset from the original data set based on a set of nondegenerate observation points in the first stage. In the second stage, we use the k-means clustering algorithm to cluster the selected subset and make these cluster centers as the true cluster centers of the original data set. The theoretical analysis and experimental verification demonstrate the feasibility and effectiveness of proposed clustering algorithm. The future studies will be focused on three directions. First, we will try to use the k-nearest neighbors (kNN) method to improve the selection of observation points. Second, we will seek the real applications for the two-stage k-means clustering algorithm. Third, we will extend our proposed algorithm to cluster big data based on the random sample partition model [24].



---

*Source: 3650926-2020-03-30.xml*