∎

11institutetext: Tran Hung Cuong 22institutetext: Department of Computer
Science, Faculty of Information Technology, Hanoi University of Industry, 298
Cau Dien Road, Bac Tu Liem District, Hanoi, Vietnam
22email<EMAIL_ADDRESS>33institutetext: Jen-Chih Yao
44institutetext: Center for General Education, China Medical University,
Taichung 40402, Taiwan
44email<EMAIL_ADDRESS>55institutetext: Nguyen Dong
Yen66institutetext: Institute of Mathematics Vietnam Academy of Science and
Technology, 18 Hoang Quoc Viet, Hanoi 10307, Vietnam
66email<EMAIL_ADDRESS>
# Qualitative Properties of the Minimum Sum-of-Squares Clustering Problem

Tran Hung Cuong  Jen-Chih Yao  Nguyen Dong Yen

(Received: date / Accepted: date)

###### Abstract

A series of basic qualitative properties of the minimum sum-of-squares
clustering problem are established in this paper. Among other things, we
clarify the solution existence, properties of the global solutions,
characteristic properties of the local solutions, locally Lipschitz property
of the optimal value function, locally upper Lipschitz property of the global
solution map, and the Aubin property of the local solution map. We prove that
the problem in question always has a global solution and, under a mild
condition, the global solution set is finite and the components of each global
solution can be computed by an explicit formula. Based on a newly introduced
concept of nontrivial local solution, we get necessary conditions for a system
of centroids to be a nontrivial local solution. Interestingly, we are able to
show that these necessary conditions are also sufficient ones. Finally, it is
proved that the optimal value function is locally Lipschitz, the global
solution map is locally upper Lipschitz, and the local solution map has the
Aubin property, provided that the original data points are pairwise distinct.
Thanks to the obtained complete characterizations of the nontrivial local
solutions, one can understand better the performance of not only the $k$-means
algorithm, but also of other solution methods for the minimum sum-of-squares
clustering problem.

###### Keywords:

Clustering problem solution existence local solution optimality condition
stability property

##  1 Introduction

Clustering is an important task in data mining and it is a powerful tool for
automated analysis of data. Cluster is a subset of the data set. The elements
of a cluster are similar in some sense (see, e.g., (Aggarwal 2014, , p. 32)
and (Kantardzic 2011, , p. 250)).

Clustering is an unsupervised technique dealing with problems of organizing a
collection of patterns into clusters based on similarity. Cluster analysis is
applied in different areas such as image segmentation, information retrieval,
pattern recognition, pattern classification, network analysis, vector
quantization and data compression, data mining and knowledge discovery
business, document clustering and image processing (see, e.g., (Aggarwal 2014,
, p. 32) and Kumar 2017 ).

There are many kinds of clustering problems, where different criteria are
used. Among these criteria, the Minimum Sum-of-Squares Clustering (MSSC for
short) criterion is one of the most used Bock_1998 ; Brusco 2006 ; Costa 2017
; Du Merle 2000 ; Kumar 2017 ; PhamDinh_LeThi_2009 ; Peng 2005 ; Shereli 2005
. Biding by this criterion, one tries to make the sum of the squared Euclidean
distances from each data point to the centroid of its cluster as small as
possible. The MSSC problem requires to partition a finite data set into a
given number of clusters in order to minimize the just mentioned sum.

The importance of the MSSC problem was noticed by researchers long time ago
and they have developed many algorithms to solve it (see, e.g., Bagirov_2008 ;
Bagirov_2014 ; Bagirov_Rubinov_2003 ; Bagirov_2011 ; Bagirov_Yearwood_2006 ;
Likas_2003 ; Ordin_Bagirov_2015 ; Xie_2011 , and the references therein).
Since this is a NP-hard problem Aloise 2009 ; Mahajan_NPHard_2009 , the
effective existing algorithms reach at most local solutions. These algorithms
may include certain techniques for improving the current data partition to
seek better solutions. For example, in Ordin_Bagirov_2015 , the authors
proposed a method to find good starting points that is based on the DCA
(Difference-of Convex-functions Algorithms). The latter has been applied to
the MSSC problem in Bagirov_2014 ; PhamDinh_LeThi_2007 ; PhamDinh_LeThi_2009 .

It is well known that a deep understanding on qualitative properties of an
optimization problem is very helpful for its numerical solution. To our
knowledge, apart from the fundamental necessary optimality condition given
recently by Ordin and Bagirov Ordin_Bagirov_2015 , qualitative properties of
the MSSC problem have not been addressed in the literature until now.

The first aim of the present paper is to prove some basic properties of the
above problem. We begin with clarifying the equivalence between the mixed
integer programming formulation and the unconstrained nonsmooth nonconvex
optimization formulation of the problem, that were given in Ordin_Bagirov_2015
. Then we prove that the MSSC problem always has a global solution and, under
a mild condition, the global solution set is finite and the components of each
global solution can be computed by an explicit formula.

The second aim of our paper is to characterize the local solutions of the MSSC
problem. Based on the necessary optimality condition in DC programming
Dem'yanov_Vasil'ev_1985 , some arguments of Ordin_Bagirov_2015 , and a newly
introduced concept of nontrivial local solution, we get necessary conditions
for a system of centroids to be a nontrivial local solution. Interestingly, we
are able to prove that these necessary conditions are also sufficient ones.
Since the known algorithms for solving the MSSC problem focus on the local
solutions, our characterizations may lead to a better understanding and
further refinements of the existing algorithms. Here, by constructing a
suitable example, we investigate the performance of the $k$-means algorithm,
which can be considered as a basic solution method for the MSSC problem.

The third aim of this paper is to analyze the changes of the optimal value,
the global solution set, and the local solution set of the MSSC problem with
respect to small changes in the data set. Three principal stability properties
will be established. Namely, we will prove that the optimal value function is
locally Lipschitz, the global solution map is locally upper Lipschitz, and the
local solution map has the Aubin property, provided that the original data
points are pairwise distinct.

The remainder of the paper consists of three sections. Section 2 describes the
MSSC problem and its basic properties. The $k$-means algorithm is recalled in
this section. In Section 3, we characterize the local solutions and
investigate the performance of the just mentioned algorithm. Section 4
presents a comprehensive stability analysis of the MSSC problem.

##  2 Basic Properties of the MSSC Problem

Let $A=\\{{\mathbf{a}}^{1},...,{\mathbf{a}}^{m}\\}$ be a finite set of points
(representing the data points to be grouped) in the $n$-dimensional Euclidean
space $\mathbb{R}^{n}$ equipped with the scalar product
$\langle\mathbf{x},\mathbf{y}\rangle=\displaystyle\sum_{i=1}^{n}x_{i}y_{i}$
and the norm
$\|\mathbf{x}\|=\left(\displaystyle\sum_{i=1}^{n}x_{i}^{2}\right)^{1/2}$.
Given a positive integer $k$ with $k\leq m$, one wants to partition $A$ into
disjoint subsets $A^{1},\dots,A^{k},$ called clusters, such that a clustering
criterion is optimized.

If one associates to each cluster $A^{j}$ a center (or centroid), denoted by
${\mathbf{x}}^{j}\in\mathbb{R}^{n}$, then the following well-known variance or
SSQ (Sum-of-Squares) clustering criterion (see, e.g., (Bock_1998, , p. 266))

$\psi(\mathbf{x},\alpha):=\frac{1}{m}\sum_{i=1}^{m}\left(\sum_{j=1}^{k}\alpha_{ij}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)\
\;\longrightarrow\ \;\min,$

where $\alpha_{ij}=1$ if ${\mathbf{a}}^{i}\in A^{j}$ and $\alpha_{ij}=0$
otherwise, is used. Thus, the above partitioning problem can be formulated as
the constrained optimization problem

$\begin{array}[]{rl}\min\Big{\\{}\psi(\mathbf{x},\alpha)\;\mid&\in\mathbb{R}^{n\times
k},\ \mathbb{\alpha}=(\alpha_{ij})\in\mathbb{R}^{m\times k},\
\alpha_{ij}\in\\{0,1\\},\\\ &\displaystyle\sum_{j=1}^{k}\alpha_{ij}=1,\
i=1,\dots,m,\ j=1,\dots,k\Big{\\}},\end{array}$ (1)

where the centroid system
$\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})$ and the incident matrix
$\mathbb{\alpha}=(\alpha_{ij})$ are to be found.

Since (1) is a difficult mixed integer programming problem, instead of it one
usually considers (see, e.g., (Ordin_Bagirov_2015, , p. 344)) next
unconstrained nonsmooth nonconvex optimization problem

$\min\Big{\\{}f(\mathbf{x}):=\frac{1}{m}\sum_{i=1}^{m}\left(\min_{j=1,\dots,k}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)\,\mid\,\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}\Big{\\}}.$ (2)

Both models (1) and (2) are referred to as the minimum sum-of-squares
clustering problem (the MSSC problem). As the decision variables of (1) and
(2) belong to different Euclidean spaces, the equivalence between these
minimization problems should be clarified. For our convenience, put
$I=\\{1,\dots,m\\}$ and $J=\\{1,\dots,k\\}$.

Given a vector
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$, we inductively construct $k$ subsets $A^{1},\dots,A^{k}$ of $A$ in the
following way. Put $A^{0}=\emptyset$ and

$A^{j}=\left\\{{\mathbf{a}}^{i}\in
A\setminus\left(\bigcup_{p=0}^{j-1}A^{p}\right)\;\mid\;\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|=\displaystyle\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\right\\}$ (3)

for $j\in J$. This means that, for every $i\in I$, the data point
${\mathbf{a}}^{i}$ belongs to the cluster $A^{j}$ if and only if the distance
$\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|$ is the minimal one among the
distances $\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|$, $q\in J$, and
${\mathbf{a}}^{i}$ does not belong to any cluster $A^{p}$ with $1\leq p\leq
j-1$. We will call this family $\\{A^{1},\dots,A^{k}\\}$ the natural
clustering associated with $\bar{\mathbf{x}}$.

###### Definition 1

Let
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$. We say that the component $\bar{\mathbf{x}}^{j}$ of $\bar{\mathbf{x}}$ is
attractive with respect to the data set $A$ if the set

$A[\bar{\mathbf{x}}^{j}]:=\left\\{{\mathbf{a}}^{i}\in
A\;\mid\;\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|=\displaystyle\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\right\\}$

is nonempty. The latter is called the attraction set of
$\bar{\mathbf{x}}^{j}$.

Clearly, the cluster $A^{j}$ in (3) can be represented as follows:

$A^{j}=A[\bar{\mathbf{x}}^{j}]\setminus\left(\bigcup_{p=1}^{j-1}A^{p}\right).$

###### Proposition 1

If $(\bar{\mathbf{x}},\bar{\alpha})$ is a solution of (1), then
$\bar{\mathbf{x}}$ is a solution of (2). Conversely, if $\bar{\mathbf{x}}$ is
a solution of (2), then the natural clustering defined by (3) yields an
incident matrix $\bar{\alpha}$ such that $(\bar{\mathbf{x}},\bar{\alpha})$ is
a solution of (1).

###### Proof

First, suppose that $(\bar{\mathbf{x}},\bar{\alpha})$ is a solution of (1). As
$\psi(\bar{\mathbf{x}},\bar{\alpha})\leq\psi(\bar{\mathbf{x}},\alpha)$ for
every $\alpha=(\alpha_{ij})\in\mathbb{R}^{m\times k}$ with
$\alpha_{ij}\in\\{0,1\\}$, $\sum_{j=1}^{k}\alpha_{ij}=1$ for all $i\in I$ and
$j\in J$, one must have

$\sum_{j=1}^{k}\bar{\alpha}_{ij}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|^{2}=\min_{j\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|^{2}\quad(\forall i\in I).$

Hence, $\psi(\bar{\mathbf{x}},\bar{\alpha})=f(\bar{\mathbf{x}})$. If
$\bar{\mathbf{x}}$ is not a solution of (2), then one can find some
$\tilde{\mathbf{x}}=(\tilde{\mathbf{x}}^{1},\dots,\tilde{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ such that $f(\tilde{\mathbf{x}})<f(\bar{\mathbf{x}})$. Let
$\\{A^{1},\dots,A^{k}\\}$ be the natural clustering associated with
$\tilde{\mathbf{x}}$. For any $i\in I$ and $j\in J$, set
$\tilde{\alpha}_{ij}=1$ if ${\mathbf{a}}^{i}\in A^{j}$ and
$\tilde{\alpha}_{ij}=0$ if ${\mathbf{a}}^{i}\notin A^{j}$. Let
$\tilde{\alpha}=(\tilde{\alpha}_{ij})\in\mathbb{R}^{m\times k}$. From the
definition of natural clustering and the choice of $\tilde{\alpha}$ it follows
that $\psi(\tilde{\mathbf{x}},\tilde{\alpha})=f(\tilde{\mathbf{x}})$. Then, we
have

$\psi(\bar{\mathbf{x}},\bar{\alpha})=f(\bar{\mathbf{x}})>f(\tilde{\mathbf{x}})=\psi(\tilde{\mathbf{x}},\tilde{\alpha}),$

contrary to the fact that $(\bar{\mathbf{x}},\bar{\alpha})$ is a solution of
(1).

Now, suppose that $\bar{\mathbf{x}}$ is a solution of (2). Let
$\\{A^{1},\dots,A^{k}\\}$ be the natural clustering associated with
$\bar{\mathbf{x}}$. Put $\bar{\alpha}=(\bar{\alpha}_{ij})$, where
$\bar{\alpha}_{ij}=1$ if ${\mathbf{a}}^{i}\in A^{j}$ and $\bar{\alpha}_{ij}=0$
if ${\mathbf{a}}^{i}\notin A^{j}$. It is easy to see that
$\psi(\bar{\mathbf{x}},\bar{\alpha})=f(\bar{\mathbf{x}})$. If there is a
feasible point $(\mathbf{x},\alpha)$ of (1) such that
$\psi(\mathbf{x},\alpha)<\psi(\bar{\mathbf{x}},\bar{\alpha})$ then, by
considering the natural clustering $\\{\tilde{A}^{1},\dots,\tilde{A}^{k}\\}$
be associated with $\mathbf{x}$ and letting
$\tilde{\alpha}=(\tilde{\alpha}_{ij})$ with $\tilde{\alpha}_{ij}=1$ if
${\mathbf{a}}^{i}\in\tilde{A}^{j}$ and $\tilde{\alpha}_{ij}=0$ if
${\mathbf{a}}^{i}\notin\tilde{A}^{j}$, we have
$f(\mathbf{x})=\psi(\mathbf{x},\tilde{\alpha})\leq\psi(\mathbf{x},\alpha)$.
Then we get

$f(\mathbf{x})\leq\psi(\mathbf{x},\alpha)<\psi(\bar{\mathbf{x}},\bar{\alpha})=f(\bar{\mathbf{x}}),$

contrary to the global optimality of $\bar{\mathbf{x}}$ for (2). One has thus
proved that $(\bar{\mathbf{x}},\bar{\alpha})$ is a solution of (1).
$\hfill\Box$

###### Proposition 2

If ${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are pairwise distinct points and
$\\{A^{1},\dots,A^{k}\\}$ is the natural clustering associated with a global
solution $\bar{\mathbf{x}}$ of (2), then $A^{j}$ is nonempty for every $j\in
J$.

###### Proof

Indeed, if there is some $j_{0}\in J$ with $A^{j_{0}}=\emptyset$, then the
assumption $k\leq m$ implies the existence of an index $j_{1}\in J$ such that
$A^{j_{1}}$ contains at least two points. Since the elements of $A^{j_{1}}$
are pairwise distinct, one could find ${\mathbf{a}}^{i_{1}}\in A^{j_{1}}$ with
${\mathbf{a}}^{i_{1}}\neq\bar{\mathbf{x}}^{j_{1}}$. Setting
$\tilde{\mathbf{x}}^{j}=\bar{\mathbf{x}}^{j}$ for $j\in J\setminus\\{j_{0}\\}$
and $\tilde{\mathbf{x}}^{j_{0}}={\mathbf{a}}^{i_{1}}$, one can easily show
that

$f(\tilde{\mathbf{x}})-f(\bar{\mathbf{x}})\leq-\frac{1}{m}\|{\mathbf{a}}^{i_{1}}-\bar{\mathbf{x}}^{j_{1}}\|^{2}<0.$

This contradicts the assumption saying that $\bar{\mathbf{x}}$ is a global
solution of (2). $\hfill\Box$

###### Theorem 2.1

Both problems (1), (2) have solutions. If
${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are pairwise distinct points, then the
solution sets are finite. Moreover, in that case, if
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ is a global solution of (2), then the attraction set
$A[\bar{\mathbf{x}}^{j}]$ is nonempty for every $j\in J$ and one has

$\bar{\mathbf{x}}^{j}=\frac{1}{|I(j)|}\displaystyle\sum_{i\in
I(j)}{\mathbf{a}}^{i},$ (4)

where $I(j):=\left\\{i\in I\mid{\mathbf{a}}^{i}\in
A[\bar{\mathbf{x}}^{j}]\right\\}$ with $|\Omega|$ denoting the number of
elements of $\Omega$.

###### Proof

a) Solution existence: By the second assertion of Proposition 1, it suffices
to show that (2) has a solution. Since the minimum of finitely many continuous
functions is a continuous function, the objective function of (2) is
continuous on $\mathbb{R}^{n\times k}$. If $k=1$, then the formula for $f$ can
be rewritten as
$f({\mathbf{x}}^{1})=\displaystyle\frac{1}{m}\displaystyle\sum_{i=1}^{m}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{1}\|^{2}$.
This smooth, strongly convex function attains its unique global minimum on
$\mathbb{R}^{n}$ at the point $\bar{\mathbf{x}}^{1}={\mathbf{a}}^{0}$, where

${\mathbf{a}}^{0}:=\frac{1}{m}\displaystyle\sum_{i\in I}{\mathbf{a}}^{i}$ (5)

is the barycenter of the data set $A$ (see, e.g., (LeeTamYen_book, , pp.
24–25) for more details). To prove the solution existence of (2) for any
$k\geq 2$, put $\rho=\displaystyle\max_{i\in
I}\|{\mathbf{a}}^{i}-{\mathbf{a}}^{0}\|$, where ${\mathbf{a}}^{0}$ is defined
by (5). Denote by $\bar{B}({\mathbf{a}}^{0},2\rho)$ the closed ball in
$\mathbb{R}^{n}$ centered at ${\mathbf{a}}^{0}$ with radius $2\rho$, and
consider the optimization problem

$\min\Big{\\{}f(\mathbf{x})\,\mid\,\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k},\;{\mathbf{x}}^{j}\in\bar{B}({\mathbf{a}}^{0},2\rho),\;\forall j\in
J\Big{\\}}.$ (6)

By the Weierstrass theorem, (6) has a solution
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})$ with
$\bar{\mathbf{x}}^{j}$ satisfying the inequality
$\|\bar{\mathbf{x}}^{j}-{\mathbf{a}}^{0}\|\leq 2\rho$ for all $j\in J$. Take
an arbitrary point
$\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ and notice by the choice of $\bar{\mathbf{x}}$ that
$f(\bar{\mathbf{x}})\leq f(\mathbf{x})$ if
$\|{\mathbf{x}}^{j}-{\mathbf{a}}^{0}\|\leq 2\rho$ for all $j\in J$. If there
exists at least one index $j\in J$ with
$\|{\mathbf{x}}^{j}-{\mathbf{a}}^{0}\|>2\rho$, then denote the set of such
indexes by $J_{1}$ and define a vector
$\tilde{\mathbf{x}}=(\tilde{\mathbf{x}}^{1},\dots,\tilde{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ by putting $\tilde{\mathbf{x}}^{j}={\mathbf{x}}^{j}$ for every $j\in
J\setminus J_{1}$, and $\tilde{\mathbf{x}}^{j}={\mathbf{a}}^{0}$ for all $j\in
J_{1}$. For any $i\in I$, it is clear that
$\|{\mathbf{a}}^{i}-\tilde{\mathbf{x}}^{j}\|=\|{\mathbf{a}}^{i}-{\mathbf{a}}^{0}\|\leq\rho<\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|$
for every $j\in J_{1}$, and
$\|{\mathbf{a}}^{i}-\tilde{\mathbf{x}}^{j}\|=\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|$
for every $j\in J\setminus J_{1}$. So, we have $f(\tilde{\mathbf{x}})\leq
f(\mathbf{x})$. As $f(\bar{\mathbf{x}})\leq f(\tilde{\mathbf{x}})$, this
yields $f(\bar{\mathbf{x}})\leq f(\mathbf{x})$. We have thus proved that
$\bar{\mathbf{x}}$ is a solution of (2).

b) Finiteness of the solution sets and formulae for the solution components:
Suppose that ${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are pairwise distinct
points,
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ is a global solution of (2), and $\\{A^{1},\dots,A^{k}\\}$ is the natural
clustering associated with $\bar{\mathbf{x}}$. By Proposition 2,
$A^{j}\neq\emptyset$ for all $j\in J$. Since

$A^{j}\subset\left\\{{\mathbf{a}}^{i}\in
A\mid\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|=\displaystyle\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\right\\}$

and $A^{j}\neq\emptyset$ for every $j\in J$, we see that $|I(j)|\geq 1$ for
every $j\in J$. This implies that right-hand-side of (4) is well defined for
each $j\in J$. To justify that formula, we can argue as follows. Fix any $j\in
J$. Since

$\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|>\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\quad\forall i\notin I(j),$

there exists $\varepsilon>0$ such that

$\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|>\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\quad\forall i\notin I(j)$ (7)

for any ${\mathbf{x}}^{j}\in\bar{B}(\bar{\mathbf{x}}^{j},\varepsilon)$. For
each ${\mathbf{x}}^{j}\in\bar{B}(\bar{\mathbf{x}}^{j},\varepsilon)$, put
$\tilde{\mathbf{x}}=(\tilde{\mathbf{x}}^{1},\dots,\tilde{\mathbf{x}}^{k})$
with $\tilde{\mathbf{x}}^{q}:=\bar{\mathbf{x}}^{q}$ for every $q\in
J\setminus\\{j\\}$ and $\tilde{\mathbf{x}}^{j}:={\mathbf{x}}^{j}$. From the
inequality $f(\bar{\mathbf{x}})\leq f(\tilde{\mathbf{x}})$ and the validity of
(7) we can deduce that

$\begin{array}[]{rl}f(\bar{\mathbf{x}})&=\displaystyle\frac{1}{m}\displaystyle\sum_{i=1}^{m}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|^{2}\right)\\\
&=\displaystyle\frac{1}{m}\left[\displaystyle\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|^{2}+\displaystyle\sum_{i\in
I\setminus I(j)}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|^{2}\right)\right]\\\ &\leq
f(\tilde{\mathbf{x}})\\\
&=\displaystyle\frac{1}{m}\left[\displaystyle\sum_{i\in I(j)}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\tilde{\mathbf{x}}^{q}\|^{2}\right)+\displaystyle\sum_{i\in
I\setminus I(j)}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\tilde{\mathbf{x}}^{q}\|^{2}\right)\right]\\\
&=\displaystyle\frac{1}{m}\left[\displaystyle\sum_{i\in I(j)}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\tilde{\mathbf{x}}^{q}\|^{2}\right)+\displaystyle\sum_{i\in
I\setminus I(j)}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|^{2}\right)\right]\\\
&\leq\displaystyle\frac{1}{m}\left[\displaystyle\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}+\displaystyle\sum_{i\in
I\setminus I(j)}\left(\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|^{2}\right)\right].\end{array}$ (8)

Consider the function
$\varphi({\mathbf{x}}^{j}):=\displaystyle\frac{1}{m}\displaystyle\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}$,
${\mathbf{x}}^{j}\in\mathbb{R}^{n}$. Comparing the expression on the second
line of (8) with the one on the sixth line yields
$\varphi(\bar{\mathbf{x}}^{j})\leq\varphi({\mathbf{x}}^{j})$ for every
${\mathbf{x}}^{j}\in\bar{B}(\bar{\mathbf{x}}^{j},\varepsilon)$. Hence
$\varphi$ attains its local minimum at $\bar{\mathbf{x}}^{j}$. By the Fermat
Rule we have $\nabla\varphi(\bar{\mathbf{x}}^{j})=0$, which gives
$\displaystyle\sum_{i\in I(j)}({\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j})=0$. This
equality implies (4). Since there are only finitely many nonempty subsets
$\Omega\subset I$, the set ${\mathcal{B}}$ of vectors $\mathbf{b}_{\Omega}$
defined by formula
$\mathbf{b}_{\Omega}=\displaystyle\frac{1}{|\Omega|}\displaystyle\sum_{i\in\Omega}\mathbf{a}^{i}$
is finite. (Note that $\mathbf{b}_{\Omega}$ is the barycenter of the subsystem
$\\{\mathbf{a}^{i}\in A\mid i\in\Omega\\}$ of $A$.) According to (4), each
component of a global solution
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})$ of (2)
must belongs to ${\mathcal{B}}$, we can assert that the solution set of (2) is
finite, provided that ${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are pairwise
distinct points. By Proposition 1, if $(\bar{\mathbf{x}},\bar{\alpha})$ is a
solution of (1), then $\bar{\mathbf{x}}$ is a solution of (2). Since
$\bar{\alpha}=(\bar{\alpha}_{ij})\in\mathbb{R}^{m\times k}$ must satisfy the
conditions $\bar{\alpha}_{ij}\in\\{0,1\\}$ and
$\displaystyle\sum_{j=1}^{k}\bar{\alpha}_{ij}=1$ for all $i\in I$, $j\in J$,
it follows that the solution set of (1) is also finite. $\hfill\Box$

###### Proposition 3

If
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ is a global solution of (2), then the components of $\bar{\mathbf{x}}$ are
pairwise distinct, i.e.,
$\bar{\mathbf{x}}^{j_{1}}\neq\bar{\mathbf{x}}^{j_{2}}$ whenever $j_{2}\neq
j_{1}$.

###### Proof

On the contrary, suppose that there are distinct indexes $j_{1},j_{2}\in J$
satisfying $\bar{\mathbf{x}}^{j_{1}}=\bar{\mathbf{x}}^{j_{2}}$. As $k\leq n$,
one has $k-1<n$. So, there must exist $j_{0}\in J$ such that
$|A[\bar{\mathbf{x}}^{j_{0}}]|\geq 2$. Therefore, one can find a data point
${\mathbf{a}}^{i_{0}}\in A[\bar{\mathbf{x}}^{j_{0}}]$ with
${\mathbf{a}}^{i_{0}}\neq\bar{\mathbf{x}}^{j_{0}}$. Setting
$\tilde{\mathbf{x}}=(\tilde{\mathbf{x}}^{1},\dots,\tilde{\mathbf{x}}^{k})$
with $\tilde{\mathbf{x}}^{j}=\bar{\mathbf{x}}^{j}$ for every $j\in
J\setminus\\{j_{2}\\}$ and $\tilde{\mathbf{x}}^{j_{2}}={\mathbf{a}}^{i_{0}}$.
The construction of $\tilde{\mathbf{x}}$ yields

$f(\tilde{\mathbf{x}})-f(\bar{\mathbf{x}})\leq-\frac{1}{m}\|{\mathbf{a}}^{i_{0}}-\bar{\mathbf{x}}^{j_{0}}\|^{2}<0,$

which is impossible because $\bar{\mathbf{x}}$ is a global solution of (2).
$\hfill\Box$

###### Remark 1

If the points ${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are not pairwise
distinct, then the conclusions of Theorem 2.1 do not hold in general. Indeed,
let $A=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{2}\\}\subset\mathbb{R}^{2}$ with
${\mathbf{a}}^{1}={\mathbf{a}}^{2}$. For $k:=2$, let
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\bar{\mathbf{x}}^{2})$ with
$\bar{\mathbf{x}}^{1}={\mathbf{a}}^{1}$ and
$\bar{\mathbf{x}}^{2}\in\mathbb{R}^{2}$ being chosen arbitrarily. Since
$f(\bar{\mathbf{x}})=0$, we can conclude that $\bar{\mathbf{x}}$ is a global
solution of (2). So, the problem has an unbounded solution set. Similarly, for
a data set
$A=\\{{\mathbf{a}}^{1},\dots,{\mathbf{a}}^{4}\\}\subset\mathbb{R}^{2}$ with
${\mathbf{a}}^{1}={\mathbf{a}}^{2}$, ${\mathbf{a}}^{3}={\mathbf{a}}^{4}$, and
${\mathbf{a}}^{2}\neq{\mathbf{a}}^{3}$. For $k:=3$, let
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\bar{\mathbf{x}}^{2},\bar{\mathbf{x}}^{3})$
with $\bar{\mathbf{x}}^{1}={\mathbf{a}}^{1}$,
$\bar{\mathbf{x}}^{2}={\mathbf{a}}^{3}$, and
$\bar{\mathbf{x}}^{3}\in\mathbb{R}^{2}$ being chosen arbitrarily. By the
equality $f(\bar{\mathbf{x}})=0$ we can assert that $\bar{\mathbf{x}}$ is a
global solution of (2). This shows that the solution set of (2) is unbounded.
Notice also that, if
$\bar{\mathbf{x}}^{3}\notin\\{\bar{\mathbf{x}}^{1},\bar{\mathbf{x}}^{2}\\}$,
then formula (4) cannot be applied to $\bar{\mathbf{x}}^{3}$, because the
index set $I(3)=\left\\{i\in I\mid{\mathbf{a}}^{i}\in
A[\bar{\mathbf{x}}^{3}]\right\\}=\left\\{i\in
I\mid\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{3}\|=\displaystyle\min_{q\in
J}\|{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\right\\}$ is empty.

Formula (4) is effective for computing certain components of any given local
solution of (2). The precise statement of this result is as follows.

###### Theorem 2.2

If
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ is a local solution of (2), then (4) is valid for all $j\in J$ whose index
set $I(j)$ is nonempty, i.e., the component $\bar{\mathbf{x}}^{j}$ of
$\bar{\mathbf{x}}$ is attractive w.r.t. the data set $A$.

###### Proof

It suffices to re-apply the arguments described in the second part of the
proof of Theorem 2.1, noting that $f(\bar{\mathbf{x}})\leq
f(\tilde{\mathbf{x}})$ if ${\mathbf{x}}^{j}$ (the $j$-th component of
$\tilde{\mathbf{x}}$) is taken from
$\bar{B}(\bar{\mathbf{x}}^{j},\varepsilon^{\prime})$ with
$\varepsilon^{\prime}\in(0,\varepsilon)$ being small enough. $\hfill\Box$

As in the proof of Theorem 2.1, if
$\Omega=\\{{\mathbf{a}}^{i_{1}},\dots,{\mathbf{a}}^{i_{r}}\\}\subset A$ is a
nonempty subset, then we put
$\mathbf{b}_{\Omega}=\displaystyle\frac{1}{r}\sum_{l=1}^{r}{\mathbf{a}}^{i_{l}}$.
Recall that the set of such points $\mathbf{b}_{\Omega}$ has been denoted by
${\mathcal{B}}$.

###### Remark 2

Theorem 2.1 shows that if the points ${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$
are pairwise distinct, then every component of a global solution must belong
to ${\mathcal{B}}$. It is clear that ${\mathcal{B}}\subset{\rm co}A$, where
${\rm co}A$ abbreviates the convex hull of $A$. Looking back to the proof of
Theorem 2.1, we see that the set $A$ lies in the ball
$\bar{B}({\mathbf{a}}^{0},\rho)$. Hence ${\mathcal{B}}\subset{\rm
co}A\subset\bar{B}({\mathbf{a}}^{0},\rho)$. It follows that the global
solutions of (2) are contained in the set

$\Big{\\{}\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}\,\mid\,{\mathbf{x}}^{j}\in\bar{B}({\mathbf{a}}^{0},\rho),\;\forall j\in
J\Big{\\}},$

provided the points ${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are pairwise
distinct. Similarly, Theorem 2.2 assures that each attractive component of a
local solution of (2) belongs to ${\mathcal{B}}$, where
${\mathcal{B}}\subset{\rm co}A\subset\bar{B}({\mathbf{a}}^{0},\rho)$.

###### Remark 3

If
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},\dots,\bar{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ is a global solution (resp., a local solution) of (2) then, for any
permutation $\sigma$ of $J$, the vector
$\bar{\mathbf{x}}^{\sigma}:=(\bar{\mathbf{x}}^{\sigma(1)},\dots,\bar{\mathbf{x}}^{\sigma(k)})$
is also a global solution (resp., a local solution) of (2). This observation
follows easily from the fact that $f(\mathbf{x})=f({\mathbf{x}}^{\sigma})$,
where
$\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ and
${\mathbf{x}}^{\sigma}:=({\mathbf{x}}^{\sigma(1)},\dots,{\mathbf{x}}^{\sigma(k)})$.

To understand the importance of the above results and those to be established
in next two sections, let us recall first the $k$-means clustering algorithm.

The $k$-means Algorithm: Despite its ineffectiveness, the $k$-means clustering
algorithm (see, e.g., (Aggarwal 2014, , pp. 89–90), Jain 2010 , (Kantardzic
2011, , pp. 263–266), and MacQueen 1967 ) is one of the most popular solution
methods for (2). One starts with selecting $k$ points
${\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k}$ in $\mathbb{R}^{n}$ as the initial
centroids. Then one inductively constructs $k$ subsets $A^{1},\dots,A^{k}$ of
the data set $A$ by putting $A^{0}=\emptyset$ and using the rule (3), where
${\mathbf{x}}^{j}$ plays the role of $\bar{\mathbf{x}}^{j}$ for all $j\in J$.
This means that $\\{A^{1},\dots,A^{k}\\}$ is the natural clustering associated
with $\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})$. Once the clusters
are formed, for each $j\in J$, if $A^{j}\neq\emptyset$ then the centroid
${\mathbf{x}}^{j}$ is updated by the rule

${\mathbf{x}}^{j}\ \;\leftarrow\
\;\widetilde{\mathbf{x}}^{j}:=\frac{1}{|I(A^{j})|}\,\displaystyle\sum_{i\in
I(A^{j})}{\mathbf{a}}^{i}$ (9)

with $I(A^{j}):=\left\\{i\in I\mid{\mathbf{a}}^{i}\in A^{j}\right\\}$; and
${\mathbf{x}}^{j}$ does not change otherwise. The algorithm iteratively
repeats the procedure until the centroid system
$\\{{\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k}\\}$ is stable, i.e.,
$\widetilde{\mathbf{x}}^{j}={\mathbf{x}}^{j}$ for all $j\in J$ with
$A^{j}\neq\emptyset$. The computation procedure is described as follows.

Input: The data set $A=\\{{\mathbf{a}}^{1},...,{\mathbf{a}}^{m}\\}$ and a
constant $\varepsilon\geq 0$ (tolerance).
Output: The set of $k$ centroids
$\\{{\mathbf{x}}^{1},...,{\mathbf{x}}^{k}\\}$.
Step 1. Select initial centroids ${\mathbf{x}}^{j}\in\mathbb{R}^{n}$ for all
$j\in J$.
Step 2. Compute $\alpha_{i}=\min\\{\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|\mid
j\in J\\}$ for all $i\in I$.
Step 3. Form the clusters $A^{1},\dots,A^{k}$:
\- Find the attraction sets

$A[{\mathbf{x}}^{j}]=\left\\{{\mathbf{a}}^{i}\in
A\;\mid\;\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|=\alpha_{i}\right\\}\quad\;(j\in
J);$

\- Set $A^{1}=A[{\mathbf{x}}^{1}]$ and

$A^{j}=A[{\mathbf{x}}^{j}]\setminus\left(\bigcup_{p=1}^{j-1}A^{p}\right)\quad\;(j=2,\dots,k).$
(10)

Step 4. Update the centroids ${\mathbf{x}}^{j}$ satisfying
$A^{j}\neq\emptyset$ by the rule (9), keeping other centroids unchanged.
Step 5. Check the convergence condition: If
$\|\widetilde{\mathbf{x}}^{j}-{\mathbf{x}}^{j}\|\leq\varepsilon$ for all $j\in
J$ with $A^{j}\neq\emptyset$ then stop, else go to Step 2.

The following example is designed to show how the algorithm is performed in
practice.

###### Example 1

Choose $m=3$, $n=2$, and $k=2$. Let
$A=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{2},{\mathbf{a}}^{3}\\}$, where
${\mathbf{a}}^{1}=(0,0),\,{\mathbf{a}}^{2}=(1,0)$, ${\mathbf{a}}^{3}=(0,1)$.
Apply the $k$-means algorithm to solve the problem (2) with the tolerance
$\varepsilon=0$.

(a) With the starting centroids ${\mathbf{x}}^{1}={\mathbf{a}}^{1}$,
${\mathbf{x}}^{2}={\mathbf{a}}^{2}$, one obtains the clusters
$A^{1}=A[{\mathbf{x}}^{1}]=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{3}\\}$ and
$A^{2}=A[{\mathbf{x}}^{2}]=\\{{\mathbf{a}}^{2}\\}$. The updated centroids are
${\mathbf{x}}^{1}=(0,\frac{1}{2})$, ${\mathbf{x}}^{2}={\mathbf{a}}^{2}$. Then,
the new clusters $A^{1}$ and $A^{2}$ coincide with the old ones. Thus,
$\|\widetilde{\mathbf{x}}^{j}-{\mathbf{x}}^{j}\|=0$ for all $j\in J$ with
$A^{j}\neq\emptyset$. So, the computation terminates. For
${\mathbf{x}}^{1}=(0,\frac{1}{2})$, ${\mathbf{x}}^{2}={\mathbf{a}}^{2}$, one
has $f(\mathbf{x})=\frac{1}{6}$.

(b) Starting with the points ${\mathbf{x}}^{1}=(\frac{1}{4},\frac{3}{4})$ and
${\mathbf{x}}^{2}=(2,3)$, one gets the clusters
$A^{1}=A[{\mathbf{x}}^{1}]=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{2},{\mathbf{a}}^{3}\\}$
and $A^{2}=A[{\mathbf{x}}^{2}]=\emptyset$. The algorithm gives the centroid
system ${\mathbf{x}}^{1}=(\frac{1}{3},\frac{1}{3})$, ${\mathbf{x}}^{2}=(2,3)$,
and $f(\mathbf{x})=\frac{1}{3}$.

(c) Starting with ${\mathbf{x}}^{1}=(0,1)$ and ${\mathbf{x}}^{2}=(0,0)$, by
the algorithm we are led to
$A^{1}=A[{\mathbf{x}}^{1}]=\\{{\mathbf{a}}^{3}\\}$,
$A^{2}=A[{\mathbf{x}}^{2}]=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{2}\\}$,
${\mathbf{x}}^{1}=(0,1)$, and ${\mathbf{x}}^{2}=(\frac{1}{2},0)$. The
corresponding value of objective function is $f(\mathbf{x})=\frac{1}{6}$.

(d) Starting with ${\mathbf{x}}^{1}=(0,0)$ and
${\mathbf{x}}^{2}=(\frac{1}{2},\frac{1}{2})$, by the algorithm one gets the
results $A^{1}=A[{\mathbf{x}}^{1}]=\\{{\mathbf{a}}^{1}\\}$,
$A^{2}=A[{\mathbf{x}}^{2}]=\\{{\mathbf{a}}^{2},{\mathbf{a}}^{3}\\}$,
${\mathbf{x}}^{1}=(0,0)$, and ${\mathbf{x}}^{2}=(\frac{1}{2},\frac{1}{2})$.
The corresponding value of objective function is $f(\mathbf{x})=\frac{4}{9}$.

(e) With ${\mathbf{x}}^{1}=(\frac{1}{3},\frac{1}{3})$ and
${\mathbf{x}}^{2}=(1+\frac{\sqrt{5}}{3},0)$ as the initial centroids, one
obtains the results
$A^{1}=A[{\mathbf{x}}^{1}]=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{2},{\mathbf{a}}^{3}\\}$,
$A^{2}=A[{\mathbf{x}}^{2}]=\emptyset$,
${\mathbf{x}}^{1}=(\frac{1}{3},\frac{1}{3})$,
${\mathbf{x}}^{2}=(1+\frac{\sqrt{5}}{3},0)$, and $f(\mathbf{x})=\frac{4}{9}$.

Based on the existing knowledge on the MSSC problem and the $k$-means
clustering algorithm, one cannot know whether the five centroid systems
obtained in the items (a)–(e) of Example 1 contain a global optimal solution
of the clustering problem, or not. Even if one knows that the centroid systems
obtained in (a) and (c) are global optimal solutions, one still cannot say
definitely whether the centroid systems obtained in the items (b), (d), (e)
are local optimal solutions of (2), or not.

The theoretical results in this and the forthcoming sections allow us to
clarify the following issues related to the MSSC problem in Example 1:

\- The structure of the global solution set (see Example 2 below);

\- The structure of the local solution set (see Example 3);

\- The performance of the $k$-means algorithm (see Example 4).

In particular, it will be shown that the centroid systems in (a) and (c) are
global optimal solutions, the centroid systems in (b) and (d) are local-
nonglobal optimal solutions, while the centroid system in (e) is not a local
solution (despite the fact that the centroid systems generated by the
$k$-means algorithm converge to it, and the value of the objective function at
it equals to the value given by the centroid system in (d)).

##  3 Characterizations of the Local Solutions

In order to study the local solution set of (2) in more details, we will
follow Ordin and Bagirov Ordin_Bagirov_2015  to consider the problem in light
of a well-known optimality condition in DC programming. For every
$\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times k}$,
we have

$\begin{array}[]{rl}f(\mathbf{x})&=\displaystyle\frac{1}{m}\sum_{i\in
I}\left(\min_{j\in J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)\\\
&=\displaystyle\frac{1}{m}\sum_{i\in I}\left[\left(\sum_{j\in
J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)-\left(\max_{j\in
J}\sum_{q\in
J\setminus\\{j\\}}\|{\mathbf{a}}^{q}-{\mathbf{x}}^{j}\|^{2}\right)\right].\end{array}$
(11)

Hence, the objective function $f$ of (2) can be expressed (Ordin_Bagirov_2015,
, p. 345) as the difference of two convex functions

$\displaystyle f(\mathbf{x})=f^{1}(\mathbf{x})-f^{2}(\mathbf{x}),$ (12)

where

$\displaystyle f^{1}(\mathbf{x}):=\displaystyle\frac{1}{m}\sum_{i\in
I}\left(\sum_{j\in J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)$ (13)

and

$\displaystyle f^{2}(\mathbf{x}):=\frac{1}{m}\sum_{i\in I}\left(\max_{j\in
J}\sum_{q\in
J\setminus\\{j\\}}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|^{2}\right).$ (14)

It is clear that $f^{1}$ is a convex linear-quadratic function. In particular,
it is differentiable. As the sum of finitely many nonsmooth convex functions,
$f^{2}$ is a nonsmooth convex function, which is defined on the whole space
$\mathbb{R}^{n\times k}$. The subdifferentials of $f^{1}(\mathbf{x})$ and
$f^{2}(\mathbf{x})$ can be computed as follows. First, one has

$\begin{array}[]{rl}\partial f^{1}(\mathbf{x})=\\{\nabla
f^{1}(\mathbf{x})\\}&=\left\\{\displaystyle\frac{2}{m}\sum_{i\in
I}\left({\mathbf{x}}^{1}-{\mathbf{a}}^{i},\dots,{\mathbf{x}}^{k}-{\mathbf{a}}^{i}\right)\right\\}\\\
&=\\{2({\mathbf{x}}^{1}-{\mathbf{a}}^{0},\dots,{\mathbf{x}}^{k}-{\mathbf{a}}^{0})\\}\end{array}$

where, as before, ${\mathbf{a}}^{0}=\mathbf{b}_{A}$ is the barycenter of the
system $\\{{\mathbf{a}}^{1},\dots,{\mathbf{a}}^{m}\\}$. Set

$\displaystyle\varphi_{i}(\mathbf{x})=\max_{j\in J}\,h_{i,j}(\mathbf{x})$ (15)

with $h_{i,j}(\mathbf{x}):=\displaystyle\sum_{q\in
J\setminus\\{j\\}}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|^{2}$ and

$\displaystyle J_{i}(\mathbf{x})=\left\\{j\in
J\,\mid\,h_{i,j}(\mathbf{x})=\varphi_{i}(\mathbf{x})\right\\}.$ (16)

###### Proposition 4

One has

$\displaystyle J_{i}(\mathbf{x})=\left\\{j\in J\,\mid\,{\mathbf{a}}^{i}\in
A[{\mathbf{x}}^{j}]\right\\}.$ (17)

###### Proof

From the formula of $h_{i,j}(\mathbf{x})$ it follows that

$\displaystyle h_{i,j}(\mathbf{x})=\left(\displaystyle\sum_{q\in
J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|^{2}\right)-\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}.$

Therefore, by (15), we have

$\displaystyle\begin{array}[]{rl}\varphi_{i}(\mathbf{x})&=\displaystyle\max_{j\in
J}\,\left[\left(\displaystyle\sum_{q\in
J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|^{2}\right)-\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right]\\\
&=\left(\displaystyle\sum_{q\in
J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|^{2}\right)+\displaystyle\max_{j\in
J}\,\left(-\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)\\\
&=\left(\displaystyle\sum_{q\in
J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|^{2}\right)-\displaystyle\min_{j\in
J}\,\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}.\end{array}$

Thus, the maximum in (15) is attained when the minimum
$\displaystyle\min_{j\in J}\,\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}$ is
achieved. So, by (16),

$\displaystyle J_{i}(\mathbf{x})=\left\\{j\in
J\,\mid\,\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|=\displaystyle\min_{q\in
J}\,\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|\right\\}.$

This implies (17). $\hfill\Box$

Invoking the subdifferential formula for the maximum function (see
(Clarke_1990, , Proposition 2.3.12) and note that the Clarke generalized
gradient coincides with the subdifferential of convex analysis if the
functions in question are convex), we have

$\displaystyle\partial\varphi_{i}(\mathbf{x})={\rm co}\left\\{\nabla
h_{i,j}(\mathbf{x})\,|\,j\in J_{i}(\mathbf{x})\right\\}={\rm
co}\left\\{2\left(\widetilde{\mathbf{x}}^{j}-\widetilde{\mathbf{a}}^{i,j}\right)\,\mid\,j\in
J_{i}(\mathbf{x})\right\\},$ (19)

where

$\displaystyle\widetilde{\mathbf{x}}^{j}=\left({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{j-1},\mathbf{0}_{\mathbb{R}^{n}},{\mathbf{x}}^{j+1},\dots,{\mathbf{x}}^{k}\right)$
(20)

and

$\displaystyle\widetilde{\mathbf{a}}^{i,j}=\left({\mathbf{a}}^{i},\dots,{\mathbf{a}}^{i},\underbrace{\mathbf{0}_{\mathbb{R}^{n}}}_{j-{\rm
th\ position}},{\mathbf{a}}^{i},\dots,{\mathbf{a}}^{i}\right).$ (21)

By the Moreau-Rockafellar theorem (Rockafellar_1970, , Theorem 23.8), one has

$\displaystyle\partial f^{2}(\mathbf{x})=\frac{1}{m}\sum_{i\in
I}\partial\varphi_{i}(\mathbf{x})$ (22)

with $\partial\varphi_{i}(\mathbf{x})$ being computed by (19).

Now, suppose
$\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times k}$
is a local solution of (2). By the necessary optimality condition in DC
programming (see, e.g., NTVHang_Yen_JOTA2016  and PhamDinh_LeThi_AMV97 ),
which can be considered as a consequence of the optimality condition obtained
by Dem’yanov et al. in quasidifferential calculus (see, e.g.,
(Dem'yanov_Rubinov_1995, , Theorem 3.1) and (Dem'yanov_Vasil'ev_1985, ,
Theorem 5.1)), we have

$\partial f^{2}(\mathbf{x})\subset\partial f^{1}(\mathbf{x}).$ (23)

Since $\partial f^{1}(\mathbf{x})$ is a singleton, $\partial
f^{2}(\mathbf{x})$ must be a singleton too. This happens if and only if
$\partial\varphi_{i}(\mathbf{x})$ is a singleton for every $i\in I$. By (19),
if $|J_{i}(\mathbf{x})|=1$, then $|\partial\varphi_{i}(\mathbf{x})|=1$. In the
case where $|J_{i}(\mathbf{x})|>1$, we can select two elements $j_{1}$ and
$j_{2}$ from $J_{i}(\mathbf{x})$, $j_{1}<j_{2}$. As
$\partial\varphi_{i}(\mathbf{x})$ is a singleton, by (19) one must have
$\widetilde{\mathbf{x}}^{j_{1}}-\widetilde{\mathbf{a}}^{i,j_{1}}=\widetilde{\mathbf{x}}^{j_{2}}-\widetilde{\mathbf{a}}^{i,j_{2}}$.
Using (20) and (21), one sees that the latter occurs if and only if
${\mathbf{x}}^{j_{1}}={\mathbf{x}}^{j_{2}}={\mathbf{a}}^{i}$. To proceed
furthermore, we need to introduce the following condition on the local
solution $\mathbf{x}$.

(C1) The components of $\mathbf{x}$ are pairwise distinct, i.e.,
${\mathbf{x}}^{j_{1}}\neq{\mathbf{x}}^{j_{2}}$ whenever $j_{2}\neq j_{1}$.

###### Definition 2

A local solution $\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})$ of (2)
that satisfies (C1) is called a nontrivial local solution.

###### Remark 4

Proposition 3 shows that every global solution of (2) is a nontrivial local
solution.

The following fundamental facts have the origin in (Ordin_Bagirov_2015, , pp.
346). Here, a more precise and complete formulation is presented. In
accordance with (17), the first assertion of next theorem means that if
$\mathbf{x}$ is a nontrivial local solution, then for each data point
${\mathbf{a}}^{i}\in A$ there is a unique component ${\mathbf{x}}^{j}$ of
$\mathbf{x}$ such that ${\mathbf{a}}^{i}\in A[{\mathbf{x}}^{j}]$.

###### Theorem 3.1

(Necessary conditions for nontrivial local optimality) Suppose that
$\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})$ is a nontrivial local
solution of (2). Then, for any $i\in I$, $|J_{i}(\mathbf{x})|=1$. Moreover,
for every $j\in J$ such that the attraction set $A[{\mathbf{x}}^{j}]$ of
${\mathbf{x}}^{j}$ is nonempty, one has

${\mathbf{x}}^{j}=\frac{1}{|I(j)|}\displaystyle\sum_{i\in
I(j)}{\mathbf{a}}^{i},$ (24)

where $I(j)=\left\\{i\in I\mid{\mathbf{a}}^{i}\in
A[{\mathbf{x}}^{j}]\right\\}$. For any $j\in J$ with
$A[{\mathbf{x}}^{j}]=\emptyset$, one has

${\mathbf{x}}^{j}\notin{\mathcal{A}}[\mathbf{x}],$ (25)

where ${\mathcal{A}}[\mathbf{x}]$ is the union of the balls
$\bar{B}({\mathbf{a}}^{p},\|{\mathbf{a}}^{p}-{\mathbf{x}}^{q}\|)$ with $p\in
I$, $q\in J$ satisfying $p\in I(q)$.

###### Proof

Suppose $\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})$ is a nontrivial
local solution of (2). Given any $i\in I$, we must have
$|J_{i}(\mathbf{x})|=1$. Indeed, if $|J_{i}(\mathbf{x})|>1$ then, by the
analysis given before the formulation of the theorem, there exist indexes
$j_{1}$ and $j_{2}$ from $J_{i}(\mathbf{x})$ such that
${\mathbf{x}}^{j_{1}}={\mathbf{x}}^{j_{2}}={\mathbf{a}}^{i}$. This contradicts
the nontriviality of the local solution $\mathbf{x}$. Let
$J_{i}(\mathbf{x})=\\{j(i)\\}$ for $i\in I$, i.e., $j(i)\in J$ is the unique
element of $J_{i}(\mathbf{x})$.

For each $i\in I$, observe by (15) that

$h_{i,j}(\mathbf{x})<h_{i,j(i)}(\mathbf{x})=\varphi_{i}(\mathbf{x})\ \;\forall
j\in J\setminus\\{j(i)\\}.$

Hence, by the continuity of the functions $h_{i,j}(\mathbf{x})$, there exists
an open neighborhood $U_{i}$ of $\mathbf{x}$ such that

$h_{i,j}(\mathbf{y})<h_{i,j(i)}(\mathbf{y})\ \;\forall j\in
J\setminus\\{j(i)\\},\;\forall\mathbf{y}\in U_{i}.$

It follows that

$\varphi_{i}(\mathbf{y})=h_{i,j(i)}(\mathbf{y})\ \;\forall\mathbf{y}\in
U_{i}.$ (26)

So, $\varphi_{i}(.)$ is continuously differentiable on $U_{i}$. Put
$U=\displaystyle\bigcap_{i\in I}U_{i}$. From (14) and (26) one can deduce that

$\displaystyle f^{2}(\mathbf{y})=\frac{1}{m}\sum_{i\in
I}\varphi_{i}(\mathbf{y})=\frac{1}{m}\sum_{i\in
I}h_{i,j(i)}(\mathbf{y})\quad\forall\mathbf{y}\in U.$

Therefore, $f^{2}(\mathbf{y})$ is continuously differentiable function on $U$.
Moreover, the formulas (19)–(21) yield

$\displaystyle\nabla f^{2}(\mathbf{y})=\frac{2}{m}\sum_{i\in
I}(\widetilde{\mathbf{y}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)})\quad\forall\mathbf{y}\in
U,$ (27)

where

$\displaystyle\widetilde{\mathbf{y}}^{j(i)}=\big{(}\mathbf{y}^{1},...,\mathbf{y}^{j(i)-1},\mathbf{0}_{\mathbb{R}^{n}},\mathbf{y}^{j(i)+1},...,\mathbf{y}^{k}\big{)}$

and

$\displaystyle\widetilde{\mathbf{a}}^{i,j(i)}=\left({\mathbf{a}}^{i},\dots,{\mathbf{a}}^{i},\underbrace{\mathbf{0}_{\mathbb{R}^{n}}}_{j(i)-{\rm
th\ position}},{\mathbf{a}}^{i},\dots,{\mathbf{a}}^{i}\right).$

Substituting $\mathbf{y}=\mathbf{x}$ into (27) and combining the result with
(23), we obtain

$\displaystyle\sum_{i\in
I}(\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)})=m({\mathbf{x}}^{1}-{\mathbf{a}}^{0},...,{\mathbf{x}}^{k}-{\mathbf{a}}^{0}).$
(28)

Now, fix an index $j\in J$ with $A[{\mathbf{x}}^{j}]\neq\emptyset$ and
transform the left-hand side of (28) as follows:

$\displaystyle\sum_{i\in
I}(\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)})$
$\displaystyle=\sum_{i\in
I,\,j(i)=j}(\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)})+\sum_{i\in
I,\,j(i)\neq
j}(\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)})$
$\displaystyle=\sum_{i\in
I,\,j(i)=j}(\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)})+\sum_{i\notin
I(j)}(\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)}).$

Clearly, if $j(i)=j$, then the $j$-th component of the vector
$\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)}$, that belongs
to $\mathbb{R}^{n\times k}$, is $\mathbf{0}_{\mathbb{R}^{n}}$. If $j(i)\neq
j$, then the $j$-th component of the vector
$\widetilde{\mathbf{x}}^{j(i)}-\widetilde{\mathbf{a}}^{i,j(i)}$ is
${\mathbf{x}}^{j}-{\mathbf{a}}^{i}$. Consequently, (28) gives us

$\displaystyle\sum_{i\notin
I(j)}({\mathbf{x}}^{j}-{\mathbf{a}}^{i})=m({\mathbf{x}}^{j}-{\mathbf{a}}^{0}).$

Since $m{\mathbf{a}}^{0}={\mathbf{a}}^{1}+\dots+{\mathbf{a}}^{m}$, this yields
$\displaystyle\sum_{i\in I(j)}{\mathbf{a}}^{i}=|I(j)|{\mathbf{x}}^{j}.$ Thus,
formula (24) is valid for any $j\in J$ satisfying
$A[{\mathbf{x}}^{j}]\neq\emptyset$.

For any $j\in J$ with $A[{\mathbf{x}}^{j}]=\emptyset$, one has (25). Indeed,
suppose to the contrary that there exits $j_{0}\in J$ with
$A[{\mathbf{x}}^{j_{0}}]=\emptyset$ such that for some $p\in I$, $q\in J$, one
has $p\in I(q)$ and
${\mathbf{x}}^{j_{0}}\in\bar{B}({\mathbf{a}}^{p},\|{\mathbf{a}}^{p}-{\mathbf{x}}^{q}\|)$.
If
$\|{\mathbf{a}}^{p}-{\mathbf{x}}^{j_{0}}\|=\|{\mathbf{a}}^{p}-{\mathbf{x}}^{q}\|$,
then $J_{p}(\mathbf{x})\supset\\{q,j_{0}\\}$. This is impossible due to the
first claim of the theorem. Now, if
$\|{\mathbf{a}}^{p}-{\mathbf{x}}^{j_{0}}\|<\|{\mathbf{a}}^{p}-{\mathbf{x}}^{q}\|$,
then $p\notin I(q)$. We have thus arrived at a contradiction.

The proof is complete. $\hfill\Box$

Roughly speaking, the necessary optimality condition given in the above
theorem is a sufficient one. Therefore, in combination with Theorem 3.1, next
statement gives a complete description of the nontrivial local solutions of
(2).

###### Theorem 3.2

(Sufficient conditions for nontrivial local optimality) Suppose that a vector
$\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times k}$
satisfies condition (C1) and $|J_{i}(\mathbf{x})|=1$ for every $i\in I$. If
(24) is valid for any $j\in J$ with $A[{\mathbf{x}}^{j}]\neq\emptyset$ and
(25) is fulfilled for any $j\in J$ with $A[{\mathbf{x}}^{j}]=\emptyset$, then
$\mathbf{x}$ is a nontrivial local solution of (2).

###### Proof

Let $\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ be such that (C1) holds, $J_{i}(\mathbf{x})=\\{j(i)\\}$ for every $i\in
I$, (24) is valid for any $j\in J$ with $A[{\mathbf{x}}^{j}]\neq\emptyset$,
and (25) is satisfied for any $j\in J$ with $A[{\mathbf{x}}^{j}]=\emptyset$.
Then, for all $i\in I$ and $j^{\prime}\in J\setminus\\{j(i)\\}$, one has

$\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j(i)}\|<\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j^{\prime}}\|.$

So, there exist $\varepsilon>0$, $q\in J$, such that

$\displaystyle\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j(i)}\|<\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j^{\prime}}\|\quad\forall
i\in I,\;\forall j^{\prime}\in J\setminus\\{j(i)\\},$ (29)

whenever vector
$\widetilde{\mathbf{x}}=(\widetilde{\mathbf{x}}^{1},...,\widetilde{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ satisfies the condition
$\|\widetilde{\mathbf{x}}^{q}-{\mathbf{x}}^{q}\|<\varepsilon$ for all $q\in
J$. By (25) and by the compactness of ${\mathcal{A}}[\mathbf{x}]$, reducing
the positive number $\varepsilon$ (if necessary) we have

$\widetilde{\mathbf{x}}^{j}\notin{\mathcal{A}}[\widetilde{\mathbf{x}}]$ (30)

whenever vector
$\widetilde{\mathbf{x}}=(\widetilde{\mathbf{x}}^{1},...,\widetilde{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ satisfies the condition
$\|\widetilde{\mathbf{x}}^{q}-{\mathbf{x}}^{q}\|<\varepsilon$ for all $q\in
J$, where ${\mathcal{A}}[\widetilde{\mathbf{x}}]$ is the union of the balls
$\bar{B}({\mathbf{a}}^{p},\|{\mathbf{a}}^{p}-\widetilde{\mathbf{x}}^{q}\|)$
with $p\in I$, $q\in J$ satisfying $p\in I(q)=\left\\{i\in
I\mid{\mathbf{a}}^{i}\in A[{\mathbf{x}}^{q}]\right\\}$.

Fix an arbitrary vector
$\widetilde{\mathbf{x}}=(\widetilde{\mathbf{x}}^{1},...,\widetilde{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ with the property that
$\|\widetilde{\mathbf{x}}^{q}-{\mathbf{x}}^{q}\|<\varepsilon$ for all $q\in
J$. Then, by (29) and (30), $J_{i}(\widetilde{\mathbf{x}})=\\{j(i)\\}$. So,

$\min_{j\in
J}\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j}\|^{2}=\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j(i)}\|^{2}.$

Therefore, one has

$\begin{array}[]{rl}f(\widetilde{\mathbf{x}})&=\displaystyle\frac{1}{m}\sum_{i\in
I}\left(\min_{j\in
J}\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j}\|^{2}\right)\\\
&=\displaystyle\frac{1}{m}\sum_{i\in
I}\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j(i)}\|^{2}\\\
&=\displaystyle\frac{1}{m}\sum_{j\in J}\left(\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j(i)}\|^{2}\right)\\\
&=\displaystyle\frac{1}{m}\sum_{j\in J}\left(\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j}\|^{2}\right)\\\
&\geq\displaystyle\frac{1}{m}\sum_{j\in J}\left(\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right)=f(\mathbf{x}),\end{array}$

where the inequality is valid because (24) obviously yields

$\sum_{i\in I(j)}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\leq\sum_{i\in
I(j)}\|{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j}\|^{2}$

for every $j\in J$ such that the attraction set $A[{\mathbf{x}}^{j}]$ of
${\mathbf{x}}^{j}$ is nonempty. (Note that ${\mathbf{x}}^{j}$ is the
barycenter of $A[{\mathbf{x}}^{j}]$.)

The local optimality of $\mathbf{x}=({\mathbf{x}}^{1},...,{\mathbf{x}}^{k})$
has been proved. Hence, $\mathbf{x}$ is a nontrivial local solution of
(2).$\hfill\Box$

###### Example 2

(A local solution need not be a global solution) Consider the clustering
problem described in Example 1. Here, $I=\\{1,2,3\\}$ and $J=\\{1,2\\}$. By
Theorem 2.1, problem (2) has a global solution. Moreover, if
$\mathbf{x}=({\mathbf{x}}^{1},{\mathbf{x}}^{2})\in\mathbb{R}^{2\times 2}$ is a
global solution then, for every $j\in J$, the attraction set
$A[{\mathbf{x}}^{j}]$ is nonempty. Thanks to Remark 4, we know that
$\mathbf{x}$ is a nontrivial local solution. So, by Theorem 3.1, the
attraction sets $A[{\mathbf{x}}^{1}]$ and $A[{\mathbf{x}}^{2}]$ are disjoint.
Moreover, the barycenter of each one of these sets can be computed by formula
(24). Clearly, $A=A[{\mathbf{x}}^{1}]\cup A[{\mathbf{x}}^{2}]$. Since
$A[{\mathbf{x}}^{j}]\subset
A=\\{{\mathbf{a}}^{1},{\mathbf{a}}^{2},{\mathbf{a}}^{3}\\}$, allowing
permutations of the components of each vector
$\mathbf{x}=({\mathbf{x}}^{1},{\mathbf{x}}^{2})\in\mathbb{R}^{2\times 2}$ (see
Remark 3), we can assert that the global solution set of our problem is
contained in the set

$\left\\{\bar{\mathbf{x}}:=\big{(}(\frac{1}{2},\frac{1}{2}),\,(0,0)\big{)},\
\hat{\mathbf{x}}:=\big{(}(0,\frac{1}{2}),\,(1,0)\big{)},\
\widetilde{\mathbf{x}}:=\big{(}(\frac{1}{2},0),(0,1)\big{)}\right\\}.$ (31)

Since $f(\bar{\mathbf{x}})=\frac{1}{3}$ and
$f(\hat{\mathbf{x}})=f(\widetilde{\mathbf{x}})=\frac{1}{6}$, we infer that
$\hat{\mathbf{x}}$ and $\widetilde{\mathbf{x}}$ are global solutions of our
problem. Using Theorem 3.2, we can assert that $\bar{\mathbf{x}}$ is a local
solution. Thus, $\bar{\mathbf{x}}$ is a local solution which does not belong
to the global solution set, i.e., $\bar{\mathbf{x}}$ is a local-nonglobal
solution of our problem.

###### Example 3

(Complete description of the set of nontrivial local solutions) Again,
consider the MSSC problem given in Example 1. Allowing permutations of the
components of each vector in $\mathbb{R}^{2\times 2}$, by Theorems 3.1 and 3.2
we find that the set of nontrivial local solutions consists of the three
vectors described in (31) and all the vectors of the form
$\mathbf{x}=({\mathbf{x}}^{1},{\mathbf{x}}^{2})\in\mathbb{R}^{2\times 2}$,
where ${\mathbf{x}}^{1}=(\frac{1}{3},\frac{1}{3})$ and

${\mathbf{x}}^{2}\notin\bar{B}({\mathbf{a}}^{1},\|{\mathbf{a}}^{1}-{\mathbf{x}}^{1}\|)\cup\bar{B}({\mathbf{a}}^{2},\|{\mathbf{a}}^{2}-{\mathbf{x}}^{1}\|)\cup\bar{B}({\mathbf{a}}^{3},\|{\mathbf{a}}^{3}-{\mathbf{x}}^{1}\|).$

This set of nontrivial local solutions is unbounded and non-closed.

###### Example 4

(Convergence analysis of the $k$-means algorithm) Consider once again the
problem described in Example 1. By the results given in Example 3, the
centroid systems in items (a), (b), (c) and (d) of Example 1 are local
solutions. In addition, by Example 2, the centroid systems in the just
mentioned items (a) and (c) are global solutions. Concerning the centroid
system in item (e) of Example 1, remark that
$\mathbf{x}:=\big{(}(\frac{1}{3},\frac{1}{3}),(1+\frac{\sqrt{5}}{3},0)\big{)}$
is not a local solution by Theorem 3.1, because ${\mathbf{a}}^{2}\in
A[{\mathbf{x}}^{1}]\cap A[{\mathbf{x}}^{2}]$, i.e.,
$J_{2}(\mathbf{x})=\\{1,2\\}$ (see Figure 1). In general, with
${\mathbf{x}}^{1}=(\frac{1}{3},\frac{1}{3})$ and
${\mathbf{x}}^{2}\in\mathbb{R}^{2\times 2}$ belonging to the boundary of the
set

$\bar{B}({\mathbf{a}}^{1},\|{\mathbf{a}}^{1}-{\mathbf{x}}^{1}\|)\cup\bar{B}({\mathbf{a}}^{2},\|{\mathbf{a}}^{2}-{\mathbf{x}}^{1}\|)\cup\bar{B}({\mathbf{a}}^{3},\|{\mathbf{a}}^{3}-{\mathbf{x}}^{1}\|),$

$\mathbf{x}:=({\mathbf{x}}^{1},{\mathbf{x}}^{2})$ is not a local solution of
the MSSC problem under consideration. The above analysis shows that the
$k$-means algorithm is very sensitive to the choice of starting centroids. The
algorithm may give a global solution, a local-nonglobal solution, as well as a
centroid system which is not a local solution. In other words, the quality of
the obtained result greatly depends on the initial centroid system.

Figure 1: The centroids in item (e) of Example 1

##  4 Stability Properties

This section is devoted to establishing the local Lipschitz property of the
optimal value function, the local upper Lipschitz property of the global
solution map, and the local Lipschitz-like property of the local solution map
of (2).

Now, let the data set $A=\\{{\mathbf{a}}^{1},...,{\mathbf{a}}^{m}\\}$ of the
problem (2) be subject to change. Put
$\mathbf{a}=({\mathbf{a}}^{1},...,{\mathbf{a}}^{m})$ and observe that
$\mathbf{a}\in\mathbb{R}^{n\times m}$. Denoting by $v(\mathbf{a})$ the optimal
value of (2), one has

$v(\mathbf{a})=\min\\{f(\mathbf{x})\,\mid\,\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}\\}.$ (32)

The global solution set of (2), denoted by $F(\mathbf{a})$, is given by

$F(\mathbf{a})=\left\\{\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}\,\mid\,f(\mathbf{x})=v(\mathbf{a})\right\\}.$

Let us abbreviate the local solution set of (2) to $F_{1}(\mathbf{a})$. Note
that $F(\mathbf{a})\subset F_{1}(\mathbf{a})$, and the inclusion may be
strict.

###### Definition 3

A family $\\{I(j)\mid j\in J\\}$ of pairwise distinct, nonempty subsets of $I$
is said to be a partition of $I$ if $\displaystyle\bigcup_{j\in J}I(j)=I$.

From now on, let
$\bar{\mathbf{a}}=(\bar{\mathbf{a}}^{1},...,\bar{\mathbf{a}}^{m})\in\mathbb{R}^{n\times
m}$ be a fixed vector with the property that
$\bar{\mathbf{a}}^{1},...,\bar{\mathbf{a}}^{m}$ are pairwise distinct.

###### Theorem 4.1

(Local Lipschitz property of the optimal value function) The optimal value
function $v:\mathbb{R}^{n\times m}\to\mathbb{R}$ is locally Lipschitz at
$\bar{\mathbf{a}}$, i.e., there exist $L_{0}>0$ and $\delta_{0}>0$ such that

$|v(\mathbf{a})-v(\mathbf{a}^{\prime})|\leq
L_{0}\|\mathbf{a}-\mathbf{a}^{\prime}\|$

for all $\mathbf{a}$ and $\mathbf{a^{\prime}}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$ and
$\|\mathbf{a}^{\prime}-\bar{\mathbf{a}}\|<\delta_{0}$.

###### Proof

Denote by $\Omega$ the set of all the partitions of $I$. Every element
$\omega$ of $\Omega$ is a family $\\{I_{\omega}(j)\mid j\in J\\}$ of pairwise
distinct, nonempty subsets of $I$ with $\displaystyle\bigcup_{j\in
J}I_{\omega}(j)=I$. We associate to each pair $(\omega,\mathbf{a})$, where
$\mathbf{a}=({\mathbf{a}}^{1},...,{\mathbf{a}}^{m})\in\mathbb{R}^{n\times m}$
and $\omega\in\Omega$, a vector
$\mathbf{x}_{\omega}(\mathbf{a})=({\mathbf{x}}^{1}_{\omega}(\mathbf{a}),\dots,{\mathbf{x}}^{k}_{\omega}(\mathbf{a}))\in\mathbb{R}^{n\times
k}$ with

${\mathbf{x}}^{j}_{\omega}(\mathbf{a})=\frac{1}{|I_{\omega}(j)|}\displaystyle\sum_{i\in
I_{\omega}(j)}{\mathbf{a}}^{i}$ (33)

for every $j\in J$. By Theorem 2.1, problem (2) has solutions and the number
of the global solutions is finite, i.e., $F(\bar{\mathbf{a}})$ is nonempty and
finite. Moreover, for each
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},...,\bar{\mathbf{x}}^{k})\in
F(\bar{\mathbf{a}})$, one can find some $\omega\in\Omega$ satisfying
$\bar{\mathbf{x}}^{j}={\mathbf{x}}^{j}_{\omega}(\bar{\mathbf{a}})$ for all
$j\in J$. Let $\Omega_{1}=\left\\{\omega_{1},\dots,\omega_{r}\right\\}$ be the
set of the elements of $\Omega$ corresponding the global solutions. Then,

$f(\mathbf{x}_{\omega_{1}}(\bar{\mathbf{a}}),\bar{\mathbf{a}})<f(\mathbf{x}_{\omega}(\bar{\mathbf{a}}),\bar{\mathbf{a}})\quad\
(\forall\omega\in\Omega\setminus\Omega_{1}),$ (34)

where

$f(\mathbf{x},\mathbf{a})=\frac{1}{m}\sum_{i\in I}\left(\min_{j\in
J}\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}\right).$ (35)

For each pair $(i,j)\in I\times J$, the rule
$(\mathbf{x},\mathbf{a})\mapsto\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|^{2}$
defines a polynomial function on $\mathbb{R}^{n\times
k}\times\mathbb{R}^{n\times m}$. In particular, this function is locally
Lipschitz on its domain. So, by (Clarke_1990, , Prop. 2.3.6 and 2.3.12) we can
assert that the function $f(\mathbf{x},\mathbf{a})$ in (35) is locally
Lipschitz on $\mathbb{R}^{n\times k}\times\mathbb{R}^{n\times m}$.

Now, observe that for any $\omega\in\Omega$ and $j\in J$, the vector function
${\mathbf{x}}^{j}_{\omega}(.)$ in (33), which maps $\mathbb{R}^{n\times m}$ to
$\mathbb{R}^{n}$, is continuously differentiable. In particular, it is locally
Lipschitz on $\mathbb{R}^{n\times m}$.

For every $\omega\in\Omega$, from the above observations we can deduce that
the function
$g_{\omega}(\mathbf{a}):=f(\mathbf{x}_{\omega}(\mathbf{a}),\mathbf{a})$ is
locally Lipschitz on $\mathbb{R}^{n\times m}$. Rewriting (34) as

$g_{\omega_{1}}(\bar{\mathbf{a}})<g_{\omega}(\bar{\mathbf{a}})\quad\
(\forall\omega\in\Omega\setminus\Omega_{1})$

and using the continuity of the functions $g_{\omega}(.)$, we can find a
number $\delta_{0}>0$ such that

$g_{\omega_{1}}(\mathbf{a})<g_{\omega}(\mathbf{a})\quad\
(\forall\omega\in\Omega\setminus\Omega_{1})$ (36)

for all $\mathbf{a}$ satisfying $\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$.
Since $\bar{\mathbf{a}}^{1},...,\bar{\mathbf{a}}^{m}$ are pairwise distinct,
without loss of generality, we may assume that
${\mathbf{a}}^{1},...,{\mathbf{a}}^{m}$ are pairwise distinct for any
$\mathbf{a}=({\mathbf{a}}^{1},...,{\mathbf{a}}^{m})$ with
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$.

Now, consider a vector $\mathbf{a}=({\mathbf{a}}^{1},...,{\mathbf{a}}^{m})$
satisfying $\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$. By (36),
$f(\mathbf{x}_{\omega_{1}}(\mathbf{a}),\mathbf{a})<f(\mathbf{x}_{\omega}(\mathbf{a}),\mathbf{a})$
for all $\omega\in\Omega\setminus\Omega_{1}$. Since $f(.,\mathbf{a})$ is the
objective function of (2), this implies that the set
$\\{x_{\omega}(\mathbf{a})\mid\omega\in\Omega\setminus\Omega_{1}\\}$ does not
contain any global solution of the problem. Thanks to Theorem 2.1, we know
that the global solution set $F(\mathbf{a})$ of (2) is contained in the set
$\\{x_{\omega}(\mathbf{a})\mid\omega\in\Omega_{1}\\}$. Hence,

$F(\mathbf{a})\subset\\{x_{\omega}(\mathbf{a})\mid\omega\in\Omega_{1}\\}=\\{x_{\omega_{1}}(\mathbf{a}),\dots,x_{\omega_{r}}(\mathbf{a})\\}.$
(37)

Since $F(\mathbf{a})\neq\emptyset$, by (37) one has

$v(\mathbf{a})=\min\,\\{f(\mathbf{x},\mathbf{a})\mid\mathbf{x}\in
F(\mathbf{a})\\}=\min\,\\{f(\mathbf{x}_{\omega_{\ell}}(\mathbf{a}),\mathbf{a})\mid\ell=1,\dots,r\\}.$

Thus, we have proved that

$v(\mathbf{a})=\min\,\\{g_{\omega_{\ell}}(\mathbf{a})\mid\ell=1,\dots,r\\}$
(38)

for all $\mathbf{a}$ satisfying $\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$.
As it has been noted, the functions $g_{\omega}$, $\omega\in\Omega$, are
locally Lipschitz on $\mathbb{R}^{n\times m}$. Hence, applying (Clarke_1990, ,
Prop. 2.3.6 and 2.3.12) to the minimum function in (38), we can assert that
$v$ is locally Lipschitz at $\bar{\mathbf{a}}$.

The proof is complete. $\hfill\Box$

###### Theorem 4.2

(Local upper Lipschitz property of the global solution map) The global
solution map $F:\mathbb{R}^{n\times m}\rightrightarrows\mathbb{R}^{n\times k}$
is locally upper Lipschitz at $\bar{\mathbf{a}}$, i.e., there exist $L>0$ and
$\delta>0$ such that

$F({\mathbf{a}})\subset
F(\bar{\mathbf{a}})+L\|\mathbf{a}-\bar{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}$ (39)

for all $\mathbf{a}$ satisfying $\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta$. Here

$\bar{B}_{\mathbb{R}^{n\times
k}}:=\left\\{\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}\;\mid\;\displaystyle\sum_{j\in J}\|{\mathbf{x}}^{j}\|\leq 1\right\\}$

denotes the closed unit ball of the product space $\mathbb{R}^{n\times k}$,
which is equipped with the sum norm $\|\mathbf{x}\|=\displaystyle\sum_{j\in
J}\|{\mathbf{x}}^{j}\|$.

###### Proof

Let $\Omega$, $\Omega_{1}=\left\\{\omega_{1},\dots,\omega_{r}\right\\}$,
$\mathbf{x}_{\omega}(\mathbf{a})=({\mathbf{x}}^{1}_{\omega}(\mathbf{a}),\dots,{\mathbf{x}}^{k}_{\omega}(\mathbf{a}))\in\mathbb{R}^{n\times
k}$, and $\delta_{0}$ be constructed as in the proof of the above theorem. For
any $\omega\in\Omega$, the vector function $\mathbf{x}_{\omega}(.)$, which
maps $\mathbb{R}^{n\times m}$ to $\mathbb{R}^{n\times k}$, is continuously
differentiable. Hence, there exist $L_{\omega}>0$ and $\delta_{\omega}>0$ such
that

$\|\mathbf{x}_{\omega}(\mathbf{a})-x_{\omega}(\widetilde{\mathbf{a}})\|\leq
L_{\omega}\|\mathbf{a}-\widetilde{\mathbf{a}}\|$ (40)

for any $\mathbf{a},\,\widetilde{\mathbf{a}}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{\omega}$ and
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{\omega}$. Set

$L=\max\\{L_{\omega_{1}},\dots,L_{\omega_{r}}\\}\ \;{\rm and}\
\;\delta=\min\\{\delta_{0},\delta_{\omega_{1}}\dots,\delta_{\omega_{r}}\\}.$

Then, for every $\mathbf{a}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta$, by (37) and (40) one has

$\displaystyle\begin{array}[]{rcl}F(\mathbf{a})\subset\\{x_{\omega_{1}}(\mathbf{a}),\dots,x_{\omega_{r}}(\mathbf{a})\\}&\subset&\\{x_{\omega_{1}}(\bar{\mathbf{a}}),\dots,x_{\omega_{r}}(\bar{\mathbf{a}})\\}+L\|\mathbf{a}-\bar{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}\\\
&=&F(\bar{\mathbf{a}})+L\|\mathbf{a}-\bar{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}.\end{array}$

Hence, inclusion (39) is valid for every $\mathbf{a}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta$. $\hfill\Box$

###### Theorem 4.3

(Aubin property of the local solution map) Let
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},...,\bar{\mathbf{x}}^{k})$ be an
element of $F_{1}(\bar{\mathbf{a}})$ satisfying condition (C1), that is,
$\bar{\mathbf{x}}^{j_{1}}\neq\bar{\mathbf{x}}^{j_{2}}$ whenever $j_{2}\neq
j_{1}$. Then, the local solution map $F_{1}:\mathbb{R}^{n\times
m}\rightrightarrows\mathbb{R}^{n\times k}$ has the Aubin property at
$(\bar{\mathbf{a}},\bar{\mathbf{x}})$, i.e., there exist $L_{1}>0$,
$\varepsilon>0$, and $\delta_{1}>0$ such that

$F_{1}(\mathbf{a})\cap B(\bar{\mathbf{x}},\varepsilon)\subset
F_{1}(\widetilde{\mathbf{a}})+L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}$ (42)

for all $\mathbf{a}$ and $\widetilde{\mathbf{a}}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{1}$ and
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{1}$.

###### Proof

Suppose that
$\bar{\mathbf{x}}=(\bar{\mathbf{x}}^{1},...,\bar{\mathbf{x}}^{k})\in
F_{1}(\bar{\mathbf{a}})$ and
$\bar{\mathbf{x}}^{j_{1}}\neq\bar{\mathbf{x}}^{j_{2}}$ for all $j_{1},j_{2}\in
J$ with $j_{2}\neq j_{1}$. Denote by $J_{1}$ the set of the indexes $j\in J$
such that $\bar{\mathbf{x}}^{j}$ is attractive w.r.t. the data set
$\\{\bar{\mathbf{a}}^{1},\dots,\bar{\mathbf{a}}^{m}\\}$. Put $J_{2}=J\setminus
J_{1}$. For every $j\in J_{1}$, by Theorem 3.1 one has

$\|\bar{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{j}\|<\|\bar{\mathbf{a}}^{i}-\bar{\mathbf{x}}^{q}\|\quad\;(\forall
i\in I(j),\;\forall q\in J\setminus\\{j\\}).$ (43)

In addition, the following holds:

$\bar{\mathbf{x}}^{j}=\frac{1}{|I(j)|}\displaystyle\sum_{i\in
I(j)}\bar{\mathbf{a}}^{i},$ (44)

where $I(j)=\left\\{i\in I\mid\bar{\mathbf{a}}^{i}\in
A[\bar{\mathbf{x}}^{j}]\right\\}$. For every $j\in J_{2}$, by Theorem 3.1 one
has

$\|\bar{\mathbf{x}}^{q}-\bar{\mathbf{a}}^{p}\|<\|\bar{\mathbf{x}}^{j}-\bar{\mathbf{a}}^{p}\|\quad\;(\forall
q\in J_{1},\;\forall p\in I(q)).$ (45)

Let $\varepsilon_{0}>0$ be such that
$\|\bar{\mathbf{x}}^{j_{1}}-\bar{\mathbf{x}}^{j_{2}}\|>\varepsilon_{0}$ for
all $j_{1},j_{2}\in J$ with $j_{2}\neq j_{1}$.

By (43) and (45), there exist $\delta_{0}>0$ and
$\varepsilon\in\left(0,\displaystyle\frac{\varepsilon_{0}}{4}\right)$ such
that

$\|{\mathbf{a}}^{i}-{\mathbf{x}}^{j}\|<\|{\mathbf{a}}^{i}-{\mathbf{x}}^{q}\|\quad\;(\forall
j\in J_{1},\;\forall i\in I(j),\;\forall q\in J\setminus\\{j\\})$ (46)

and

$\|{\mathbf{x}}^{q}-{\mathbf{a}}^{p}\|<\|{\mathbf{x}}^{j}-{\mathbf{a}}^{p}\|\quad\;(\forall
j\in J_{2},\;\forall q\in J_{1},\;\forall p\in I(q))$ (47)

for all
$\mathbf{a}=({\mathbf{a}}^{1},...,{\mathbf{a}}^{m})\in\mathbb{R}^{n\times m}$
and
$\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ satisfying $\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$ and
$\|\mathbf{x}-\bar{\mathbf{x}}\|<2k\varepsilon$. As
$\bar{\mathbf{x}}^{j_{1}}\neq\bar{\mathbf{x}}^{j_{2}}$ for all $j_{1},j_{2}\in
J$ with $j_{2}\neq j_{1}$, by taking a smaller $\varepsilon>0$ (if necessary),
for any
$\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})\in\mathbb{R}^{n\times
k}$ satisfying $\|\mathbf{x}-\bar{\mathbf{x}}\|<2k\varepsilon$ we have
${\mathbf{x}}^{j_{1}}\neq{\mathbf{x}}^{j_{2}}$ for all $j_{1},j_{2}\in J$ with
$j_{2}\neq j_{1}$.

For every $j\in J_{1}$ and
$\mathbf{a}=({\mathbf{a}}^{1},...,{\mathbf{a}}^{m})\in\mathbb{R}^{n\times m}$,
define

${\mathbf{x}}^{j}(\mathbf{a})=\frac{1}{|I(j)|}\displaystyle\sum_{i\in
I(j)}{\mathbf{a}}^{i}.$ (48)

Comparing (48) with (44) yields
${\mathbf{x}}^{j}(\bar{\mathbf{a}})=\bar{\mathbf{x}}^{j}$ for all $j\in
J_{1}$. Then, by the continuity of the vector functions ${\mathbf{x}}^{j}(.)$,
where $j\in J_{1}$, we may assume that
$\|{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})-\bar{\mathbf{x}}^{j}\|<\varepsilon$
for all $j\in J_{1}$ and
$\widetilde{a}=(\widetilde{\mathbf{a}}^{1},...,\widetilde{\mathbf{a}}^{m})\in\mathbb{R}^{n\times
m}$ satisfying $\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{0}$ (one
can take a smaller $\delta_{0}>0$, if necessary).

Since the vector functions ${\mathbf{x}}^{j}(.)$, $j\in J_{1}$, are
continuously differentiable, there exist $L_{1}>0$ such that

$\|{\mathbf{x}}^{j}(\mathbf{a})-{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})\|\leq\frac{1}{k}L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|$
(49)

for any $\mathbf{a},\,\widetilde{\mathbf{a}}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{0}$ and
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{0}$ (one can take a
smaller $\delta_{0}>0$, if necessary). Choose $\delta_{1}\in(0,\delta_{0})$ as
small as $\frac{2}{k}L_{1}\delta_{1}<\varepsilon$.

With the chosen constants $L_{1}>0$, $\varepsilon>0$, and $\delta_{1}>0$, let
us show that the inclusion (42) is fulfilled for all $\mathbf{a}$ and
$\widetilde{\mathbf{a}}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{1}$ and
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{1}$.

Let $\mathbf{a}$ and $\widetilde{\mathbf{a}}$ be such that
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{1}$ and
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{1}$. Select an arbitrary
element $\mathbf{x}=({\mathbf{x}}^{1},\dots,{\mathbf{x}}^{k})$ of the set
$F_{1}(\mathbf{a})\cap B(\bar{\mathbf{x}},\varepsilon)$. Put
$\widetilde{\mathbf{x}}^{j}={\mathbf{x}}^{j}(\widetilde{\mathbf{a}})$ for all
$j\in J_{1}$, where ${\mathbf{x}}^{j}(\mathbf{a})$ is given by (48). For any
$j\in J_{2}$, set $\widetilde{\mathbf{x}}^{j}={\mathbf{x}}^{j}$.

Claim 1. The vector
$\widetilde{\mathbf{x}}=(\widetilde{\mathbf{x}}^{1},\dots,\widetilde{\mathbf{x}}^{k})$
belongs to $F_{1}(\widetilde{\mathbf{a}})$.

Indeed, the inequalities $\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{1}$ and
$\|\mathbf{x}-\bar{\mathbf{x}}\|<\varepsilon$ imply that both properties (46)
and (47) are available. From (46) it follows that, for every $j\in J_{1}$, the
attraction set $A[{\mathbf{x}}^{j}]$ is $\\{{\mathbf{a}}^{i}\mid i\in
I(j)\\}$. Since $I(j)\neq\emptyset$ for each $j\in J_{1}$ and $\mathbf{x}\in
F_{1}(\mathbf{a})$, by Theorem 3.1 we have

${\mathbf{x}}^{j}=\frac{1}{|I(j)|}\displaystyle\sum_{i\in
I(j)}{\mathbf{a}}^{i}.$ (50)

Comparing (50) with (48) yields
${\mathbf{x}}^{j}={\mathbf{x}}^{j}(\mathbf{a})$ for all $j\in J_{1}$. By (47)
we see that, for every $j\in J_{2}$, the attraction set $A[{\mathbf{x}}^{j}]$
is empty. Moreover, one has

${\mathbf{x}}^{j}\notin{\mathcal{A}}[\mathbf{x}]\quad\;(\forall j\in J_{2})$
(51)

where ${\mathcal{A}}[\mathbf{x}]$ is the union of the balls
$\bar{B}({\mathbf{a}}^{p},\|{\mathbf{a}}^{p}-{\mathbf{x}}^{q}\|)$ with $p\in
I$, $q\in J$ satisfying $p\in I(q)$.

For each $j\in J_{1}$, using (49) we have

$\begin{array}[]{rl}\|{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})-\bar{\mathbf{x}}^{j}\|&\leq\|{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})-{\mathbf{x}}^{j}(\mathbf{a})\|+\|{\mathbf{x}}^{j}(\mathbf{a})-\bar{\mathbf{x}}^{j}\|\\\
&\leq\frac{1}{k}L_{1}\|\widetilde{\mathbf{a}}-\mathbf{a}\|+\varepsilon\\\
&\leq\frac{1}{k}L_{1}\left(\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|+\|\bar{\mathbf{a}}-\mathbf{a}\|\right)+\varepsilon\\\
&\leq\frac{2}{k}L_{1}\delta_{1}+\varepsilon<2\varepsilon.\end{array}$

Besides, for each $j\in J_{2}$, we have
$\|{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})-\bar{\mathbf{x}}^{j}\|=\|{\mathbf{x}}^{j}-\bar{\mathbf{x}}^{j}\|<\varepsilon$.
Therefore,

$\|\widetilde{\mathbf{x}}-\bar{\mathbf{x}}\|=\sum_{j\in
J_{1}}\|{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})-\bar{\mathbf{x}}^{j}\|+\sum_{j\in
J_{2}}\|{\mathbf{x}}^{j}-\bar{\mathbf{x}}^{j}\|<2k\varepsilon.$

In combination with the inequality
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{1}$, this assures that
the properties (46) and (47), where $\widetilde{\mathbf{a}}$ and
$\mathbf{x}(\widetilde{\mathbf{a}})$ respectively play the roles of
$\mathbf{a}$ and $\mathbf{x}$, hold. In other words, one has

$\|\widetilde{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{j}\|<\|\widetilde{\mathbf{a}}^{i}-\widetilde{\mathbf{x}}^{q}\|\quad\;(\forall
j\in J_{1},\;\forall i\in I(j),\;\forall q\in J\setminus\\{j\\})$ (52)

and

$\|\widetilde{\mathbf{x}}^{q}-\widetilde{\mathbf{a}}^{p}\|<\|\widetilde{\mathbf{x}}^{j}-\widetilde{\mathbf{a}}^{p}\|\quad\;(\forall
j\in J_{2},\;\forall q\in J_{1},\;\forall p\in I(q)).$ (53)

So, similar to the above case of $\mathbf{x}$, for every $j\in J_{1}$, the
attraction set $A[\widetilde{\mathbf{x}}^{j}]$ is
$\\{\widetilde{\mathbf{a}}^{i}\mid i\in I(j)\\}$. Recall that
$I(j)\neq\emptyset$ for each $j\in J_{1}$ and $\widetilde{\mathbf{x}}^{j}$ was
given by

$\widetilde{\mathbf{x}}^{j}={\mathbf{x}}^{j}(\widetilde{\mathbf{a}})=\frac{1}{|I(j)|}\displaystyle\sum_{i\in
I(j)}\widetilde{\mathbf{a}}^{i}.$ (54)

In addition, for every $j\in J_{2}$, the attraction set
$A[\widetilde{\mathbf{x}}^{j}]$ is empty and one has

$\widetilde{\mathbf{x}}^{j}\notin{\mathcal{A}}[\widetilde{\mathbf{x}}]\quad\;(\forall
j\in J_{2}),$ (55)

where ${\mathcal{A}}[\widetilde{\mathbf{x}}]$ is the union of the balls
$\bar{B}(\widetilde{\mathbf{a}}^{p},\|\widetilde{\mathbf{a}}^{p}-\widetilde{\mathbf{x}}^{q}\|)$
with $p\in I$, $q\in J$ satisfying $p\in I(q)$. Besides, from (52) and (53) it
follows that $|J_{i}(\widetilde{\mathbf{x}})|=1$ for every $i\in I$. Since
$\|\widetilde{\mathbf{x}}-\bar{\mathbf{x}}\|<2k\varepsilon$, we have
$\widetilde{\mathbf{x}}^{j_{1}}\neq\widetilde{\mathbf{x}}^{j_{2}}$ for all
$j_{1},j_{2}\in J$ with $j_{2}\neq j_{1}$. Due to the last two properties and
(54), (55), by Theorem 3.2 we conclude that $\widetilde{\mathbf{x}}\in
F_{1}(\widetilde{\mathbf{a}})$.

Claim 2. One has
$\mathbf{x}\in\widetilde{\mathbf{x}}+L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}$.

Indeed, since ${\mathbf{x}}^{j}={\mathbf{x}}^{j}(\mathbf{a})$ for all $j\in
J_{1}$, $\widetilde{\mathbf{x}}^{j}={\mathbf{x}}^{j}$ for any $j\in J_{2}$, by
(54) and (49) we have

$\displaystyle\|\mathbf{x}-\widetilde{\mathbf{x}}\|$ $\displaystyle=\sum_{j\in
J_{1}}\|{\mathbf{x}}^{j}-\widetilde{\mathbf{x}}^{j}\|+\sum_{j\in
J_{2}}\|{\mathbf{x}}^{j}-\widetilde{\mathbf{x}}^{j}\|$
$\displaystyle=\sum_{j\in
J_{1}}\|{\mathbf{x}}^{j}(\mathbf{a})-{\mathbf{x}}^{j}(\widetilde{\mathbf{a}})\|$
$\displaystyle\leq
k.\frac{1}{k}L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|=L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|.$

It follows that
${\mathbf{x}}\in\widetilde{\mathbf{x}}+L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}$.

Combining Claim 2 with Claim 1, we have ${\mathbf{x}}\in
F_{1}(\widetilde{\mathbf{a}})+L_{1}\|\mathbf{a}-\widetilde{\mathbf{a}}\|\bar{B}_{\mathbb{R}^{n\times
k}}$. Thus, property (42) is valid for all $\mathbf{a}$ and
$\widetilde{\mathbf{a}}$ satisfying
$\|\mathbf{a}-\bar{\mathbf{a}}\|<\delta_{1}$ and
$\|\widetilde{\mathbf{a}}-\bar{\mathbf{a}}\|<\delta_{1}$. $\hfill\Box$

## References

  * (1) Aggarwal, C.C., Reddy, C.K.: Data Clustering Algorithms and Applications. Chapman $\&$ Hall/CRC Press, Boca Raton, Florida (2014)
  * (2) Aloise, D., Deshpande, A., Hansen, P., Popat P.:NP-hardness of Euclidean sum-of-squares clustering. Mach. Learn. 75, 245–248 (2009)
  * (3) Bagirov, A.M.: Modified global $k$-means algorithm for minimum sum-of-squares clustering problems. Pattern Recognit. 41, 3192–3199 (2008)
  * (4) Bagirov, A.M.: An incremental DC algorithm for the minimum sum-of-squares clustering. Iranian J. Oper. Res. 5, 1–14 (2014)
  * (5) Bagirov, A.M., Rubinov, A.M., Soukhoroukova, N.V., Yearwood, J.: Unsupervised and supervised data classification via nonsmooth and global optimization. TOP 11, 1–93 (2003)
  * (6) Bagirov, A.M., Ugon, J., Webb, D.: Fast modified global $k$-means algorithm for incremental cluster construction. Pattern Recognit. 44, 866–876 (2011)
  * (7) Bagirov, A.M., Yearwood, J.: A new nonsmooth optimization algorithm for minimum sum-of-squares clustering problems. European J. Oper. Res. 170, 578–596 (2006)
  * (8) Bock, H.H.: Clustering and neural networks. In: Advances in Data Science and Classification, pp. 265–277, Springer, Berlin (1998)
  * (9) Brusco, M.J.: A repetitive branch-and-bound procedure for minimum within-cluster sum of squares partitioning. Psychometrika 71, 347–363 (2006)
  * (10) Clarke, F.H.: Optimization and Nonsmooth Analysis. SIAM, Philadelphia (1990)
  * (11) Costa, L.R., Aloise, D., Mladenović, N.: Less is more: basic variable neighborhood search heuristic for balanced minimum sum-of-squares clustering. Inform. Sci. 415/416, 247–253 (2017)
  * (12) Dem’yanov, V.F., Rubinov, A.M.: Constructive Nonsmooth Analysis. Peter Lang Verlag, Frankfurt am Main (1995)
  * (13) Dem’yanov, V.F., Vasil’ev, L.V.: Nondifferentiable Optimization. Translated from the Russian by T. Sasagawa. Optimization Software Inc., New York (1985)
  * (14) Du Merle, O., Hansen, P., Jaumard, B., Mladenović, N.: An interior point algorithm for minimum sum of squares clustering. SIAM J. Sci. Comput. 21, 1485–1505 (2000)
  * (15) Hang, N.T.V., Yen, N.D.: On the problem of minimizing a difference of polyhedral convex functions under linear constraints. J. Optim. Theory Appl. 171, 617–642 (2016)
  * (16) Jain, A.K., Data clustering: 50 years beyond $k$-means, Pattern Recognit. Lett. 31, 651–666 (2010)
  * (17) Kantardzic, M.: Data Mining Concepts, Models, Methods, and Algorithms, Second Edition. John Wiley $\&$ Sons, Hoboken, New Jersey (2011)
  * (18) Kumar, K.M., Reddy, A.R.M: An efficient $k$-means clustering filtering algorithm using density based initial cluster centers. Inform. Sci. 418/419, 286–301 (2017)
  * (19) Le Thi, H.A., Belghiti, M.T., Pham Dinh, T.: A new efficient algorithm based on DC programming and DCA for clustering. J. Global Optim. 37, 593–608 (2007)
  * (20) Le Thi, H.A., Pham Dinh, T.: Minimum sum-of-squares clustering by DC programming and DCA. ICIC 2009, LNAI 5755, 327–340 (2009)
  * (21) Lee, G.M., Tam, N.N., Yen, N.D.: Quadratic Programming and Affine Variational Inequalities: A Qualitative Study. Springer Verlag, New York (2005)
  * (22) Likas, A., Vlassis, N., Verbeek, J.J.: The global $k$-means clustering algorithm. Pattern Recognit. 36, 451–461 (2003)
  * (23) MacQueen, J.: Some methods for classification and analysis of multivariate observations. In: Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability, pp. 281–297 (1967)
  * (24) Mahajan, M., Nimbhorkar, P., Varadarajan, K.: The planar $k$-means problem is NP-hard. Theoret. Comput. Sci. 442, 13–21 (2012)
  * (25) Ordin, B., Bagirov, A.M.: A heuristic algorithm for solving the minimum sum-of-squares clustering problems. J. Global Optim. 61, 341–361 (2015)
  * (26) Peng, J., Xiay, Y.: A cutting algorithm for the minimum sum-of-squared error clustering. In: Proceedings of the SIAM International Data Mining Conference (2005)
  * (27) Pham Dinh, T., Le Thi, H.A.: Convex analysis approach to d.c. programming: theory, algorithms and applications. Acta Math. Vietnam. 22, 289–355 (1997)
  * (28) Rockafellar, R.T.: Convex Analysis. Princeton University Press, Princeton (1970)
  * (29) Sherali, H.D., Desai, J.: A global optimization RLT-based approach for solving the hard clustering problem. J. Global Optim. 32, 281–306 (2005)
  * (30) Xie, J., Jiang, S., Xie, W., Gao, X.: An efficient global $k$-means clustering algorithm. J. Comput. 6, 271–279 (2011)