# Asymptotics for The $k$-means

Tonglin Zhang 111Department of Statistics, Purdue University, 250 North
University Street,West Lafayette, IN 47907-2066, Email<EMAIL_ADDRESS>
###### Abstract

The $k$-means is one of the most important unsupervised learning techniques in
statistics and computer science. The goal is to partition a data set into many
clusters, such that observations within clusters are the most homogeneous and
observations between clusters are the most heterogeneous. Although it is well
known, the investigation of the asymptotic properties is far behind, leading
to difficulties in developing more precise $k$-means methods in practice. To
address this issue, a new concept called clustering consistency is proposed.
Fundamentally, the proposed clustering consistency is more appropriate than
the previous criterion consistency for the clustering methods. Using this
concept, a new $k$-means method is proposed. It is found that the proposed
$k$-means method has lower clustering error rates and is more robust to small
clusters and outliers than existing $k$-means methods. When $k$ is unknown,
using the Gap statistics, the proposed method can also identify the number of
clusters. This is rarely achieved by existing $k$-means methods adopted by
many software packages.

Key Words: Clustering Consistency; Criterion Consistency; Gaussian Mixture
Models; $k$-means; Sum-of-Squares Criterion; Unsupervised Learning.

AMS Classification: 62H30; 62J12;

##  1 Introduction

Clustering is one of the most important unsupervised learning techniques for
understanding the underlying data structures. The goal is to partition a data
set into many subsets, called clusters, such that the observations within the
subsets are the most homogeneous and the observations between the subsets are
the most heterogeneous. Clustering is usually carried out by specifying a
similarity or dissimilarity measure between observations. Examples include the
$k$-means [17, 19, 29, 37], the $k$-medians [3], the $k$-modes [5], and the
generalized $k$-means [2, 31, 45], as well as many of their modifications [21,
24, 42]. Among those, the $k$-means has been considered as one of the most
straightforward and popular methods since it was proposed sixty years ago [23,
36]. Although it is well known, the investigation of the theoretical
properties is still far behind, leading to difficulties in developing more
precise $k$-means methods in practice. The goal of the present research is to
propose a new concept called clustering consistency for the asymptotics of the
$k$-means with a resulting clustering method better than the existing
$k$-means methods adopted by many software packages, including those adopted
by R and Python.

The research has three main contributions. The first is the development of
clustering consistency, which is more appropriate than the previous criterion
consistency used for decades for asymptotics of clustering methods. In the
literature, theoretical properties of the $k$-means are investigated by the
performance of the sum-of-squares (SSQ) criterion function, which is
classified as criterion consistency by the article. The difference is that
clustering consistency focuses on the asymptotic properties of the partition
resulted by the $k$-means but not criterion consistency focuses on the
asymptotic properties of the SSQ criterion used by the $k$-means. The second
is the derivation of a better $k$-means method using clustering consistency.
It is better because any local solution of the proposed $k$-means method is
clustering consistent, implying that we successfully avoid the NP-hardness
issue [30] for the global solution of the $k$-means problem. This property is
violated if other $k$-means methods are used. The third is estimation of the
number of clusters when $k$ is unknown. This is a major concern in practice.
Our method can be used to estimate the number of clusters if it is combined
with the well known Gap statistics [39].

Among the large amount of the literature on the categories of unsupervised
clustering learning, the $k$-means is considered as the simplest and also
stands out to be the most popular one. It belongs to the category of
partitioning clustering, which can be interpreted by the centroidal Voronoi
tessellation method in mathematics [8]. Besides the $k$-means, partitioning
clustering also includes the $k$-medians, the $k$-modes, and the generalized
$k$-means. Other clustering categories proposed in the literature include
hierarchical clustering [46], fuzzy clustering [40], density-based clustering
[9, 25], and model-based clustering [15, 26].

The $k$-means is an iterative method that tries to partition a data set into
$k$ distinct non-overlapping and non-empty clusters. In particular, let ${\cal
D}=\\{{\bm{x}}_{1},\dots,{\bm{x}}_{n}\in\mathbb{R}^{p}\\}$ be the data set
composed by $n$ observations and
$\mathscr{C}_{k}=\\{C_{1},\dots,C_{k}:\bigcup_{r=1}^{k}C_{r}={\cal
D},C_{r}\not=\emptyset,C_{r}\cap C_{r^{\prime}}=\emptyset,r\not=r^{\prime}\\}$
be a size-$k$ partition of ${\cal D}$. The $k$-means estimates the partition
by

$\widehat{\mathscr{C}_{k}}=\\{\hat{C}_{1},\dots,\hat{C}_{k}\\}=\mathop{\arg\\!\min}_{\mathscr{C}_{k}}SSQ(\mathscr{C}_{k}),$
(1)

where

$SSQ(\mathscr{C}_{k})=\sum_{r=1}^{k}\sum_{i\in{C_{r}}}\|{\bm{x}}_{i}-{\bm{c}}_{r}\|^{2},$
(2)

is the SSQ criterion and ${\bm{c}}_{r}\in\mathbb{R}^{p}$ is the centroid of
$C_{r}$. The $k$ centroids ${\bm{c}}_{1},\dots,{\bm{c}}_{k}$ are treated as
unknown parameters. In practice, $\widehat{\mathscr{C}_{k}}$ and
$\hat{\bm{c}}_{1},\dots,\hat{\bm{c}}_{k}$ (i.e., the corresponding estimates
of the $k$ centroids) are computed iteratively. A $k$-means method starts with
initial guesses of $\hat{\bm{c}}_{1},\dots,\hat{\bm{c}}_{k}$. They are used to
compute the optimal partition of the data set. The centroids are updated after
that. By iterations, the final estimates of the partition and the centroids
are derived. It has been pointed out by many researchers that the final answer
may not be the global solution because the output depends on the initial
guesses of the centroids. To avoid a local optimizer, many researchers
recommend performing the algorithm multiple times with different initial
guesses of the centroids [10, 17]. It has been pointed out that this may
induce thousands of replications in searching the global optimizer even when
$n$ is a few thousand [35].

The theoretical properties of $k$-means algorithms highly depend on the
initialization of the centroids. The global solution of (1) is NP-hard [30].
To avoid this issue, many researchers investigate local solutions. The
performance of the basic (i.e., the earliest) $k$-means algorithm is bad
because it uses randomly generated initial centroids. To overcome the
difficulty, methods using initial centroids far away from each other have been
proposed. Among those, the $k$-means++ [1], which is the default of python, is
considered as one of the best methods. Although it is not the default of R,
the $k$-means++ has been incorporated by library LICORS. Although the
performance of the $k$-means++ is better than that of the basic $k$-means, we
identify a better $k$-means method using clustering consistency.

The article is organized as follows. In Section 2, we review previous
criterion consistency for the $k$-means. In Section 3, we propose clustering
consistency for the $k$-means and use it to develop a new $k$-means method
under the $k$-means model. In Section 4, we extend our findings beyond the
$k$-means model. In Section 5, we evaluate the performance of our method with
the comparison to the previous $k$-means methods by Monte Carlo simulations.
In Section 6, we apply our method to a real world data set. In Section 7, we
provide a discussion. We put all of the proofs in the Appendix.

##  2 Related Work

In the literature, asymptotics for the $k$-means is evaluated by the
properties of the SSQ-criterion. It is classified as criterion consistency in
the article. Rather than the proposed clustering consistency which focuses on
the local solution of the $k$-means, the previous criterion consistency
focuses on the the global solutions only. The interest is to show the final
$SSQ(\widehat{\mathscr{C}_{k}})$ (assume that it is the global solution)
satisfies

${1\over
n}[SSQ(\widehat{\mathscr{C}_{k}})-SSQ(\mathscr{C}_{k,T})]\stackrel{{\scriptstyle
P}}{{\rightarrow}}0,{\rm as}\ n\rightarrow\infty,$ (3)

where
$SSQ(\mathscr{C}_{k,T})=\sum_{r=1}^{k}\sum_{i\in{C_{r,T}}}\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|^{2}$,
${\bm{\mu}}_{r}$ is the true center of the $r$th cluster, and
$\mathscr{C}_{k.T}$ is the true partition.

Based on criterion consistency, as early as in 1967, MacQueen [29] obtained
weak consistency for the $k$-means. He showed that the corresponding SSQ
criterion converges almost surely, but he did not prove consistency of the
cluster centroids. In 1978, Hartigan [18] obtained consistency and asymptotic
normality of the SSQ criterion in the one-dimensional case when the
observations are splitted into two clusters. In 1981, Pollard [33] provided
conditions for almost sure convergence of the cluster centers as the sample
size increases to infinity in the high-dimensional case. In 1984, Wong [43]
obtained consistency of the univariate $k$-means when $k$ goes to infinity
with $n$. In 2015, Terada [38] provided strong consistency of the SSQ
criterion when the principal component analysis (PCA) approach is applied to
the $k$-means for dimension reduction. In 2022, Chakrabary and Das [4]
established strong consistency of estimators of cluster centroids. In all of
these, the interest is consistency of the SSQ-criterion and the estimators of
the centroids but not the estimator of the partition itself.

To implement (3) for criterion consistency, it is necessary to provide a nice
definition for true clusters and true clustering. It is well known that there
is no agreed definition of what the true clustering is if only the observed
data set ${\cal D}$ is used. To interpret this concept, a previous method
treats the global solution of (1) as the true clustering, but this is
considered as an inappropriate approach [20]. In the situation of clustering
analysis, it is assumed that the true cluster labels are unknown. The goal of
a clustering method is to assign labels to all ${\bm{x}}_{i}$. Without using
the true labels, it is hard to define true partition and true clustering.
Therefore, we assume that the true cluster labels are disclosed in the
definition of $\mathscr{C}_{k,T}$ although they cannot be used to formulate a
usable clustering method, meaning that we treat true clustering as an unusable
clustering method. Using this idea, we obtain the formal definition of the
true clustering and true partition. If this is adopted, then the criterion
consistency given by (3) is well defined. We provide interpretations of the
previous criterion consistency and the proposed clustering consistency.

##  3 Main Result

We use a complete data model and an incomplete data model for the development
of our method. In the complete data model, we assume that the true cluster
labels are disclosed. Each ${\bm{x}}_{i}$ is associated with an indicator
vector ${\bm{z}}_{i}=(Z_{i1},\dots,Z_{ik})^{\top}$ with $Z_{ir}\in\\{1,0\\}$
and $\sum_{r=1}^{k}Z_{ir}=1$ for its cluster assignment. The complete data
model assums that both ${\bm{x}}_{i}$ and ${\bm{z}}_{i}$ are available,
leading to the complete date set as ${\cal
D}_{c}=\\{({\bm{x}}_{i}:{\bm{z}}_{i}):i=1,\cdots,n\\}$. The incomplete data
model assumes that the true cluster labels are not observed. It treats
${\bm{z}}_{i}$ as unobserved latent variables, leading to the incomplete data
set as ${\cal D}=\\{{\bm{x}}_{i}:i=1,\dots,n\\}$. We use ${\cal D}_{c}$ to
define the true clustering. Because a usable clustering method can only use
${\cal D}$, the true clustering is treated as an unusable method. The goal of
clustering consistency is to show that the estimate of the partition resulted
by a usable clustering method is asymptotically identical to the partition
resulted by the true clustering method, indicates that it is asymptotically
equivalent to the case when the cluster assignments are disclosed. Due to the
NP-hardness issue of the global solution of the $k$-means problem, we provide
a new $k$-means method with clustering consistency satisfied by its local
solutions, implying that the NP-hardness issue is not a concern in our method.

###  3.1 Gaussian Mixture Models

We evaluate the asymptotic properties of the general $k$-means method and
develop our new $k$-means method under the framework of the Gaussian mixture
models (GMMs). Theoretically, a GMM is a special case of mixture models. A
mixture model is defined by a hierarchical modeling approach. The first level
specifies the distributions of ${\bm{z}}_{i}$ for the ground truth. The second
level specifies the distribution of the observations given ${\bm{z}}_{i}$. The
mixture model can be interpreted by either a complete data model or an
incomplete data model. The complete data model assumes that the complete data
set ${\cal D}_{c}$ is available, leading to the underlying distribution as

${\bm{x}}_{i}|{\bm{z}}_{i}\sim^{iid}\sum_{r=1}^{k}Z_{ir}f_{r}({\bm{x}}_{i})=\prod_{r=1}^{k}[f_{r}({\bm{x}}_{i})]^{Z_{ir}}.$
(4)

and the loglikelihood function as

${\bm{\ell}}_{c}({\bm{\theta}})=\sum_{i=1}^{n}\sum_{r=1}^{k}Z_{ir}\log[f_{r}({\bm{x}}_{i})],$
(5)

where ${\bm{\theta}}$ the vector composed by all parameters involved in the
model.

The incomplete data model treats ${\bm{z}}_{i}$ as unobserved latent
variables. It is assumed that only ${\cal D}=\\{{\bm{x}}_{i}:i=1,\dots,n\\}$
is available. Suppose that ${\bm{z}}_{i}$ are iid from a Dirichlet
distribution with probability vector
${\bm{\pi}}=(\pi_{1},\dots,\pi_{k})^{\top}$. By integrating out the latent
variables ${\bm{z}}_{i}$ from (4), the underlying distribution of the
incomplete data model becomes
${\bm{x}}_{i}\sim^{iid}\sum_{r=1}^{k}\pi_{r}f_{r}({\bm{x}}_{i})$. A GMM is
derived if $f_{r}({\bm{x}}_{i})$ is specified the PDF of ${\cal
N}({\bm{\mu}}_{r},{\bm{\Sigma}}_{r})$ as

$f_{r}({\bm{x}}_{i})=\varphi({\bm{x}}_{i};{\bm{\mu}}_{r},{\bm{\Sigma}}_{r})=-{1\over(2\pi)^{p\over
2}|{\det({\bm{\Sigma}}_{r})|^{1\over 2}}}e^{-{1\over
2}({\bm{x}}_{i}-{\bm{\mu}}_{r})^{\top}{\bm{\Sigma}}_{r}^{-1}({\bm{x}}_{i}-{\bm{\mu}}_{r})}$
(6)

leading to underlying model in the incomplete data model as

${\bm{x}}_{i}\sim^{iid}\sum_{r=1}^{k}\pi_{r}{\cal
N}({\bm{\mu}}_{r},{\bm{\Sigma}}_{r}),$ (7)

where ${\bm{\mu}}_{r}$ is the mean vector (i.e., the true center) and
${\bm{\Sigma}}_{r}$ is the variance-covariance matrix of the $r$th cluster.

The mixture probability (i.e., weight) $\pi_{r}$ is positive and satisfies
$\sum_{r=1}^{k}\pi_{r}=1$. The GMM treats $\pi_{r}$, ${\bm{\mu}}_{r}$, and
${\bm{\Sigma}}_{r}$ as unknown parameters. By treating the latent variables as
missing, the MLE of the parameters can be computed by the EM-algorithm [7]. If
${\bm{\Sigma}_{r}}=\sigma^{2}{\bf I}$ for all $r=1,\dots,k$, then the GMM
becomes the $k$-means model. Because the EM algorithm is a method for the MLE
when missing values are present [44], it can also be used to study the
$k$-means problem. In this case, we need to assume that
${\bm{\mu}}_{1},\dots,{\bm{\mu}}_{k}$ are all different in (7) to be
consistent with the $k$-means problem.

We briefly review the EM algorithm for the GMM given by (4) with
$f_{r}({\bm{x}}_{i})$ specified by (6), which has bee formulated previously by
many authors [32, e.g.]. Let $\pi_{r}^{(t)}$, ${\bm{\mu}}_{r}^{(t)}$, and
${\bm{\Sigma}}_{r}^{(t)}$ be the $t$th iterative values of the MLEs of
$\pi_{r}$, ${\bm{\mu}}_{r}$, and ${\bm{\Sigma}}_{r}$ (denoted by
$\hat{\pi}_{r}$, $\hat{\bm{\mu}}_{r}$, and $\hat{\bm{\Sigma}}_{r}$,
respectively) in the EM for the GMM defined by (7). The E-step of the EM-
algorithm predicts ${\bm{z}}_{i}$ as

$Z_{ir}^{(t+1)}={\pi_{r}^{(t)}\varphi({\bm{x}}_{i}|{\bm{\mu}}_{r}^{(t)},{\bm{\Sigma}}_{r}^{(t)})\over\sum_{j=1}^{k}\pi_{j}^{(t)}\varphi({\bm{x}}_{i}|{\bm{\mu}}_{j}^{(t)},{\bm{\Sigma}}_{j}^{(t)})}.$
(8)

The M-step of the EM-algorithm estimates the parameters as

$\displaystyle{\pi_{r}^{(t+1)}=}$ | $\displaystyle{{}{1\over n}\sum_{i=1}^{n}Z_{ir}^{(t+1)},}$
---|---
$\displaystyle{{\bm{\mu}}_{r}^{(t+1)}=}$ | $\displaystyle{{}{\sum_{i=1}^{n}Z_{it}^{(t+1)}{\bm{x}}_{i}\over\sum_{i=1}^{n}Z_{ir}^{(t+1)}},}$
$\displaystyle{{\bm{\Sigma}}_{r}^{(t+1)}=}$ | $\displaystyle{{}{\sum_{i=1}^{n}Z_{it}^{(t+1)}({\bm{x}}_{i}-{\bm{\mu}}_{r}^{(t+1)})^{\top}({\bm{x}}_{i}-{\bm{\mu}}_{r}^{(t+1)})\over\sum_{i=1}^{n}Z_{ir}^{(t+1)}}.}$
(9)

In the end, the EM-algorithm estimates the partition by

$\hat{C}_{r}=\\{i:\hat{Z}_{ir}=\mathop{\arg\\!\max}_{j\in\\{1,\dots,k\\}}\hat{Z}_{ij}\\},\
r=1,\dots,k,$ (10)

where $\hat{\bm{z}}_{i}$ is the $i$th final imputed vector of the ground truth
given by iterating (8) and (9). By taking ${\bm{\Sigma}}_{r}=\sigma^{2}{\bf
I}$ in the GMM, we obtain the $k$-means model with

$f_{r}({\bm{x}}_{i})=\varphi({\bm{x}}_{i};{\bm{\mu}}_{r},\sigma^{2}{\bf
I})=-{1\over(2\pi)^{p\over
2}\sigma^{p}}e^{-{\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|^{2}\over 2\sigma^{2}}}$ (11)

in (4), leading to the EM-algorithm for the $k$-means model.

The $k$-means method does not use the EM-algorithm. Instead, it directly
computes the next iterative value of ${\bm{z}}_{i}$ given the previous
${\bm{\mu}}_{1}^{(t)},\dots,{\bm{\mu}}_{k}^{(t)}$ as

$Z_{ir}^{(t+1)}=\left\\{\begin{array}[]{cc}1,&r=\mathop{\arg\\!\min}_{j\in\\{1,\dots,k\\}}\|{\bm{x}}_{i}-{\bm{\mu}}_{j}^{(t)}\|,\cr
0,&{\rm otherwise},\end{array}\right.$ (12)

and the next iterative value of the MLE of ${\bm{\mu}}_{r}$ as

${\bm{\mu}}_{r}^{(t+1)}={\sum_{i=1}^{n}Z_{ir}^{(t+1)}{\bm{x}}_{i}\over\sum_{i=1}^{n}Z_{ir}^{(t+1)}}.$
(13)

In the end, it estimates the partition by

$\hat{C}_{r}=\\{i:\hat{Z}_{ir}=1\\},\ r=1,\dots,k,$ (14)

where $\hat{\bm{z}}_{i}=(\hat{Z}_{i1},\dots,\hat{Z}_{ik})^{\top}$ is the final
imputed vector of the ground truth given by iterating (12) and (13). Neither
the EM-algorithm nor the $k$-means uses the ground truth ${\bm{z}}_{i}$ in
estimating the clusters. Instead, they use the imputed $\hat{\bm{z}}_{i}$.
Therefore, they are usable clustering methods.

###  3.2 Clustering Consistency

The usage of (10) by the EM-algorithm or (14) by a $k$-means method indicates
that we should use ${\bm{z}}_{i}$ to define the true partition and the true
clustering. In the literature, the definition of the true partition and true
clustering is considered as a concern in the study of theoretical properties
of clustering methods. Several approaches have been proposed previously [20],
but none of them use ${\bm{z}}_{i}$ to define the true partition or the true
clustering. To overcome the difficulty, we provide the following definition.

###### Definition 1

The true partition for a mixture model is
$\mathscr{C}_{k,T}=\\{C_{1,T},\dots,C_{k,T}\\}$ with
$C_{r,T}=\\{i:Z_{ir}=1\\}$ for $r=1,\dots,k$. A clustering method is said true
clustering if it uses ${\bm{z}}_{i}$ to assign the labels of clusters. A
clustering method is said usable if it does not use ${\bm{z}}_{i}$ to assign
the clusters.

The true clustering does not make any mistakes in assigning the cluster
labels. The true clustering given by Definition 1 is not usable because it
uses ${\bm{z}}_{i}$. It is generally impossible to formulate a clustering
method which can always result $\mathscr{C}_{k,T}$ without using
${\bm{z}}_{i}$ when $n$ is finite. Therefore, it is important to compare the
clusters assigned by a usable clustering method with those assigned by the
true clustering method as $n\rightarrow\infty$, leading to the clustering
consistency developed by the current research.

###### Definition 2

We say that a $k$-means method resulting $\widehat{\mathscr{C}_{k}}$ is
clustering consistent, denoted by
$\widehat{\mathscr{C}_{k}}\stackrel{{\scriptstyle
P}}{{\rightarrow}}\mathscr{C}_{k,T}$, if

$\lim_{n\rightarrow\infty}{\rm
Pr}(\widehat{\mathscr{C}_{k}}=\mathscr{C}_{k,T})=1.$ (15)

We say that the $k$-means method is clustering inconsistent if (15) is
violated.

Clustering consistency can be extended to arbitrary (usable) clustering
methods, because the only requirement is to ensure that (15) holds. Because
the labels of clusters are not specified in the data, clustering consistency
given by Definition 2 does not mean $\hat{C}_{r}\stackrel{{\scriptstyle
P}}{{\rightarrow}}C_{r,T}$ for every $r$. Instead, it only needs that there is
a permutation $\tau$ on $\\{1,\dots,k\\}$, such that
$\hat{C}_{\tau(r)}\stackrel{{\scriptstyle P}}{{\rightarrow}}C_{r,T}$ for every
$r$. Therefore, Definition 2 contains $k!$ labeling cases. In the case when
the ground truth variables are provided but not used in the $k$-means method,
to check whether clustering consistency holds, we need to find $\tau$. This
kind of problems can be solved by the well known Hungarin algorithm in
mathematics.

To obtain clustering consistency, we need to impose regularity conditions for
${\bm{\mu}}_{1},\dots,{\bm{\mu}}_{k}$ and $\sigma^{2}$. The following example
shows that the minimum distance between ${\bm{\mu}}_{1},\dots,{\bm{\mu}}_{k}$
should go to infinity as $n\rightarrow\infty$ if $\sigma^{2}$ is fixed.

Example 1: Suppose that $\sigma^{2}=1$, $k=2$, and $\pi_{1}=\pi_{2}=1/2$ in
(11). Then, the centroid of the first cluster is ${\bm{c}}_{1}={\rm
E}[{\bm{x}}|\|{\bm{x}}-{\bm{\mu}}_{1}\|\leq\|{\bm{x}}-{\bm{\mu}}_{2}\|]$ and
the center of the second cluster is ${\bm{c}}_{2}={\rm
E}[{\bm{x}}|\|{\bm{x}}-{\bm{\mu}}_{1}\|\geq\|{\bm{x}}-{\bm{\mu}}_{2}\|]$,
which are functions of ${\bm{\mu}}_{1}$, ${\bm{\mu}}_{2}$, and $\sigma^{2}$.
Both ${\bm{c}}_{1}$ and ${\bm{c}}_{2}$ are on the line segment between
${\bm{\mu}}_{1}$ and ${\bm{\mu}}_{2}$. Let $\hat{\bm{c}}_{1}$ and
$\hat{\bm{c}}_{2}$ be the estimates of ${\bm{c}}_{1}$ and ${\bm{c}}_{2}$
provided by the global solution. Then,
$(\hat{\bm{c}}_{1}+\hat{\bm{c}}_{2})/2\stackrel{{\scriptstyle
P}}{{\rightarrow}}({\bm{\mu}_{1}}+{\bm{\mu}}_{2})/2$ as $n\rightarrow\infty$
by the asymptotic theory of M-estimation. Asymptotically, the $k$-means
assigns ${\bm{x}}_{i}$ to the first cluster if
$\|{\bm{x}}_{i}-{\bm{\mu}}_{1}\|\leq\|{\bm{x}}_{i}-{\bm{\mu}}_{2}\|$ or the
second cluster otherwise. Because ${\rm
Pr}\\{\|{\bm{x}}-{\bm{\mu}}_{1}\|\leq\|{\bm{x}}-{\bm{\mu}}_{2}\||{\bm{x}}\in
C_{2,T}\\}={\rm
Pr}\\{\|{\bm{x}}-{\bm{\mu}}_{2}\|\leq\|{\bm{x}}-{\bm{\mu}}_{1}\||{\bm{x}}\in
C_{1,T}\\}$ is positively bounded from $0$ asymptotically, we do not have
$\hat{C}_{\tau(1)}\stackrel{{\scriptstyle P}}{{\rightarrow}}C_{1,T}$ and
$\hat{C}_{\tau(2)}\stackrel{{\scriptstyle P}}{{\rightarrow}}C_{2,T}$ as
$n\rightarrow\infty$, implying that any $k$-means methods cannot be clustering
consistent under this setting.

The main issue in Example 1 is that $\|{\bm{\mu}}_{1}-{\bm{\mu}}_{2}\|$ is
bounded when $\sigma^{2}$ is fixed as $n\rightarrow\infty$, leading to
clustering inconsistency of all possible $k$-means methods. When $n$ becomes
large, the size of the first cluster grows. Sooner or later, it will cover
${\bm{\mu}}_{2}$. This also happens in the second cluster. Therefore, to
obtain clustering consistency, we need to assume that

$d=\min\\{\|{\bm{\mu}}_{r}-{\bm{\mu}}_{r^{\prime}}\|:r\not=r^{\prime}\\}\rightarrow\infty$
(16)

as $n\rightarrow\infty$. This has also been pointed out in the case when the
$k$-means is applied to group functional data in two clusters (i.e., $k=2$) by
a concept called asymptotic perfectness [6]. The difference is that asymptotic
perfectness focuses on the global solution but clustering consistency also
investigates the local solutions. Using clustering consistency, we conclude if
the speed for $d$ to approach $\infty$ is faster than that of $\log^{1/2}{n}$
(when $\sigma^{2}$ is fixed), then any location solution provided by the
proposed $k$-means method (i.e., given by Algorithm 1) is clustering
consistent. This property may not be satisfied if other $k$-means methods are
used. Thus, we can use local solutions of our method. If other $k$-means
method are used, then we should use the global solution. Then, we have the
following theorems.

###### Theorem 1

If (a) ${\bm{x}}_{1},\dots,{\bm{x}}_{n}$ are independently collected from (7)
with ${\bm{\Sigma}}_{r}=\sigma^{2}{\bf I}$, (b) the initial centroids are
provided by data points in different clusters, (c)
$\mathop{\lim\\!\inf}_{n\rightarrow\infty}d/[4\sigma(2\log{n})^{1/2}]>1$, and
(d) $\lim_{n\rightarrow\infty}p\log[d^{2}/(p\sigma^{2})]/\log{n}=0$, then the
$k$-means method is clustering consistent.

Theorem 1 can be implemented even if $\sigma^{2}$ and $p$ vary with $n$.
Condition (d) is not needed if $p$ is a constant. Therefore, Condition (c) is
more important. It has two special cases. The first case assumes that
$\sigma^{2}$ is a constant. We need that $d$ grows to infinity with $n$. The
second assumes that $d$ is a constant. We need that $\sigma^{2}$ goes to $0$
as $n$ goes to infinity. We have the following Corollary.

###### Corollary 1

Assume that Conditions (a) and (b) of Theorem 1 are satisfied, and $p$ is a
constant. If (e) (i) $\sigma^{2}$ is a constant and
$\lim\\!\inf_{n\rightarrow\infty}d/(\log{n})^{1/2}>4\sqrt{2}\sigma$ or (ii)
$d$ is a constant and
$\lim\\!\sup_{n\rightarrow\infty}\sigma^{2}\log{n}<d^{2}/32$, then the
$k$-means method is clustering consistent.

When $p$ and $\sigma^{2}$ does not vary with $n$, we can obtain a $k$-mean
method with clustering consistency if $d\rightarrow\infty$ with speed faster
than $\log^{1/2}{n}$. To ensure clustering consistency, we also need to set
the initial cluster centroids in distinct clusters. This is usually violated
in many well known $k$-means methods. For instance, the kmeans function of R
selects the initial centroids by randomness, which cannot guarantee that they
belong to distinct clusters (as $n\rightarrow\infty$). The kmeanspp of R and
the kmeans function of python selects the next centroid with a probability
inversely proportional to the squares of the minimum distance to the previous
selected centroids, which cannot guarantee that they belong to distinct
clusters either. Therefore, it is more important to propose a new $k$-means
method satisfying Condition (b) of Theorem 1.

###  3.3 A New $k$-means Method

We devise an approach to select the $k$ initial centroids by Condition (b) of
Theorem 1. Our approach contains $k$ steps. Each step provides an initial
centroid. In the end, it provides the $k$ initial centroids. In the first
step, we randomly select the first initial centroid $\tilde{\bm{c}}_{1}$. In
the second step, we select the second initial centroid $\tilde{\bm{c}}_{2}$ as
the data point with the maximum distance to the first initial centroid.
Therefore, $\tilde{\bf c}_{2}$ satisfies

$\tilde{\bf c}_{2}=\mathop{\arg\\!\max}_{{\bm{x}}\in{\cal
D},{\bm{x}}\not=\tilde{\bm{c}}_{1}}\|{\bm{x}}-\tilde{\bm{c}}_{1}\|.$ (17)

From the third step (i.e., the $r$th step for $r\geq 3$), we define the
distance measure

$D({\bm{x}})=\min_{\tilde{\bm{c}}\in\tilde{C}_{prev}}(\|{\bm{x}}-\tilde{\bm{c}}\|),$
(18)

where $\tilde{\cal
C}_{prev}=\\{\tilde{\bm{c}}_{1},\dots,\tilde{\bm{c}}_{r-1}\\}$ is the set of
the initial centroids selected by the previous steps. We select the $r$th
initial centroid by

$\tilde{\bf c}_{r}=\mathop{\arg\\!\max}_{{\bm{x}}_{i}\in{\cal
D}\backslash\tilde{C}_{prev}}D({\bm{x}}_{i}).$ (19)

Therefore, $\tilde{\bm{c}}_{r}$ satisfies the max-min principal. It is the
data point which maximizes the minimum distance to the previous selected
initial centroids. In the end, we obtain the set of $k$ initial centroids,
denoted by $\tilde{\cal C}=\\{\tilde{\bm{c}}_{1},\dots,\tilde{\bm{c}}_{k}\\}$,
leading to the proposed $k$-means method by Algorithm 1.

Algorithm 1 Algorithm for the improved $k$-means clustering

1:Input: Data set ${\cal D}=\\{{\bm{x}}_{1},\dots,{\bm{x}}_{n}\\}$ and the
number of clusters $k$

2:Output: labels $1,\dots,k$ for each ${\bm{x}}_{i}\in{\cal D}$.

3:Initialization

4:Choose the first initial centroid $\tilde{\bm{c}}_{1}$ from ${\cal D}$
randomly

5:For each ${\bm{x}}_{i}\in{\cal D}$ distinct from $\tilde{\bm{c}}_{1}$,
compute $\|{\bm{x}}_{i}-\tilde{\bm{c}}_{1}\|$. Choose the second initial
centroid $\tilde{\bm{c}}_{2}$ as the data point with the maximum
$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{1}\|$ value.

6:Starting from the third initial centroid, compute $D({\bm{x}}_{i})$ by (18)
for all ${\bm{x}}_{i}\in{\cal D}$ but distinct from the previous selected
initial centroids. Choose the current initial centroid as the data point with
the maximum $D({\bm{x}}_{i})$ value. Repeat this step until all of the $k$
initial centroids $\tilde{\bm{c}}_{1},\dots,\tilde{\bm{c}}_{k}$ are chosen.

7:Let ${\bm{c}}_{r}=\tilde{\bm{c}}_{r}$ for all $r=1,\dots,k$.

8:Begin Iteration

9:For each ${\bm{x}}_{i}\in{\cal D}$, assign $r$ as the label of
${\bm{x}}_{i}$ if
$\|{\bm{x}}_{i}-{\bm{c}}_{r}\|<\|{\bm{x}}_{i}-{\bm{c}}_{r^{\prime}}\|$ for any
$r^{\prime}\not=r$.

10:For each $r=1,\dots,k$, update ${\bm{c}}_{r}$ by the average of data points
with labels $r$.

11:End Iteration

12:Output

The final solution given by Algorithm 1 depends on the first initial centroid
as it is obtained by randomness. The main difference between Algorithm 1 and
other $k$-means algorithms is the initialization stage given by Steps 1 to 3.
The goal is to make the minimum distance between the initial centroids as
large as possible. The iteration stage given by Steps 5 and 6 is identical to
that used by other $k$-means methods. Similarly, we cannot guarantee that the
solution provided by Algorithm 1 is a global solution of (1). We can only show
that it is a local solution. We show that any local optimizer provided by
Algorithm 1 is clustering consistent, implying that local optimization is good
enough if Algorithm 1 is used.

###### Theorem 2

Any solution provided by Algorithm 1 is a local optimizer of (1).

###### Lemma 1

If Conditions (a), (c), and (d) of Theorem 1 are satisfied, then
$\tilde{\bm{c}}_{1},\dots,\tilde{\bm{c}}_{k}$ reside in distinct clusters in
probability.

###### Theorem 3

If Conditions (a), (c), and (d) of Theorem 1 are satisfied, then any local
optimizer provided by Algorithm 1 for (1) is clustering consistent.

###### Corollary 2

Assume that Condition (a) of Theorem 1 is satisfied, and $p$ is a constant. If
Condition (e) of Corollary 1 holds, then any local optimizer provided by
Algorithm 1 for (1) is clustering consistent.

Because any local optimizer given by Algorithm 1 is clustering consistent, it
does not matter whether the global optimizer of (1) can be numerically
computed in polynomial times or not. It is enough to use a local optimizer to
partition the data set. This means that our research successfully reduces the
$k$-means clustering problem to a local optimization problem, imply that the
NP-hardness issue does not affect the implementation of our method. Therefore,
our method is computationally efficient.

Corollary 2 only requires that the minimum distance between the centroids of
the clusters approaches infinite with the speed faster than $\log^{1/2}{n}$.
It does not matter whether $k$ is bounded or approaches infinite with $n$.
Therefore, Algorithm 1 can be implemented even when there exist many small
clusters. Because outliers can be treated as clusters with sizes equal to $1$,
our method can also be used even when outliers are present. This cannot be
achieved by other $k$-means methods. We put this issue in Section 5.

Although we assume that $k$ is known in Algorithm 1, our method can also be
used even when $k$ is unknown. The idea is to implement Algorithm 1 to various
candidate values of $k$ with the best value of $k$ to be chosen by a
predefined criterion. A few predefined criteria for $k$ have been proposed in
the literature. Examples include the minimum message length (MML) criterion
[12], the minimum description length (MDL) criterion [16], and the Bayesian
information criterion (BIC) [45]. Besides, the Gap statistics is considered as
one of the most popular approaches for the number of clusters [39]. Thus, we
study the Gap statistics approach. We find that it works well in our method
but not in other $k$-means methods. We also put this issue in Section 5.

##  4 Extension

We extend our method to the GMM defined by (7) with arbitrary
${\bm{\Sigma}}_{r}$ and $\pi_{r}$ satisfying $\sum_{r=1}^{k}\pi_{r}=1$.
Besides the $k$-means model, this also includes the linear discriminant
analysis (LDA) model by assuming that all ${\bm{\Sigma}}_{r}$ are identical
and the quadratic discriminant analysis (QDA) model by assuming that
${\bm{\Sigma}}_{r}$ for $r=1,\dots,k$ are different. The GMM has many other
variants. Examples include the spherical and unequal volume model by setting
${\bm{\Sigma}}_{r}=\sigma_{r}^{2}{\bf I}$ with distinct $\sigma_{r}^{2}$, the
diagonal and equal volume and shape model by setting ${\bm{\Sigma}}_{r}$ all
identical and diagonal, and elliptical and equal volume and shape model by
setting ${\bm{\Sigma}}_{r}=\sigma_{r}^{2}{\bm{\Sigma}}_{0}$ for a certain
${\bm{\Sigma}}_{0}$ with distinct $\sigma_{r}^{2}$. These variants have been
incorporated in the mclust package of R.

We show that our method is robust to the specifications of
${\bm{\Sigma}}_{r}$. We do not expect that any $k$-means method is criterion
consistent due to violations of the assumptions of the $k$-means model.
However, we can still show that our method is clustering consistent under a
few weak regularity conditions. This means that criterion consistency and
clustering consistency are two different concepts. Thus, they should be
evaluated separately.

We still focus on the clustering consistency problem of the proposed $k$-means
method. To control the variations of variance-covariance matrices among
clusters, we control the largest eigenvalue of ${\bm{\Sigma}}_{r}$. We show
that the role of the largest eigenvalue used in this section is equivalent to
the role of $\sigma^{2}$ used in the previous section, leading to clustering
consistency of our $k$-means method for the general GMMs.

###### Theorem 4

Assume that ${\bm{x}}_{1},\cdots,{\bm{x}}_{n}$ are independently collected
from (7) and each cluster has its own ${\bm{\Sigma}}_{r}$. If
$\mathop{\lim\\!\inf}_{n\rightarrow\infty}d/(32\lambda_{\max}\log{n})^{1/2}>1$
and $\lim_{n\rightarrow\infty}p\log[d^{2}/(p\lambda_{\max})]/\log{n}=0$, where
$\lambda_{\max}$ is the maximum of the largest eigenvalue of
${\bm{\Sigma}}_{r}$ for $r=1,\dots,k$, then any local solution provided by
Algorithm 1 for (1) is clustering consistent.

###### Corollary 3

Suppose that the assumptions of Theorem 4 are satisfied and $p$ is fixed as
$n\rightarrow\infty$. If (i) $\lambda_{\max}$ is fixed and
$\mathop{\lim\\!\inf}_{n\rightarrow\infty}d/(32\log{n})^{1/2}>\lambda_{\max}$
or (ii) $d$ is fixed and
$\mathop{\lim\\!\sup}_{n\rightarrow\infty}\lambda_{max}\log{n}<d^{2}/32$, then
then any local solution provided by Algorithm 1 for (1) is clustering
consistent.

We have shown that our proposed $k$-mean method is clustering consistent even
when ${\bm{\Sigma}}_{r}$ for $r=1,\dots,k$ are distinct. Thus, our method can
be used to the general GMM defined by (7) for data clustering. This includes
the LDA and QDA models. As we do not specify any conditions for $k$ and
$\pi_{r}$, our method can be implemented even if the cluster sizes are very
unbalanced. This includes the case when some of the clusters are composed by
one observation, meaning that they are outliers. If $k$ is unknown, we also
use the Gap statistic to determine the number of clusters. We evaluate these
issues in our simulation studies.

##  5 Simulation

We compared our $k$-means method with the comparison to our competitors via
Monte Carlo simulations. Our competitors included the EM-algorithm for GMMs
clustering (denoted by EM-GMM), the basic $k$-means clustering (using randomly
selected initial cluster centroids), and the $k$-means++ clustering [1]. To
make initial centroids far away from each other, the $k$-means++ selects the
first initial centroid randomly from the data points. From the second initial
centroid, it defines a distance measure as the minimum distance to the
previous chosen centroids. The probability for the remaining data points to be
selected as a new initial centroid is inversely proportional to the square of
the distance measure. The procedure is repeated until all of the $k$ initial
centroids are chosen. The remaining steps of the $k$-means++ are identical to
those of the basic $k$-means.

Both the basic $k$-means and the $k$-means++ have four variants. The four
variants of the basic $k$-means are the Hartigan-Wong $k$-means [19], the
Lloyd $k$-means [27], the Forgy $k$-means [14], and the MacQueen $k$-means
[29]. They are denoted as HW, Lloyd, Forgy, and MacQueen, respectively. The
four variants of the $k$-means++ are derived by the combinations of those with
the $k$-means++ initialization, respectively. They are denoted as HW++,
Lloyd++, Forgy++, and MacQueen++, respectively. We implemented the EM-GMM by
the Mclust function of the mclust package, the four variants of the basic
$k$-means by the kmeans function of the base library, the four variants of the
$k$-means++ by the kmeanspp function of the LICORS package of R. Therefore, we
had $10$ methods in our simulation. The first was our proposed method. The
remaining nine was our competitors. We compared those methods in the cases
when $k$ was either known or unknown. When $k$ was unknown, we used the Gap
statistics [39] to estimate the number of clusters.

Table 1:  Simulations with $1000$ replications for the percentages of the clustering errors (i.e., $100{\rm CER}$) reported by the proposed $k$-means method with the comparison to its competitors in the cases when outliers are absent ($k=10$) and present ($k=20$ and $k_{outlier}=10$), respectively. |  $k$-means Model
---|---
|  $\phi$ for Outliers Absent |  $\phi$ for Outliers Present
Method | 0.4 | 0.6 | 0.8 | 0.4 | 0.6 | 0.8
Proposed | $6.2$ | $1.1$ | $0.4$ | $6.1$ | $1.4$ | $0.6$
EM-GMM | $12.0$ | $10.7$ | $6.1$ | $36.8$ | $33.7$ | $21.6$
HW | $20.9$ | $20.9$ | $21.7$ | $50.9$ | $49.9$ | $50.3$
Lloyd | $23.1$ | $22.0$ | $22.2$ | $51.8$ | $50.7$ | $50.4$
Forgy | $22.7$ | $22.5$ | $21.8$ | $51.9$ | $50.6$ | $50.3$
MacQueen | $21.3$ | $21.5$ | $22.2$ | $51.3$ | $50.4$ | $50.8$
HW++ | $16.0$ | $15.4$ | $15.3$ | $44.8$ | $45.1$ | $45.4$
Lloyd++ | $16.5$ | $15.4$ | $15.8$ | $45.0$ | $44.9$ | $46.0$
Lloyd++ | $16.6$ | $16.0$ | $15.8$ | $45.6$ | $45.3$ | $45.8$
MacQueen++ | $16.3$ | $16.2$ | $15.7$ | $45.3$ | $45.0$ | $46.0$
| QDA Model
Method |  $\phi$ for Outliers Absent |  $\phi$ for Outliers Present
Proposed | $10.4$ | $3.3$ | $0.8$ | $10.4$ | $3.0$ | $1.0$
EM-GMM | $20.1$ | $16.6$ | $11.0$ | $48.4$ | $45.4$ | $41.8$
HW | $23.1$ | $21.4$ | $21.6$ | $50.4$ | $48.5$ | $48.5$
Lloyd | $25.2$ | $22.8$ | $22.2$ | $51.6$ | $49.5$ | $49.3$
Forgy | $25.2$ | $22.9$ | $22.3$ | $51.0$ | $49.4$ | $49.2$
MacQueen | $23.4$ | $22.4$ | $21.9$ | $50.7$ | $49.0$ | $49.0$
HW++ | $19.0$ | $16.0$ | $15.7$ | $43.8$ | $42.4$ | $42.6$
Lloyd++ | $19.8$ | $16.3$ | $15.9$ | $43.1$ | $42.9$ | $43.2$
Forgy++ | $19.7$ | $16.2$ | $15.7$ | $43.3$ | $42.7$ | $43.4$
MacQueen++ | $19.8$ | $16.7$ | $16.0$ | $43.4$ | $43.0$ | $43.2$

We used the GMM defined by (7) with $p=5$ to generate clusters. We considered
two scenarios. In the first scenario, we chose
${\bm{\Sigma}}_{r}=\sigma^{2}{\bf I}$ for all $r=1,\dots,k$, implying that it
was the $k$-means model. In the second scenario, we chose all
${\bm{\Sigma}}_{r}$ different, implying that it was the QDA model. In both
scenarios, we compared the ten methods under the cases when outliers were
absent and present, respectively.

We carried out a three-step procedure to generate the clusters. In the first
step, we generated the cluster expected vectors ${\bm{\mu}}_{r}$ independently
by ${\cal N}({\bm{0}},\phi^{2}{\bf I})$. In the second step, we generated the
clusters sizes $n_{r}$. When the outliers were absent, we chose $k=10$. When
the outliers were present, we chose $k=20$. In both cases, we generated $5$
small clusters and $5$ large clusters. If outliers were present, then the
number of outliers was $k_{outlier}=10$. The sizes of the small clusters were
independently generated by $n_{r}~{}\sim{\cal P}(50)$. The sizes of the large
clusters were independently generated by $n_{r}\sim{\cal P}(1000)$. The
cluster sizes for outliers were $n_{r}=1$. In the third step, we generated the
observations. In the $k$-means model, we independently generated $n_{r}$
observations from ${\cal N}({\bm{\mu}}_{r},0.1^{2}{\bf I})$. In the QDA model,
for the $r$th cluster, we first independently generated $5$ eigenvalues
$d_{r1},\dots,d_{r5}$ of ${\bm{\Sigma}}_{r}$ by $d_{rj}\sim{\cal U}(0,0.2)$
for $j=1,\dots,5$. We then independently generated $5$ eigenvectors of
${\bm{\Sigma}}_{r}$ by a random orthogonal matrix ${\bf U}_{r}$. We calculated
the variance-covariance matrix of the $r$th cluster by ${\bm{\Sigma}}_{r}={\bf
U}_{r}{\bf D}_{r}{\bf U}_{r}^{\top}$ with ${\bf D}_{r}={\rm
diag}(d_{r1},\dots,d_{r_{5}})$. After all ${\bm{\Sigma}}_{r}$ were derived, we
independently generated $n_{r}$ observations of the $r$th cluster by ${\cal
N}({\bm{\mu}}_{r},{\bm{\Sigma}}_{r})$. The total number of observations was
$n=\sum_{r=1}^{k}n_{r}$.

We used the percentages of cluster errors to compare our proposed method with
our competitors. Due to the permutation issue, the clustering labels are
equivalent if the result of one method can be changed to that of the other
method by a permutation $\tau$ on $\\{1,\dots,k\\}$. Using this criterion, we
compute the percentage of clustering errors. In particular, let
$\\{\hat{n}_{1},\dots,\hat{n}_{k}\\}$ be the sizes of clusters given by a
clustering method. The estimated labels may not match the true labels. To
address the issue, we fix the true labels and change the estimated labels by a
permutation, leading to the clustering error rate as

${\rm CER}=1-{1\over n}\max_{\tau}\sum_{r=1}^{k}\hat{n}_{r,\tau(r)},$ (20)

where $\hat{n}_{r,\tau(r)}$ is the number of data points with label $\tau(r)$
assigned by the clustering method and label $r$ assigned by the truth. The
percentage of the clustering errors is $100{\rm CER}$. The computation for
${\rm CER}$ can be easily solved by the well known Hungarian algorithm.

We generated $1,\\!000$ data sets from (7) for each selected $\phi$ in the
$k$-means and the QDA models, respectively. We applied all of the ten methods
to each data set. We computed the percentages of clustering errors (Table 1).
We compared the performance of the ten methods under the $k$-means and the QDA
models without and with the outliers, respectively. We found that our method
was consistently and significantly better than our competitors. In the
comparison between our competitors, we found that the EM-GMM method was better
than the four basic $k$-means and the four $k$-means++ methods. There were no
significantly differences within the four basic $k$-means and the four
$k$-means++ methods, respectively. The four $k$-means++ methods were better
than the four basic $k$-means methods. This supports the previous statement
that the initialization adopted by the $k$-means++ improves $k$-means
clustering. Because the $k$-means++ was significantly worse than our proposed
method, we conclude that the initialization adopted by the $k$-means++ is far
from the optimal. The proposed method significantly improves the quality of
$k$-means clustering.

We then compared the performance of each of ten methods between the $k$-means
and the QDA models in the cases when outliers were absent and present,
respectively. We found that the performance of our method was not
significantly affected by the presence of the outliers, but the performance of
our competitor was. This means that our method is robust to outliers, but our
competitors are not. In the QDA model, the performance of our method was still
better than that of the EM-GMM method, implying that the proposed $k$-means
method is also appropriate for general GMMs.

Table 2:  Simulations with $1000$ replications for the average number of outliers reported by the proposed $k$-means method with the comparison to its competitors, where $k=20$ and $k_{outlier}=10$. |  $\phi$ for $k$-means Model |  $\phi$ for QDA Model
---|---|---
Method | $0.4$ | $0.6$ | $0.8$ | $0.4$ | $0.6$ | $0.8$
Proposed | $9.94$ | $9.93$ | $9.94$ | $9.97$ | $9.97$ | $9.98$
EM-GMM | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
HW | $0.07$ | $0.06$ | $0.10$ | $0.20$ | $0.25$ | $0.27$
Lloyd | $0.17$ | $0.16$ | $0.16$ | $0.31$ | $0.37$ | $0.35$
Forgy | $0.14$ | $0.17$ | $0.16$ | $0.31$ | $0.35$ | $0.34$
MacQueen | $0.21$ | $0.19$ | $0.18$ | $0.45$ | $0.44$ | $0.44$
HW++ | $0.38$ | $0.23$ | $0.18$ | $0.83$ | $0.63$ | $0.53$
Lloyd++ | $0.81$ | $0.53$ | $0.36$ | $1.67$ | $1.15$ | $0.89$
Forgy++ | $0.84$ | $0.50$ | $0.32$ | $1.63$ | $1.16$ | $0.85$
MacQueen++ | $0.85$ | $0.53$ | $0.37$ | $1.68$ | $1.17$ | $0.88$

We studied whether the ten methods could be used to identify outliers when
they were present. We claimed a cluster as an outlier if its size was one. The
simulation results showed that only our method successfully identified all of
the outliers. All of our competitors failed to identify most of the outliers.
In all of the cases that we simulated, the EM-GMM did not discover any
outliers. In most of the cases, the basic $k$-means the $k$-means++ methods
identified zero or one outlier. Therefore, we conclude that the proposed
method can identify outliers but our competitors cannot.

Figure 1: Percentage of correctly identified number of clusters (i.e., % of
$\hat{k}=k$) reported by the combination of the Gap statistics with the
proposed $k$-means, the EM-GMM, the HW, and the HW++ methods, respectively.

In the end, we compared the performance of our method with that of our
competitors in the case when $k$ was unknown. In practice, the determination
of the number of clusters is considered as one of the most important issues
for the $k$-means. Many approaches have been proposed in the literature.
Examples include the MML criteria [12], the MDL criterion [16], the BIC [45],
and the Bayesian non-parametric [11, 34] approaches. Among those, the Gap
statistics [39] is considered as one of the most popular approaches.
Therefore, we decided to carry out the comparison based on the Gap statistics.
In particular, we estimated the number of clusters by the combination of each
of the ten methods with the Gap statistic. We denoted $\hat{k}$ as the
estimated number of clusters and $k$ as the true number of clusters. We
concluded that the true number of clusters was correctly identified if
$\hat{k}=k$ or misidentified otherwise. We compared the ten methods using the
percentage of $\hat{k}=k$ (Figure 1). Because there was no significant
difference among the four basic $k$-means and the four $k$-means++ methods, we
only report the results for HW and HW++. Our simulations showed that both the
proposed $k$-means and the EM-GMM colud identify the true number of clusters
if they were combined with the Gap statistics. None of the basic $k$-means or
$k$-means++ methods could identify the number of clusters.

In summary, our simulations show that our proposed $k$-means method
outperforms our competitors in all of the aspects that we have studied. It has
low clustering error rates. It can also detect outliers. Our method can
identify the number of clusters in the case when the true $k$ is unknown. Our
study generally supports the Gap statistics approach for the determination of
the number of clusters. It strongly suggests not to implement the Gap
statistics approach to the basic $k$-means or $k$-means++ methods for the
determination of the number of clusters.

##  6 Application

We implemented our method with the comparison to our competitors to the Swiss
Banknote data set. The data set was originally provided by Tables 1 and 2 on
Page 5–8 of [13]. It is now available in the mclust package of R. The Swiss
Banknote data set contains six measurements (denoted by $x_{1},\dots,x_{6}$,
respectively) made on $100$ genuine and $100$ counterfeit old-Swiss 1000-franc
bank notes of the Swiss National Bank in 1980s. Thus, $n=200$. It also
contains an additional status variable for the status of the banknotes as
genuine or counterfeit. Thus, status is the ground truth variable. It was used
in the computation of CER.

Figure 2:  The original observed cluster pattern by plots (a) $x_{4}$ vs
$x_{6}$ and (b) $x_{5}$ vs $x_{6}$, and the generated randomly stamped cluster
pattern displayed by plots (c) $x_{4}$ vs $x_{6}$ and (d) $x_{5}$ vs $x_{6}$.

We studied two circumstances. In the first circumstance, we implemented all of
the ten methods directly to the observed data. We computed the percentage of
cluster errors (i.e., $100{\rm CER}$) for each method using (20), where we
assumed that $k=2$ was known. The percentage of the proposed method, the basic
HW, Loyod, Forgy, and MacQueen methods, and the HW++, Loyod++, Forgy++, and
MacQueen++ methods were all $0\%$. The percentage of the EM-GMM was $0.5\%$.
Therefore, we conclude that the performances of all of the methods were close
in the direct implementation to the Swiss Banknote data set.

To emphasize the importance of clustering consistency, we studied the second
circumstance. We used the random stamps approach. The random stamps approach
treats the observed pattern as the base stamp. The base stamp is duplicated
first and them shifted randomly by a distribution. In the end, the random
stamps approach generates a pattern by duplicating the base stamp with a
number of location shifts. The goal for us to consider the random stamps
approach was to investigate the influence of well separated cluster patterns
on a selected clustering method.

To implement the random stamps approach, we first selected positive integer
$m$ for the number of stamps. We then generated $m$ vectors
${\bm{\xi}}_{v}=(\xi_{v1},\dots,\xi_{v6})$ identically and independently from
the multivariate normal distribution with mean zero and variance-covariance
matrix ${\bf A}={\rm diag}(10s_{1},\dots,10s_{6})$, where $s_{v}$ was the
standard error of $x_{j}$ for $j=1,\dots,6$ provided by the base stamp. We
used ${\bm{\xi}}_{v}$ to shift the base stamp pattern individually. The goal
was to ensure that all of the stamps were well separated by the location
shifts. To make it clearer, we illustrate a realization of the random stamps
approach by Figure 2 with $m=5$ stamps. Note that the Swiss Banknote data set
contains $6$ variables. We need to generate a six-dimensional vector to shift
the base stamp each time. Based on the original observed cluster patterns
displayed by the plot using $x_{4}$ and $x_{6}$ given by Figure 2(a), we
shifted $x_{i4}$ to $x_{i4}+\xi_{v4}$ and $x_{i6}$ to $x_{i6}+\xi_{v6}$,
$i=1,\dots,n$. We obtained Figure 2(c). Similarly, we obtained Figure 2(d)
using Figure 2(b). Figures 2(c) and 2(d) were obtained from the same data set
using different variables for visualization. Thus, the number of clusters was
$k=2m=10$. The stamps were well separated. Each pair of clusters contained by
a stamp was identical to the observed cluster pattern contained by the base
stamp. Therefore, the performance of a nice clustering method should not be
significantly influenced by the random stamps approach used in the
application.

Figure 3:  Percentage (%) of cluster errors for the generated patterns derived
by randomly stamping the Swiss Banknote data set, where the number of clusters
is $k=2m$.

We varied $m$ from $1$ to $15$ and obtained $k=2$ to $k=30$ clusters,
respectively, by the random stamps approach for the Swiss Banknote data. We
implemented the proposed $k$-means, the four basic $k$-means, the four
$k$-means++, and the EM-GMM methods. We compared these methods using the
percentage of clustering errors (Figure 3). As the shifts were randomly
generated, we repeated $1000$ times for each selected $m$. As the curves for
the four basic $k$-means and the four $k$-means methods were close,
respectively, we only put the curves for the HW basic $k$-means and HW++
$k$-means++ methods in Figure 3. We found that only the percentages of
clustering errors reported by the proposed $k$-means method were consistently
low. This means that our competitors are not reliable to the increase of the
number of clusters even if they are far away from each other. This problem has
been overcome by our method. Therefore, we conclude that the performance of
the proposed $k$-means method is better than our competitors in the
application to the Swiss Banknote data set.

##  7 Discussion

We study asymptotic properties for the $k$-means clustering. We point out that
the proposed clustering consistency is more appropriate than the previous
criterion consistency for the $k$-means, when it is applied to the recent
unsupervised learning problems arisen in statistics, machine learning, and
artificial intelligence. Using the proposed clustering consistency, we
investigate existing $k$-means methods, including the basic $k$-means and the
$k$-means++ adopted by many software packages, such as R and python. We find
that all of the existing $k$-means methods are not optimal. To ensure a
$k$-means method to be clustering consistent, it is more appropriate to choose
the initial $k$-centroids as far as possible, leading to the proposed
$k$-means method. We show that any local solution provided by the proposed
$k$-means method is clustering consistent, implying that it is better than the
existing $k$-means methods adopted by R and python. More interesting, we find
that the combination of the $k$-means with the well known Gap statistics can
be used to determine the number of clusters when $k$ is unknown, which cannot
be achieved by the basic $k$-means and the $k$-means++ methods.

The performance of the basic $k$-means and the $k$-means++ becomes worse as
$k$ becomes large. This can be interpreted by Theorem 1 of the article. For
any $k$-means method, if it can guarantee that the centroids reside in
distinct clusters sooner or later in the iterations, then the $k$-means method
is clustering consistent. If a basic $k$-means or a $k$-means++ method is used
to the case when $k$ is small (e.g., $k=2$ or $3$), then the centroids may
belong to distinct clusters after a few steps of iterations. This can be
hardly achieved in the case when $k$ is large (e.g., $k=10$ or higher).
Therefore, it is more important to evaluate a clustering method under the case
when $k$ is not very small. In addition, we find that the performance of the
proposed $k$-means is better than that of the EM-GMM in all of the cases that
we have studied, implying that the proposed $k$-means is comparable to the GMM
method. If a basic $k$-means or a $k$-means++ method is used, then we cannot
draw a similar conclusion. Therefore, the derivation of an optimal $k$-means
method is critical in the comparison between the $k$-means and its
competitors.

We treat clustering consistency as the minimum requirement of a clustering
method. It means that if the distances between individual clusters increases,
then the clustering error rates should be small. We show that clustering
consistency can be achieved by local solutions of our method, indicating that
the NP-hardness is not an issue in our method. In practice, it is enough to
use our method to provide a local optimizer. Note that clustering consistency
is an asymptotic property and the sample size is always finite in practice. We
may still gain benefits if we want to determine a better solution by
replications. However, this does not affect the asymptotic properties
formulated under clustering consistency proposed in the current research.

##  Appendix A Proofs

Proof of Theorem 1. Let $\tilde{\bm{c}}_{1},\dots,\tilde{\bm{c}}_{k}$ be the
initial centroids of the clusters selected by the $k$-means. Without the loss
of generality, we can assume that their true cluster labels are $1,\cdots,k$,
respectively. In the first iteration, the $k$-means assigns ${\bm{x}}_{i}$ to
the $r$th cluster if
$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r}\|<\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|$
for any $r^{\prime}\not=r$. Let $\check{r}_{i}$ be the true cluster label of
${\bm{x}}_{i}$. By the implmentation of the triangle inequality to
$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r}\|=\|({\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}})+({\bm{\mu}}_{\check{r}_{i}}-{\bm{\mu}}_{r})+({\bm{\mu}}_{r}-\tilde{\bm{c}}_{r})\|$,
we have

$\displaystyle{\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r}\|\geq}$ | $\displaystyle{{}\|{\bm{\mu}}_{\check{r}_{i}}-{\bm{\mu}}_{r}\|-\|({\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}})+({\bm{\mu}}_{r}-\tilde{\bm{c}}_{r})\|}$
---|---
$\displaystyle{\geq}$ | $\displaystyle{{}\|{\bm{\mu}}_{\check{r}_{i}}-{\bm{\mu}}_{r}\|-(\|{\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}}\|+\|\tilde{\bm{c}}_{r}-{\bm{\mu}}_{r}\|).}$

Similarly, by
$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{\check{r}_{i}}\|=\|({\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}})+({\bm{\mu}}_{\check{r}_{i}}-\tilde{\bm{c}}_{\check{r}_{i}})\|$,
we have

$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{\check{r}_{i}}\|\leq\|{\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}}\|+\|\tilde{\bm{c}}_{\check{r}_{i}}-{\bm{\mu}}_{\check{r}_{i}}\|.$

If $r\not=\check{r}_{i}$, then
$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r}\|<\|{\bm{x}}_{i}-\tilde{\bm{c}}_{\check{r}_{i}}\|$,
leading to

$\displaystyle{\|{\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}}\|\geq}$ | $\displaystyle{{}{1\over 2}\|{\bm{\mu}}_{\check{r}_{i}}-{\bm{\mu}}_{r}\|-{1\over 2}(\|\tilde{\bm{c}}_{\check{r}_{i}}-{\bm{\mu}}_{\check{r}_{i}}\|+\|\tilde{\bm{c}}_{r}-{\bm{\mu}}_{r}\|)}$
---|---
$\displaystyle{\geq}$ | $\displaystyle{{}{d\over 2}-{1\over 2}(\|\tilde{\bm{c}}_{\check{r}_{i}}-{\bm{\mu}}_{\check{r}_{i}}\|+\|\tilde{\bm{c}}_{r}-{\bm{\mu}}_{r}\|).}$

It is derived by the combination of the above two inequalities with the
definition of $d$ given by (16). Note that all of
$\|{\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}}\|^{2}$,
$\|\tilde{\bm{c}}_{\check{r}_{i}}-{\bm{\mu}}_{\check{r}_{i}}\|^{2}$, and
$\|\tilde{\bm{c}}_{r}-{\bm{\mu}}_{r}\|^{2}$ are $\sigma^{2}\chi_{p}^{2}$
random variables. We can use the formulation for the upper tail probability of
the $\chi_{q}^{2}$-distribution given by [22] as

${\rm Pr}\\{\chi_{q}^{2}\geq
u\\}\leq{u\over\sqrt{\pi}(u-q+2)}\exp\left\\{-{1\over
2}[u-q-(q-2)\log(u/q)+\log{q}]\right\\},$ (21)

for any $q\geq 2$ and $u\geq q-2$. If $p\geq 2$, then

| $\displaystyle{{}{\rm Pr}({\bm{x}}_{i}\not\in C_{r,T}\ \exists
i\in\\{1,\dots,n\\})}$
---|---
$\displaystyle{\leq}$ | $\displaystyle{{}\sum_{i=1}^{n}{\rm Pr}({\bm{x}}_{i}\not\in C_{r,T})}$
$\displaystyle{=}$ | $\displaystyle{{}\sum_{i=1}^{n}{\rm Pr}\left\\{\|{\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}}\|+{1\over 2}\|\tilde{\bm{c}}_{\check{r}_{i}}-{\bm{\mu}}_{\check{r}_{i}}\|+{1\over 2}\|\tilde{\bm{c}}_{r}-{\bm{\mu}}_{r}\|\geq{d\over 2}\right\\}}$
$\displaystyle{\leq}$ | $\displaystyle{{}\sum_{i=1}^{n}\left\\{{\rm Pr}\left(\|{\bm{x}}_{i}-{\bm{\mu}}_{\check{r}_{i}}\|\geq{d\over 4}\right)+{\rm Pr}\left(\|\tilde{\bm{c}}_{\check{r}_{i}}-{\bm{\mu}}_{\check{r}_{i}}\|\geq{d\over 4}\right)+{\rm Pr}\left(\|\tilde{\bm{c}}_{r}-{\bm{\mu}}_{r}\|\geq{d\over 4}\right)\right\\}}$
$\displaystyle{\leq}$ | $\displaystyle{{}3n{\rm Pr}(\chi_{p}^{2}\geq{d^{2}\over 16\sigma^{2}}).}$

By taking $q=p$ and $u=d^{2}/(16\sigma^{2})$ in (21), we have

| $\displaystyle{{}{\rm Pr}({\bm{x}}_{i}\not\in C_{r,T}\ \exists
i\in\\{1,\dots,n\\})}$
---|---
$\displaystyle{\leq}$ | $\displaystyle{{}{3d^{2}\over\sqrt{\pi}[d^{2}-16\sigma^{2}(p-2)]}\exp\left\\{-{1\over 2}\left[{d^{2}\over 16\sigma^{2}}-p-(p-2)\log{d^{2}\over 16p\sigma^{2}}+\log{p}-2\log{n}\right]\right\\}.}$
(22)

Using Condition
$\lim_{n\rightarrow\infty}p\log[d^{2}/(p\sigma^{2})]/\log{n}=0$, we ignore the
second, third, and fourth terms in the exponent when we evaluate the limit of
the above probability. We only need to compare the first and the last terms.
We find that the probability goes to $0$ if
$\lim_{n\rightarrow\infty}[d^{2}/(16\sigma^{2})-2\log{n}]=\infty$, which can
be achieved by Condition
$\lim\\!\inf_{n\rightarrow\infty}d/[4\sigma(2\log{n})^{1/2}]>1$ of the
theorem. If $p=1$, then by the same approach, we have

${\rm Pr}({\bm{x}}_{i}\not\in C_{r,T}\ \exists i\in\\{1,\dots,n\\})\leq 3n{\rm
Pr}(\chi_{1}^{2}\geq{d^{2}\over 16\sigma^{2}}).$ (23)

By

$\lim_{t\rightarrow\infty}{{\rm Pr}[{\cal N}(0,1)\geq
t]\over({\sqrt{2\pi}t})^{-1}e^{-{t^{2}\over 2}}}=1$

and the connection between the $\chi_{1}^{2}$ and ${\cal N}(0,1)$
distributions, we have

${\rm Pr}\left(\chi_{1}^{2}\geq{d^{2}\over 16\sigma^{2}}\right)\leq
e^{-{d^{2}\over 32\sigma^{2}}}$

for sufficiently large $d^{2}/\sigma^{2}$. Then, we obtain

${\rm Pr}({\bm{x}}_{i}\not\in C_{r,T}\ \exists i\in\\{1,\dots,n\\})\leq
e^{-({d^{2}\over 32\sigma^{2}}-\log{n})}.$

It also approaches $0$ as $n\rightarrow\infty$ under the assumptions of the
theorem. Asymptotically, $\check{r}_{i}$ and $r$ must be identical for all
$i$. We draw the conclusion. $\diamondsuit$

Proof of Corollary 1. By treating $p$ as a constant in the evaluation of the
asymptotic properties of (22) and (23) as $n\rightarrow\infty$, we can show
that the $k$-means method is clustering consistent if
$\lim\\!\inf_{n\rightarrow\infty}d^{2}/(32\sigma^{2}\log{n})>1$, which is
satisfied when either (i) or (ii) of (e) holds. $\diamondsuit$

Proof of Theorem 2. It is enough to show that the SSQ criterion always
decreases in Steps 5 and 6, respectively, in Algorithm 1. In Step 5, each
${\bm{x}}_{i}$ is assigned a label by lowest $\|{\bm{x}}_{i}-{\bm{c}}_{r}\|$
value. The $SSQ(\widehat{\mathscr{C}}_{k})$ value does not increase. In Step
6, each ${\bm{c}}_{r}$ is the least square estimator of the cluster centroid
based on the given label. The $SSQ(\widehat{\mathscr{C}}_{k})$ value does not
increase either. Therefore, any solution provided by Algorithm 1 is a local
optimizer of (1). $\diamondsuit$

Proof of Lemma 1. Let $\tilde{\cal
C}_{prev}=\\{\tilde{\bm{c}}_{1},\cdots,\tilde{\bm{c}}_{r-1}\\}$ be the initial
centroids selected by the previous steps and assume that the current step is
to select $\tilde{\bm{c}}_{r}$. For distinct
${\bm{x}}_{i},{\bm{x}}_{i^{\prime}}\in{\cal D}\backslash\tilde{\cal C}_{prev}$
satisfying that $r$ does not reside in the clusters containing $\tilde{\cal
C}_{prev}$ but $r^{\prime}$ does, where the true labels of ${\bm{x}}_{i}$ and
${\bm{x}}_{i^{\prime}}$ are $r$ and $r^{\prime}$, respectively, we implement
the triangle inequality to
$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|=\|({\bm{x}}_{i}-{\bm{\mu}}_{r})+({\bm{\mu}}_{r}-{\bm{\mu}}_{r^{\prime}})+({\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}})\|$
and obtain

$\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|\geq\|{\bm{\mu}}_{r}-{\bm{\mu}}_{r^{\prime}}\|-(\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|).$

We implement the triangle inequality to
$\|{\bm{x}}_{i^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|=\|({\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}})+({\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}})\|$
and obtain

$\|{\bm{x}}_{i^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|\leq\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|.$

Combining the above two inequalities, we obtain

$\displaystyle{\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|-\|{\bm{x}}_{i^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|\geq}$ | $\displaystyle{{}\|{\bm{\mu}}_{r}-{\bm{\mu}}_{r^{\prime}}\|-(\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|}$
---|---
| $\displaystyle{{}\hskip
56.9055pt+\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|),}$

leading to

$\displaystyle{\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|-\|{\bm{x}}_{i^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|\geq}$ | $\displaystyle{{}\min_{r\not=r^{\prime}}\|{\bm{\mu}}_{r}-{\bm{\mu}}_{r^{\prime}}\|-\max\\{(\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|}$
---|---
| $\displaystyle{{}\hskip
56.9055pt+\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|)\\}}$
$\displaystyle{=}$ | $\displaystyle{{}d-\max\\{(\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|}$
| $\displaystyle{{}\hskip
28.45274pt+\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|)\\}.}$

Because $\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|^{2}$,
$\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|^{2}$,
$\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|^{2}$, and
$\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|^{2}$ are
$\sigma^{2}\chi_{p}^{2}$ random variables, we have

| $\displaystyle{{}{\rm
Pr}\left[\min(\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|-\|{\bm{x}}_{i^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|)>0\right]}$
---|---
$\displaystyle{\geq}$ | $\displaystyle{{}1-{\rm Pr}\left[\max(\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|+\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|)>d\right]}$
$\displaystyle{\geq}$ | $\displaystyle{{}1-\sum_{i=1}^{n}{\rm Pr}\left[(\|{\bm{x}}_{i}-{\bm{\mu}}_{r}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|+\|{\bm{x}}_{i^{\prime}}-{\bm{\mu}}_{r^{\prime}}\|+\|{\bm{\mu}}_{r^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|)>d\right]}$
$\displaystyle{\geq}$ | $\displaystyle{{}1-4n{\rm Pr}(\chi_{p}^{2}>{d^{2}\over 16\sigma^{2}}).}$

If $p\geq 2$, then by (21), we obtain

| $\displaystyle{{}{\rm
Pr}\left[\min(\|{\bm{x}}_{i}-\tilde{\bm{c}}_{r^{\prime}}\|-\|{\bm{x}}_{i^{\prime}}-\tilde{\bm{c}}_{r^{\prime}}\|)>0\right]}$
---|---
$\displaystyle{\geq}$ | $\displaystyle{{}1-{4d^{2}\over\sqrt{\pi}[d^{2}-16\sigma^{2}(p-2)]}\exp\left\\{-{1\over 2}\left[{d^{2}\over 16\sigma^{2}}-p-(p-1)\log{d^{2}\over 16p\sigma^{2}}+\log{p}-2\log{n}\right]\right\\},}$
(24)

which goes to $1$ as $n\rightarrow\infty$. This means that
${\bm{x}}_{r^{\prime}}$ cannot be the $r$ initial centroid asymptotically.
Note that ${\bm{x}}_{i}$ and ${\bm{x}}_{i^{\prime}}$ represent general points
satisfying the condition. We draw the conclusion for the case when $p\geq 2$.
Similarly, we draw the conclusion for the case when $p=1$. This completes the
proof. $\diamondsuit$.

Proof of Theorem 3. The conclusion can be directly implied by the Combination
of Theorem 1 and Lemma 1. $\diamondsuit$

Proof of Corollary 2. The conclusion can be directly implied by the
Combination of Corollary 1 and Lemma 1. $\diamondsuit$

Proof of Theorem 4. We almost repeat the entire proofs of Theorems 1 and Lemma
1 except the distribution of $\|{\bm{x}}_{i}-{\bm{\mu}}_{r_{i}}\|^{2}$ when
the true label of ${\bm{x}}_{i}$ is $r_{i}$. As ${\bm{\Sigma}}_{r_{i}}$ is a
general variance-covariance matrix of a multivariate normal distribution, the
distribution of $\|{\bm{x}}_{i}-{\bm{\mu}}_{r_{i}}\|^{2}$ cannot be expressed
by a constant times a chi-square distribution. We only have
$({\bm{x}}_{i}-{\bm{\mu}}_{r_{i}})^{\top}{\bm{\Sigma}}_{r_{i}}^{-1}({\bm{x}}_{i}-{\bm{\mu}}_{r_{i}})=\|{\bm{z}}_{i}\|^{2}\sim\chi_{p}^{2}$
for each $i\in\\{1,\dots,n\\}$, where
${\bm{z}}_{i}={\bm{\Sigma}}_{r}^{-1/2}({\bm{x}}_{i}-{\bm{\mu}}_{r_{i}})$. By
${\bm{x}}_{i}-{\bm{\mu}}_{r_{i}}={\bm{\Sigma}}_{r_{i}}^{1/2}{\bm{z}}_{i}$, we
have
$\|{\bm{x}}_{i}-{\bm{\mu}}_{r_{i}}\|\leq\lambda_{\max}^{1/2}\|{\bm{z}}_{i}\|$,
leading to
$\|{\bm{x}}_{i}-{\bm{\mu}}_{r_{i}}\|^{2}\leq\lambda_{\max}\|{\bm{z}}_{i}\|^{2}$,
implying that ${\rm Pr}(\|{\bm{x}}_{i}-{\bm{\mu}}_{r_{i}}\|^{2}\leq t)\leq{\rm
Pr}(\lambda_{\max}\chi_{q}^{2}\leq t)$ for any $t>0$. We next obtain two
inequalities similar to (22) and (24) using $\lambda_{\max}$ to replace
$\sigma^{2}$, respectively. Then, we draw the conclusion for the case when
$p\geq 2$. The proof for the case when $p=1$ is similar. By the combination of
the two cases, we draw the conclusion. $\diamondsuit$

Proof of Corollary 3. The conclusion can be directly implied by Theorem 4.
$\diamondsuit$

## References

  * [1] Arthur, D.; Vassilvitskii, S. (2007). k-means++: the advantages of careful seeding. Proceedings of the 18th annual ACM-SIAM symposium on Discrete algorithms. Society for Industrial and Applied Mathematics Philadelphia, PA, USA. pp. 1027–1035.
  * [2] Bock, H. (2008). Origins and extensions of the k-means algorithm in cluster analysis electron. Electronic Journal for History of Probability and Statistics, 4, Article 14.
  * [3] Cardot, H., Cénac, P., and Monnez, J.M. (2012). A fast and recursive algorithm for clustering large datasets with $k$-medians. Computational Statistics and Data Analysis, 6, 1434-1449.
  * [4] Chakraborty, S. and Das, S. (2022). On the strong consistency of kernel k‐means: A Rademacher complexity approach. Statistics and Probability Letters, 182, 109291.
  * [5] Chaturvedi, A., Green, P.E., and Caroll, J.D., (2001). K-modes clustering. Journal of Classification, 18, 35-55.
  * [6] Delaigle, A., Hall, P., Pham, T. (2019). Clustering functional data into groups by using projections. Journal of Royal Statistical Society Series B, 81, 271-304.
  * [7] Dempster, A., Laird, N., Rubin, D. (1997). Maximum likelihood from incomplete data via the EM algorithm (with discussion). Journal of Royal Statistical Society Series B, 39, 1-38.
  * [8] Du, Q., and Wong, T.W. (2002). Numerical studies for MacQueen’s $k$-means algorithms for computing the centroidal Voronoi tessellations. Computer and Mathematics with Applications, 44, 511-523.
  * [9] Ester, M., Kriegel, H.P., Sander, J., and Xu, X (1996). A density-based algorithm for discovering clusters in large spatial databases with noise. In kdd, 96, 226-231.
  * [10] Falkenauer, E., and Marchand, A. Using k-means? consider arrayminer. In International Conference on Mathematics and Engineering Techniques in Medicine and Biological Sciences, Las Vegas. 2001.
  * [11] Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annals of Statistics, 1, 209–230.
  * [12] Figueiredo, M. A. T., and Jain, A. K. (2002). Unsupervised learning of finite mixture models. IEEE Transactions on pattern analysis and machine intelligence, 24, 381-396.
  * [13] Flury, B. and Riedwyl, H. (1988). Multivariate Statistics: A practical approach. London: Chapman & Hall, Tables 1.1 and 1.2, pp. 5-8.
  * [14] Forgy, E. W. (1965). Cluster analysis of multivariate data: efficiency vs interpretability of classifications. Biometrics, 21, 768–769
  * [15] Fraley, C. and Raftery, A.E. (2002). Model-based clustering, discriminant analysis, and density estimation. Journal of the American Statistical Association, 97, 611-631.
  * [16] Hansen, M. H., and Yu, B. (2001). Model selection and the principle of minimum description length. Journal of the American Statistical Association, 96, 746-774.
  * [17] Hartigan, J.A. (1975). Clustering Algorithms. New York, Wiley.
  * [18] Hartigan, J.A. (1978). Asymptotic distributions for clustering criteria. Annals of Statistics, 6, 117-131.
  * [19] Hartigan, J. A. and Wong, M. A. (1979). A $k$-means clustering algorithm. Applied Statistics, 28, 100–108.
  * [20] Henning, C. (2015). What are the true clustering? Pattern Recognition Letters, 64, 53-62.
  * [21] Huang, J.Z., Ng, M.K., Rong, H, and Li, Z.(2005). Automated variable weighting in k-means type clustering. IEEE transactions on pattern analysis and machine intelligence, 5, 657-668.
  * [22] Inglot, T. (2010). Inequalities for quantiles of the chi-square distribution. Probability and Mathematical Statistics, 30, 339-351.
  * [23] Jain, A.K. (2010). Data clustering: 50 years beyond K-means. Pattern Recognition Letters, 31, 651-666.
  * [24] Jing, L., Ng, M.K., and Huang, J. (2007). An entropy weighting k-means algorithm for subspace clustering of high-dimensional sparse data. IEEE Transactions on knowledge and data engineering, 19, 1026-1041.
  * [25] Kriegel, H.P., Kröger, P., Sander, J., and Zimek, A. (2001). Density-based clustering. WIRES Data Mining Knowledge Discovery, 1, 231-240.
  * [26] Lau, J.W. and Green, P.J. (2007). Bayesian model-based clustering procedures. Journal of Computational and Graphical Statistics, 16, 526-558.
  * [27] Lloyd, S. P. (1982). Least squares quantization in PCM. Technical Note, Bell Laboratories. IEEE Transactions on Information Theory, 28, 128–137.
  * [28] Löffler, M., Zhang, A. Y., and Zhou, H. H. (2021). Optimality of spectral clustering in the Gaussian mixture model. Annals of Statistics, 49, 2506-2530.
  * [29] MacQueen, J. B. (1967). Some Methods for classification and Analysis of Multivariate Observations. Proceedings of 5th Berkeley Symposium on Mathematical Statistics and Probability. University of California Press. pp. 281–297.
  * [30] Mahajan, M., Nimbhorkar, P., and Varadarajan, K. (2012). The planar k-means problem is NP-hard. Theoretical Computer Science, 442, 13-21.
  * [31] Maranzana, F.E. (1963). On the location of supply points to minimize transportation costs. IBM Systems Journal, 2, 129-135.
  * [32] Nityasuddhi, D., and Böhning, D. (2003). Asymptotic properties of the EM algorithm estimate for normal mixture models with component specific variances. Computational Statistics and Data Analysis, 41, 591-601.
  * [33] Pollard, D. (1982a). A central limit theorem for k-means clustering. Annals of Probability, 10, 919-926.
  * [34] Rasmussen, C. (1999). The infinite Gaussian mixture model. Advances in neural information processing systems, 12, 554–560
  * [35] Steinley, D. (2003). Local optima in K-means clustering: what you don’t know may hurt you. Psychological methods, 8, 294-304.
  * [36] Steinley, D. (2006). K-means clustering: A half-century synthesis. British Journal of Mathematical and Statistical Psychology, 59, 1-34.
  * [37] Steinhaus, H. (1957). Sur la division des corps matériels en parties. Bulletin L’Académie Polonaise des Science, 4, 801-804, (In French).
  * [38] Terada, Y. (2015). Strong consistency of factorial K-means clustering. Annals of the Institute of Statistical Mathematics, 67, 335-357.
  * [39] Tibshirani, R., Walther, G., and Hastie (2001). Estimating he number of clusters in data set via the gap statistic. Journal of Royal Statistical Society Series B, 63, 411-423.
  * [40] Trauwaert, E., Kaufman, L., and Rousseeuw, P. (1991). Fuzzy clustering algorithms based on the maximum likelihood principle. Fuzzy Sets and Systems, 42, 213-227.
  * [41] van der Vaart, A.W. (1998). Asymptotic Statistics, Cambridge University Press, Cambridge, UK.
  * [42] Witten, D.M., and Tibshirani, R. (2010). A framework for feature selection in clustering. Journal of the American Statistical Association, 105, 713-726.
  * [43] Wong, M. A. (1984). Asymptotic properties of univariate sample K-means clusters. Journal of Classification, 1, 255-270.
  * [44] Wu, C. J. (1983). On the convergence properties of the EM algorithm. Annals of statistics, 11, 95-103.
  * [45] Zhang, T., and Lin, G. (2021). Generalized k-means in GLMs with applications to the outbreak of COVID-19 in the United States. Computational Statistics and Data Analysis, 159, 107217.
  * [46] Zhao, Y. and Karypis, G. (2005). Hierarchical clustering algorithms for document datasets. Data Mining and Knowledge Discovery, 10, 141-168.