# Optimized K-Means Algorithm

Authors: Samir Brahim Belhaouari; Shahnawaz Ahmed; Samer Mansour
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The localization of the region of interest (ROI), which contains the face, is the first step
in any automatic recognition system, which is a special case of the face detection. However, face
localization from input image is a challenging task due to possible variations in location, scale, pose,
occlusion, illumination, facial expressions, and clutter background. In this paper we introduce a new
optimized k-means algorithm that finds the optimal centers for each cluster which corresponds to the
global minimum of the k-means cluster. This method was tested to locate the faces in the input image
based on image segmentation. It separates the input image into two classes: faces and nonfaces. To
evaluate the proposed algorithm, MIT-CBCL, BioID, and Caltech datasets are used. The results show
significant localization accuracy.


---

## Body

## 1. Introduction

The face detection problem is to determine the face existence in the input image and then return its location if it exists. But in the face localization problem it is given that the input image contains a face and the goal is to determine the location of this face. The automatic recognition system, which is a special case of the face detection, locates, in its earliest phase, the region of interest (ROI) that contains the face. Face localization in an input image is a challenging task due to possible variations in location, scale, pose, occlusion, illumination, facial expressions, and clutter background. Various methods were proposed to detect and/or localize faces in an input image; however, there are still needs to improve the performance of localization and detection methods. A survey on some face recognition and detection techniques can be found in [1]. A more recent survey, mainly on face detection, was written by Yang et al. [2].One of the ROI detection trends is the idea of segmenting the image pixels into a number of groups or regions based on similarity properties. It gained more attention in the recognition technologies, which rely on grouping the features vector to distinguish between the image regions, and then concentrate on a particular region, which contains the face. One of the earliest surveys on image segmentation techniques was done by Fu and Mui [3]. They classify the techniques into three classes: features shareholding (clustering), edge detection, and region extraction. A later survey by N. R. Pal and S. K. Pal [4] did only concentrate on fuzzy, nonfuzzy, and colour images techniques. Many efforts had been done in ROI detection, which may be divided into two types: region growing and region clustering. The difference between the two types is that the region clustering searches for the clusters without prior information while the region growing needs initial points called seeds to detect the clusters. The main problem in region growing approach is to find the suitable seeds since the clusters will grow from the neighbouring pixels of these seeds based on a specified deviation. For seeds selection, Wan and Higgins [5] defined a number of criteria to select insensitive initial points for the subclass of region growing. To reduce the region growing time, Chang and Li [6] proposed a fast region growing method by using parallel segmentation.As mentioned before, the region clustering approach searches for clusters without prior information. Pappas and Jayant [7] generalized the K-means algorithm to group the image pixels based on spatial constraints and their intensity variations. These constraints include the following: first, the region intensity to be close to the data and, second, to have a spatial continuity. Their generalization algorithm is adaptive and includes the previous spatial constraints. Like the K-means clustering algorithm, their algorithm is iterative. These spatial constraints are included using Gibbs random field model and they concluded that each region is characterized by a slowly varying intensity function. Because of unimportant features, the algorithm is not that much accurate. Later, to improve their method, they proposed using a caricature of the original image to keep only the most significant features and remove the unimportant features [8]. Besides the problem of time costly segmentation, there are three basic problems that occur during the clustering process which are center redundancy, dead centers, and trapped centers at local minima. Moving K-mean (MKM) proposed in [9] can overcome the three basic problems by minimizing the center redundancy and the dead center as well as reducing indirectly the effect of trapped centers at local minima. In spite of that, MKM is sensitive to noise, centers are not located in the middle or centroid of a group of data, and some members of the centers with the largest fitness are assigned as members of a center with the smallest fitness. To reduce the effects of these problems. Isa et al. [10] proposed three modified versions of the MKM clustering algorithm for image segmentation: fuzzy moving k-means, adaptive moving k-means, and adaptive fuzzy moving k-means. After that Sulaiman and Isa [11] introduced a new segmentation method based on a new clustering algorithm which is Adaptive Fuzzy k-means Clustering Algorithm (AFKM). On the other hand, the Fuzzy C-Means (FCM) algorithm was used in image segmentation but it still has some drawbacks that are the lack of enough robustness to noise and outliers, the crucial parameter α being selected generally based on experience, and the time of segmenting of an image being dependent on its size. To overcome the drawbacks of the FCM algorithm, Cai et al. [12] integrated both local spatial and gray information to propose fast and robust Fuzzy C-means algorithm called Fast Generalized Fuzzy C-Means Algorithm (FGFCM) for image segmentation. Moreover, many researchers have provided a definition for some data characteristics where it has a significant impact on the K-means clustering analysis including the scattering of the data, noise, high-dimensionality, the size of the data and outliers in the data, datasets, types of attributes, and scales of attributes [13]. However, there is a need for more investigation to understand how data distributions affect K-means clustering performance. In [14] Xiong et al. provided a formal and organized study of the effect of skewed data distributions on K-means clustering. Previous research found the impact of high-dimensionality on K-means clustering performance. The traditional Euclidean notion of proximity is not very useful on real-world high-dimensional datasets, such as document datasets and gene-expression datasets. To overcome this challenge, researchers turn to make use of dimensionality reduction techniques, such as multidimensional scaling [15]. Second, K-means has difficulty in detecting the “natural” clusters with nonspherical shapes [13]. Modified K-means clustering algorithm is a direction to solve this problem. Guha et al. [16] proposed the CURE method, which makes use of multiple representative points to obtain the “natural” clusters shape information. The problem of outliers and noise in the data can also reduce clustering algorithms performance [17], especially for prototype-based algorithms such as K-means. The direction to solve this kind of problem is by combining some outlier removal techniques before conducting K-means clustering. For example, a simple method [9] of detecting outliers is based on the distance measure. On the other hand, many modified K-means clustering algorithms that work well for smaller medium-size datasets are unable to deal with large datasets. Bradley et al. in [18] introduced a discussion of scaling K-means clustering to large datasets.Arthur and Vassilvitskii implemented a preliminary version, k-means++, in C++ to evaluate the performance of k-means. K-means++ [19] uses a careful seeding method to optimize the speed and the accuracy of k-means. Experiments were performed on four datasets and results showed that that this algorithm is O
(
log
⁡
k
)-competitive with the optimal clustering. Experiments also proved its better speed (70% faster) and accuracy (potential value obtained is better by factors of 10 to 1000) than k-means. Kanungo et al. [20] also proposed an O
(
n
3
ϵ
-
d
) algorithm for the k-means problem. It is (
9
+
ϵ
) competitive but quiet slow in practice. Xiong et al. investigates the measures that best reflects the performance of K-means clustering [14]. An organized study was performed to understand the impact of data distributions on the performance of K-means clustering. Research work also has improved the traditional K-means clustering so that it can handle datasets with large variation of cluster sizes. This formal study illustrates that the entropy sometimes provided misleading information on the clustering performance. Based on this, a coefficient of variation (CV) is proposed to validate the clustering outcome. Experimental results proved that, for datasets with great variation in “true” cluster sizes, K-means lessens (less than 1.0) the variation in resultant cluster sizes. For datasets with small variation in “true” cluster sizes K-means increases (greater than 0.3) the variation in resultant cluster sizes.Scalability of the k-means for large datasets is one of the major limitations. Chiang et al. [21] proposed a simple and easy to implement algorithm for reducing the computation time of k-means. This pattern reduction algorithm compresses and/or removes iteration patterns that are not likely to change their membership. The pattern is compressed by selecting one of the patterns to be removed and setting its value to the average of all patterns removed. If D
[
r
l
i
] is the pattern to be removed then we can write it as follows:
(1)
D
[
r
l
i
]
=
∑
j
=
1
|
R
l
i
|
‍
D
[
r
i
j
l
]
,
where
r
l
i
∈
R
l
i
.This algorithm can be also applied to other clustering algorithms like population-based clustering algorithms. The computation time of many clustering algorithms can be reduced using this technique, as it is especially effective for large and high-dimensional datasets. In [22] Bradley et al. proposed Scalable k-means that uses buffering and a two-stage compression scheme to compress or remove patterns in order to improve the performance of k-means. The algorithm is slightly faster than k-means but does not show the same result always. The factors that degrade the performance of scalable k-means are the compression processes and compression ratio. Farnstrom et al. also proposed a simple single pass k-means algorithm [23], which reduces the computation time of k-means. Ordonez and Omiecinski proposed relational k-means [24] that uses the block and incremental concept to improve the stability of scalable k-means. The computation time of k-means can also be reduced using Parallel bisecting k-means [25] and triangle inequality [26] methods.Although K-means is the most popular and simple clustering algorithm, the major difficulty with this process is that it cannot ensure the global optimum results because the initial cluster centers are selected randomly. Reference [27] presents a novel technique that selects the initial cluster centers by using Voronoi diagram. This method automates the selection of the initial cluster centers. To validate the performance experiments were performed on a range of artificial and real world datasets. The quality of the algorithm is assessed using the following error rate (ER) equation:
(2)
ER
=
Number
of
misclassified
objects
Total
number
of
objects
×
100
%
.The lower the error rate the better the clustering. The proposed algorithm is able to produce better clustering results than the traditional K-means. Experimental results proved60
% reduction in iterations for defining the centroids.All previous solutions and efforts to increase the performance of K-means algorithm still need more investigation since they are all looking for a local minimum. Searching for the global minimum will certainly improve the performance of the algorithm. The rest of the paper is organized as follows. In Section2, the proposed method is introduced that finds the optimal centers for each cluster which corresponds to the global minimum of the k-means cluster. In Section 3 the results are discussed and analyzed using two sets from MIT, BioID, and Caltech datasets. Finally in Section 4, conclusions are drawn.


## 2. K-Means Clustering

K-means, an unsupervised learning algorithm first proposed by MacQueen, 1967 [28], is a popular method of cluster analysis. It is a process of organizing the specified objects into uniform classes called clusters based on similarities among objects based on certain criteria. It solves the well-known clustering problem by considering certain attributes and performing an iterative alternating fitting process. This algorithm partitions a dataset X into k disjoint clusters such that each observation belongs to the cluster with the nearest mean. Let x
-
i be the centroid of cluster c
i and let d
(
x
j
,
x
-
i
) be the dissimilarity between x
-
i and object x
j
∈
c
i. Then the function minimized by the k-means is given by the following equation:
(3)
min
⁡
x
-
1
,
…
,
x
-
k
E
=
∑
i
=
1
k
‍
∑
x
j
∈
c
i
‍
d
(
x
j
,
x
-
i
)
.K-means clustering is utilized in a vast number of applications including machine learning, fault detection, pattern recognition, image processing, statistics, and artificial intelligent [11, 29, 30]. k-means algorithm is considered as one of the fastest clustering algorithms with a number of variants that are sensitive to the selection of initial points and are intended for solving many issues of k-means like the evaluation of the clusters number [31], the method of initialization of the clusters centroids [32], and the speed of the algorithm [33].


### 2.1. Modifications in K-Means Clustering

Global k-means algorithm [34] is a proposal to improve global search properties of k-means algorithm and its performance on very large datasets by computing clusters successively. The main weakness of k-means clustering, that is, its sensitivity to initial positions of the cluster centers, is overcome by global k-means clustering algorithm which consists of a deterministic global optimization technique. It does not select initial values randomly; instead an incremental approach is applied to optimally add one new cluster center at each stage. Global k-means also proposes a method to reduce the computational load while maintain the solution quality. Experiments were performed to compare the performance of k-means, global k-means, min k-means, and fast global k-means as shown in Figure 1. Numerical results show that the global k-means algorithm considerably outperforms other k-means algorithms.Figure 1
Comparison of k-means.Bagirov [35] proposed a new version of the global k-means algorithm for minimum sum-of-squares clustering problems. He also compared three different versions of the k-means algorithm to propose the modified version of the global k-means algorithm. The proposed algorithm computes clusters incrementally and k
-
1 cluster centers from the previous iteration are used to compute k-partition of a dataset. The starting point of the kth cluster center is computed by minimizing an auxiliary cluster function. Given a finite set of points A in the n-dimensional space the global k-means compute the centroid X
1 of the set A as the following equation:
(4)
X
1
=
1
m
∑
i
=
1
m
‍
a
i
,
a
i
∈
A
,
i
=
1
,
…
,
m
.Numerical experiments performed on 14 datasets demonstrate the efficiency of the modified global k-means algorithm in comparison to the multistart k-means (MS k-means) and the GKM algorithms when the number of clustersk
>
5. Modified global k-means however requires more computational time than the global k-means algorithm. Xie and Jiang [36] also proposed a novel version of the global k-means algorithm. The method of creating the next cluster center in the global k-means algorithm is modified and that showed a better execution time without affecting the performance of the global k-means algorithm. Another extension of the standard k-means clustering is the Global Kernel k-means [37] which optimizes the clustering error in feature space by employing kernel-means as a local search procedure. A kernel function is used in order to solve the M-clustering problem and near-optimal solutions are used to optimize the clustering error in the feature space. This incremental algorithm has high ability to identify nonlinearly separable clusters in input space and no dependence on cluster initialization. It can handle weighted data points making it suitable for graph partitioning applications. Two major modifications were performed in order to reduce the computational cost with no major effect on the solution quality.Video imaging and Image segmentation are important applications for k-means clustering. Hung et al. [38] modified the k-means algorithm for color image segmentation where a weight selection procedure in the W-k-means algorithm is proposed. The evaluation function F
(
I
) is used for comparison:
(5)
F
(
I
)
=
R
∑
i
=
1
R
‍
(
e
i
/
256
)
2
A
i
,
where I = segmented image, R = the number of regions in the segmented image, A
i = the area, or the number of pixels of the ith region, and e
i = the color error of region i.e
i is the sum of the Euclidean distance of the color vectors between the original image and the segmented image. Smaller values of F
(
I
) give better segmentation results. Results from color image segmentation illustrate that the proposed procedure produces better segmentation than the random initialization. Maliatski and Yadid-Pecht [39] also propose an adaptive k-means clustering, hardware-driven approach. The algorithm uses both pixel intensity and pixel distance in the clustering process. In [40] also a FPGA implementation of real-time k-means clustering is done for colored images. A filtering algorithm is used for hardware implementation. Suliman and Isa [11] presented a unique clustering algorithm called adaptive fuzzy-K-means clustering for image segmentation. It can be used for different types of images including medical images, microscopic images, and CCD camera images. The concepts of fuzziness and belongingness are applied to produce more adaptive clustering as compared to other conventional clustering algorithms. The proposed adaptive fuzzy-K-means clustering algorithm provides better segmentation and visual quality for a range of images.


### 2.2. The Proposed Methodology

In this method, the face location is determined automatically by using an optimized K-means algorithm. First, the input image is reshaped into vectors of pixels values, followed by clustering the pixels, based on applying a certain threshold determined by the algorithm, into two classes. Pixels with values below the threshold are put in the nonface class. The rest of the pixels are put in the face class. However, some unwanted parts may get clustered to this face class. Therefore, the algorithm is applied again to the face class to obtain only the face part.


#### 2.2.1. Optimized K-Means Algorithm

The proposed K-means modified algorithm is intended to solve the limitations of the standard version by using differential equations to determine the optimum separation point. This algorithm finds the optimal centers for each cluster which corresponds to the global minimum of the k-means cluster. As an example, say we would like to cluster an image into two subclasses: face and nonface. So we look for a separator which will separate the image into two different clusters in two dimensions. Then the means of the clusters can be found easily. Figure2 shows the separation points and the mean of each class. To find these points we proposed the following modification in continuous and discrete cases.Figure 2
The separation points in two dimensions.Let{
x
1
,
x
2
,
x
3
,
…
,
x
N
} be a dataset and let x
i be an observation. Let f
(
x
) be the probability mass function (PMF) of x
(6)
f
(
x
)
=
∑
-
∞
∞
‍
1
N
δ
(
x
-
x
i
)
,
x
i
∈
R
,
where δ
(
x
) is the Dirac function.When the size of the data is very large, the question of the cost of the minimization of∑
i
=
1
k
‍
∑
x
i
∈
s
i
‍
|
x
i
-
μ
i
|
2 arises, where k is the number of clusters and μ
i is the mean of cluster s
i.Let(7)
Var_
k
=
∑
i
=
1
k
‍
∑
x
i
∈
s
i
‍
|
x
i
-
μ
i
|
2
.Then it is enough to minimizeVar_1 without losing the generality since
(8)
min
⁡
s
i
∑
i
=
1
d
‍
∑
x
j
∈
s
i
‍
|
x
j
-
μ
i
|
2
≥
∑
l
=
1
d
‍
min
⁡
s
i
(
∑
i
=
1
k
‍
∑
x
j
∈
s
i
‍
|
x
j
l
-
μ
i
l
|
2
)
,
where d is the dimension of x
j
=
(
x
j
1
,
…
,
x
j
d
) and μ
j
=
(
μ
j
1
,
…
,
μ
j
d
).


#### 2.2.2. Continuous Case

For the continuous casef
(
x
) is like the probability density function (PDF). Then for 2 classes, we need to minimize the following equation:
(9)
Var_2
=
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
.Ifd
≥
2 we can use (8)
(10)
min
⁡
Var_2
≥
∑
k
=
1
d
‍
min
⁡
(
∫
-
∞
M
‍
(
x
k
-
μ
1
k
)
2
f
(
x
k
)
d
x
+
∫
-
∞
M
‍
(
x
k
-
μ
2
k
)
2
f
(
x
k
)
d
x
)
.Letx
∈
R be any random variable with probability density function f
(
x
); we want to find μ
1 and μ
2 such that it minimizes (7) as follows:
(11)
(
μ
1
*
,
μ
2
*
)
=
arg
(
μ
1
,
μ
2
)
⁡
[
min
⁡
M
,
μ
1
,
μ
2
[
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
]
]
,
(12)
min
⁡
M
,
μ
1
,
μ
2
Var_2
=
min
⁡
M
[
min
⁡
μ
1
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
min
⁡
μ
2
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
]
,
(13)
min
⁡
M
,
μ
1
,
μ
2
Var_2
=
min
⁡
M
[
min
⁡
μ
1
∫
-
∞
M
‍
(
x
-
μ
1
(
M
)
)
2
f
(
x
)
d
x
+
min
⁡
μ
2
∫
M
∞
‍
(
x
-
μ
2
(
M
)
)
2
f
(
x
)
d
x
]
.We know that(14)
min
⁡
x
E
[
(
X
-
x
)
2
]
=
E
[
(
X
-
E
(
x
)
)
2
]
=
Var
⁡
(
x
)
.So we can findμ
1 and μ
2 in terms of M. Let us define Var
⁡
1 and Var
⁡
2 as follows:
(15)
Var
⁡
1
=
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
,
Var
⁡
2
=
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
.After simplification we get(16)
Var
⁡
1
=
∫
-
∞
M
‍
x
2
f
(
x
)
d
x
+
μ
1
2
∫
-
∞
M
‍
f
(
x
)
d
x
-
2
μ
1
∫
-
∞
M
‍
x
f
(
x
)
d
x
,
Var
⁡
2
=
∫
M
∞
‍
x
2
f
(
x
)
d
x
+
μ
2
2
∫
M
∞
‍
f
(
x
)
d
x
-
2
μ
2
∫
M
∞
‍
x
f
(
x
)
d
x
.To find the minimum,Var
⁡
1 and Var
⁡
2 are both partially differentiated with respect to μ
1 and μ
2, respectively. After simplification, we conclude that the minimization occurs when
(17)
μ
1
=
E
(
X
-
)
P
(
X
-
)
,
μ
2
=
E
(
X
+
)
P
(
X
+
)
,
where X
-
=
X
·
1
x
<
M and X
+
=
X
·
1
x
≥
M.To find the minimum ofVar_2, we need to find ∂
μ
/
∂
M
(18)
μ
1
=
(
∫
-
∞
M
‍
x
f
(
x
)
d
x
∫
-
∞
M
‍
f
(
x
)
d
x
)
⟹
∂
μ
1
∂
M
=
M
f
(
M
)
P
(
X
-
)
-
f
(
M
)
E
(
X
-
)
P
2
(
X
-
)
,
μ
2
=
(
∫
M
∞
‍
x
f
(
x
)
d
x
∫
M
∞
‍
f
(
x
)
d
x
)
⟹
∂
μ
2
∂
M
=
M
f
(
M
)
P
(
X
+
)
-
f
(
M
)
E
(
X
+
)
P
2
(
X
+
)
.After simplification we get(19)
∂
μ
1
∂
M
=
f
(
M
)
P
2
(
X
-
)
[
M
P
(
X
-
)
-
E
(
X
-
)
]
,
∂
μ
2
∂
M
=
f
(
M
)
P
2
(
X
+
)
[
M
P
(
X
+
)
-
E
(
X
+
)
]
.In order to find the minimum ofVar_2, we will find ∂
Var
⁡
1
/
∂
M and ∂
V
a
r
2
/
∂
M separately and add them as follows:
(20)
∂
∂
M
(
Var_2
)
=
∂
∂
M
(
Var
⁡
1
)
+
∂
∂
M
(
Var
⁡
2
)
∂
Var
⁡
1
∂
M
=
∂
∂
M
(
∫
-
∞
M
‍
x
2
f
(
x
)
d
x
+
μ
1
2
∫
-
∞
M
‍
f
(
x
)
d
x
-
2
μ
1
∫
-
∞
M
‍
x
f
(
x
)
d
x
)
=
M
2
f
(
M
)
+
2
μ
1
μ
1
′
P
(
X
-
)
+
μ
1
2
f
(
M
)
-
2
μ
1
′
E
(
X
-
)
-
2
μ
1
M
f
(
M
)
.After simplification we get(21)
∂
Var
⁡
1
∂
M
=
f
(
M
)
[
M
2
+
E
2
(
X
-
)
P
(
X
-
)
2
-
2
M
E
(
X
-
)
P
(
X
-
)
]
.Similarly(22)
∂
Var
⁡
2
∂
M
=
-
M
2
f
(
M
)
+
2
μ
2
μ
2
′
P
(
X
+
)
-
μ
2
2
f
(
M
)
-
2
μ
2
′
E
(
X
+
)
+
2
μ
2
M
f
(
M
)
.After simplification we get(23)
∂
Var
⁡
2
∂
M
=
f
(
M
)
[
-
M
2
+
E
2
(
X
+
)
P
(
X
+
)
2
-
2
M
E
(
X
+
)
P
(
X
+
)
]
.Adding both these differentials and after simplifying, it follows that(24)
∂
∂
M
(
Var_2
)
=
f
(
M
)
[
(
E
2
(
X
-
)
P
(
X
-
)
2
-
E
2
(
X
+
)
P
(
X
+
)
2
)
-
2
M
(
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
)
]
.Equating this with 0 gives(25)
∂
∂
M
(
Var_2
)
=
0
,
[
(
E
2
(
X
-
)
P
(
X
-
)
2
-
E
2
(
X
+
)
P
(
X
+
)
2
)
-
2
M
(
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
)
]
=
0
,
[
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
]
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
-
2
M
]
=
0
.But[
(
E
(
X
-
)
/
P
(
X
-
)
)
-
(
E
(
X
+
)
/
P
(
X
+
)
)
]
≠
0, as in that case μ
1
=
μ
2, which contradicts the initial assumption. Therefore
(26)
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
-
2
M
]
=
0
⟹
M
=
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
]
·
1
2
,
M
=
μ
1
+
μ
2
2
.Therefore, we have to find all values ofM that satisfy (26), which is easier than minimizing Var_2 directly.To find the minimum in terms ofM, it is enough to derive
(27)
∂
Var_2
(
μ
1
(
M
)
,
μ
2
(
M
)
)
∂
M
=
f
(
M
)
[
E
2
(
X
-
)
P
(
X
<
M
)
2
-
E
2
(
X
+
)
P
(
X
≥
M
)
2
-
2
M
[
E
(
X
-
)
P
(
X
≤
M
)
-
E
(
X
+
)
P
(
X
>
M
)
]
]
.Then we find(28)
M
=
[
E
(
X
+
)
P
(
X
+
)
+
E
(
X
-
)
P
(
X
-
)
]
·
1
2
.


#### 2.2.3. Finding the Centers of Some Common Probability Density Functions

(
1
)
Uniform Distribution. The probability density function of the continuous uniform distribution is
(29)
f
(
x
)
=
{
1
b
-
a
,
a
<
x
<
b
0
otherwise
,
where a and b are the two parameters of the distribution. The expected values for X
·
1
x
≥
M and X
·
1
x
<
M are
(30)
E
(
X
+
)
=
1
2
(
b
2
-
M
2
b
-
a
)
,
E
(
X
-
)
=
1
2
(
M
2
-
a
2
b
-
a
)
.Putting all these in (28) and solving, we get
(31)
M
=
1
2
(
a
+
b
)
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
1
4
(
3
a
+
b
)
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
1
4
(
a
+
3
b
)
.(
2
)
Log-Normal Distribution. The probability density function of a log-normal distribution is
(32)
f
(
x
;
μ
,
σ
)
=
1
x
σ
2
π
e
-
(
ln
⁡
x
-
μ
)
2
/
2
σ
2
,
x
>
0
.The probabilities in left and right ofM are, respectively,
(33)
P
(
X
+
)
=
1
2
-
1
2
erf
⁡
(
1
2
2
(
ln
⁡
M
-
μ
)
σ
)
,
P
(
X
-
)
=
1
2
+
1
2
erf
⁡
(
1
2
2
(
ln
⁡
M
-
μ
)
σ
)
.The expected values forX
·
1
x
≥
M and X
·
1
x
<
M are
(34)
E
(
X
+
)
=
[
1
2
e
μ
+
(
1
/
2
)
σ
2
{
1
+
erf
⁡
(
1
2
(
μ
+
σ
2
-
ln
⁡
M
)
2
σ
)
}
]
,
E
(
X
-
)
=
[
1
2
e
μ
+
(
1
/
2
)
σ
2
{
1
+
erf
⁡
(
1
2
(
-
μ
-
σ
2
+
ln
⁡
M
)
2
σ
)
}
]
.Putting all these in (28) and assuming μ
=
0 and σ
=
1 and solving, we get
(35)
M
=
4.244942583
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
1.196827816
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
7.293057348
.(
3
)
Normal Distribution. The probability density function of a normal distribution is
(36)
f
(
x
,
μ
,
σ
)
=
1
σ
2
π
e
-
(
x
-
μ
)
2
/
2
σ
2
.The probabilities in left and right ofM are, respectively,
(37)
P
(
X
+
)
=
1
2
-
1
2
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
,
P
(
X
-
)
=
1
2
+
1
2
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
.The expected values forX
·
1
x
≥
M and X
·
1
x
<
M are
(38)
E
(
X
+
)
=
-
1
2
μ
{
-
1
+
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
}
,
E
(
X
-
)
=
1
2
μ
{
1
+
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
}
.Putting all these in (28) and solving we get M
=
μ. Assuming μ
=
0 and σ
=
1 and solving, we get
(39)
M
=
0
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
-
2
π
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
2
π
.


#### 2.2.4. Discrete Case

LetX be a discrete random variable and assume we have a set of N observations from X. Let p
j be the mass density function for an observation x
j.Letx
i
-
1
<
M
i
-
1
<
x
i and x
i
<
M
i
<
x
i
+
1. We define μ
x
≤
x
i as the mean for all x
≤
x
i, and μ
x
≥
x
i as the mean for all x
≥
x
i. Define Var_2
(
M
i
-
1
) as the variance of two centers forced by M
i
-
1 as a separator
(40)
Var
⁡
(
M
i
-
1
)
=
∑
x
j
≤
x
i
-
1
‍
(
x
j
-
μ
x
j
≤
x
i
-
1
)
2
p
j
+
∑
x
j
≥
x
i
‍
(
x
j
-
μ
x
j
≥
x
i
)
2
p
j
,
(41)
Var
⁡
(
M
i
)
=
∑
x
j
≤
x
i
‍
(
x
j
-
μ
x
j
≤
x
i
)
2
p
j
+
∑
x
j
≥
x
i
+
1
‍
(
x
j
-
μ
x
j
≥
x
i
+
1
)
2
p
j
.The first part ofVar
⁡
(
M
i
-
1
) simplifies as follows:
(42)
∑
x
j
≤
x
i
-
1
‍
(
x
j
-
μ
x
j
≤
x
i
-
1
)
2
p
j
=
∑
x
j
≤
x
i
-
1
‍
x
j
2
p
j
-
2
μ
x
j
≤
x
i
-
1
×
∑
x
j
≤
x
i
-
1
‍
x
j
p
j
+
μ
x
j
≤
x
i
-
1
2
∑
x
j
≤
x
i
-
1
‍
p
j
=
E
(
X
2
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
.Similarly,(43)
∑
x
j
≥
x
i
‍
(
x
j
-
μ
x
j
≥
x
i
)
2
p
j
=
E
(
X
2
·
1
x
j
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
,
Var_2
(
M
i
-
1
)
=
E
(
X
2
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
2
·
1
x
j
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
,
Var_2
(
M
i
)
=
E
(
X
2
·
1
x
j
≤
x
i
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
+
E
(
X
2
·
1
x
j
≥
x
i
+
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
.
Define Δ
Var_2
=
Var_2
(
M
i
)
-
Var_2
(
M
i
-
1
)
(44)
Δ
Var_2
=
[
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
[
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
]
.Now, in order to simplify the calculation, we rearrangeΔ
Var_2 as follows:
(45)
Δ
Var_2
=
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
]
+
[
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
.The first part is simplified as follows:(46)
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
=
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
+
x
i
p
i
]
2
P
(
X
≤
x
i
-
1
)
+
p
i
=
1
P
(
X
≤
x
i
-
1
)
[
P
(
X
≤
x
i
-
1
)
+
p
i
]
×
{
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
)
+
p
i
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
2
x
i
p
i
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
x
i
2
p
i
2
P
(
X
≤
x
i
-
1
)
}
=
1
P
(
X
≤
x
i
-
1
)
[
P
(
X
≤
x
i
-
1
)
+
p
i
]
×
{
p
i
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
-
2
x
i
p
i
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
x
i
2
p
i
2
P
(
X
≤
x
i
-
1
)
}
.Similarly, simplifying the second part yields(47)
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
=
[
E
(
X
·
1
x
j
≥
x
i
+
1
)
+
x
i
p
i
]
2
P
(
X
≥
x
i
+
1
)
+
p
i
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
=
1
P
(
X
≥
x
i
+
1
)
[
P
(
X
≥
x
i
+
1
)
+
p
i
]
×
{
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
-
p
i
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
2
x
i
p
i
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
x
i
2
p
i
2
P
(
X
≥
x
i
+
1
)
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
}
=
1
P
(
X
≥
x
i
+
1
)
[
P
(
X
≥
x
i
+
1
)
+
p
i
]
×
{
-
p
i
+
1
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
+
2
x
i
p
i
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
x
i
2
p
i
+
1
2
P
(
X
≥
x
i
+
1
)
-
p
i
+
1
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
}
.In general, these types of optimization arise when the dataset is very big.If we considerN to be very large then p
k
≪
max
⁡
(
∑
i
≥
k
+
1
‍
P
i
;
∑
i
≤
k
-
1
‍
P
i
), P
(
X
≤
x
i
-
1
)
+
p
i
≈
P
(
X
≤
x
i
), and P
(
X
≥
x
i
+
1
)
+
p
i
≈
P
(
X
≥
x
i
+
1
):
(48)
Δ
Var_2
≈
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
p
i
P
2
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
p
i
P
2
(
X
≥
x
i
+
1
)
]
-
2
p
i
x
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
p
i
2
x
i
2
[
1
P
(
X
≤
x
i
-
1
)
-
1
P
(
X
≥
x
i
+
1
)
]
and, as lim
⁡
N
→
∞
{
p
i
2
x
i
2
[
(
1
/
P
(
X
≤
x
i
-
1
)
)
-
(
1
/
P
(
X
≥
x
i
+
1
)
)
]
}
=
0, we can further simplify Δ
Var_2 as follows:
(49)
Δ
Var_2
≈
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
p
i
P
2
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
p
i
P
2
(
X
≥
x
i
+
1
)
]
-
2
p
i
x
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
=
p
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
×
[
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
)
-
2
x
i
]
.To find the minimum ofΔ
Var_2, let us locate when Δ
Var_2
=
0:
(50)
⟹
p
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
×
[
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
)
-
2
x
i
]
=
0
.Butp
i
≠
0, and so is [
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
/
P
(
X
≤
x
i
-
1
)
)
-
(
E
(
X
·
1
x
j
≥
x
i
+
1
)
/
P
(
X
≥
x
i
+
1
)
)
]
<
0, since x
i is an increasing sequence, which implies that μ
i
-
1
<
μ
i where
(51)
μ
i
-
1
=
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
,
μ
i
=
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
.Therefore,(52)
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
2
x
i
=
0
.Let us define(53)
F
(
x
i
)
=
2
x
i
-
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
F
(
x
i
)
=
0
⟹
x
i
=
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
·
1
2
⟹
x
i
=
μ
i
-
1
+
μ
i
2
,
where μ
i
-
1 and μ
i are the means of the first and second parts, respectively.To find the minimum for the total variation, it is enough to find the good separatorM between x
i and x
i
-
1 such that F
(
x
i
)
·
F
(
x
i
+
1
)
<
0 and F
(
x
i
)
<
0 (see Figure 3).Figure 3
Discrete random variables.After finding few local minima, we can get the global minimum by finding the smallest among them.


## 2.1. Modifications in K-Means Clustering

Global k-means algorithm [34] is a proposal to improve global search properties of k-means algorithm and its performance on very large datasets by computing clusters successively. The main weakness of k-means clustering, that is, its sensitivity to initial positions of the cluster centers, is overcome by global k-means clustering algorithm which consists of a deterministic global optimization technique. It does not select initial values randomly; instead an incremental approach is applied to optimally add one new cluster center at each stage. Global k-means also proposes a method to reduce the computational load while maintain the solution quality. Experiments were performed to compare the performance of k-means, global k-means, min k-means, and fast global k-means as shown in Figure 1. Numerical results show that the global k-means algorithm considerably outperforms other k-means algorithms.Figure 1
Comparison of k-means.Bagirov [35] proposed a new version of the global k-means algorithm for minimum sum-of-squares clustering problems. He also compared three different versions of the k-means algorithm to propose the modified version of the global k-means algorithm. The proposed algorithm computes clusters incrementally and k
-
1 cluster centers from the previous iteration are used to compute k-partition of a dataset. The starting point of the kth cluster center is computed by minimizing an auxiliary cluster function. Given a finite set of points A in the n-dimensional space the global k-means compute the centroid X
1 of the set A as the following equation:
(4)
X
1
=
1
m
∑
i
=
1
m
‍
a
i
,
a
i
∈
A
,
i
=
1
,
…
,
m
.Numerical experiments performed on 14 datasets demonstrate the efficiency of the modified global k-means algorithm in comparison to the multistart k-means (MS k-means) and the GKM algorithms when the number of clustersk
>
5. Modified global k-means however requires more computational time than the global k-means algorithm. Xie and Jiang [36] also proposed a novel version of the global k-means algorithm. The method of creating the next cluster center in the global k-means algorithm is modified and that showed a better execution time without affecting the performance of the global k-means algorithm. Another extension of the standard k-means clustering is the Global Kernel k-means [37] which optimizes the clustering error in feature space by employing kernel-means as a local search procedure. A kernel function is used in order to solve the M-clustering problem and near-optimal solutions are used to optimize the clustering error in the feature space. This incremental algorithm has high ability to identify nonlinearly separable clusters in input space and no dependence on cluster initialization. It can handle weighted data points making it suitable for graph partitioning applications. Two major modifications were performed in order to reduce the computational cost with no major effect on the solution quality.Video imaging and Image segmentation are important applications for k-means clustering. Hung et al. [38] modified the k-means algorithm for color image segmentation where a weight selection procedure in the W-k-means algorithm is proposed. The evaluation function F
(
I
) is used for comparison:
(5)
F
(
I
)
=
R
∑
i
=
1
R
‍
(
e
i
/
256
)
2
A
i
,
where I = segmented image, R = the number of regions in the segmented image, A
i = the area, or the number of pixels of the ith region, and e
i = the color error of region i.e
i is the sum of the Euclidean distance of the color vectors between the original image and the segmented image. Smaller values of F
(
I
) give better segmentation results. Results from color image segmentation illustrate that the proposed procedure produces better segmentation than the random initialization. Maliatski and Yadid-Pecht [39] also propose an adaptive k-means clustering, hardware-driven approach. The algorithm uses both pixel intensity and pixel distance in the clustering process. In [40] also a FPGA implementation of real-time k-means clustering is done for colored images. A filtering algorithm is used for hardware implementation. Suliman and Isa [11] presented a unique clustering algorithm called adaptive fuzzy-K-means clustering for image segmentation. It can be used for different types of images including medical images, microscopic images, and CCD camera images. The concepts of fuzziness and belongingness are applied to produce more adaptive clustering as compared to other conventional clustering algorithms. The proposed adaptive fuzzy-K-means clustering algorithm provides better segmentation and visual quality for a range of images.


## 2.2. The Proposed Methodology

In this method, the face location is determined automatically by using an optimized K-means algorithm. First, the input image is reshaped into vectors of pixels values, followed by clustering the pixels, based on applying a certain threshold determined by the algorithm, into two classes. Pixels with values below the threshold are put in the nonface class. The rest of the pixels are put in the face class. However, some unwanted parts may get clustered to this face class. Therefore, the algorithm is applied again to the face class to obtain only the face part.


### 2.2.1. Optimized K-Means Algorithm

The proposed K-means modified algorithm is intended to solve the limitations of the standard version by using differential equations to determine the optimum separation point. This algorithm finds the optimal centers for each cluster which corresponds to the global minimum of the k-means cluster. As an example, say we would like to cluster an image into two subclasses: face and nonface. So we look for a separator which will separate the image into two different clusters in two dimensions. Then the means of the clusters can be found easily. Figure2 shows the separation points and the mean of each class. To find these points we proposed the following modification in continuous and discrete cases.Figure 2
The separation points in two dimensions.Let{
x
1
,
x
2
,
x
3
,
…
,
x
N
} be a dataset and let x
i be an observation. Let f
(
x
) be the probability mass function (PMF) of x
(6)
f
(
x
)
=
∑
-
∞
∞
‍
1
N
δ
(
x
-
x
i
)
,
x
i
∈
R
,
where δ
(
x
) is the Dirac function.When the size of the data is very large, the question of the cost of the minimization of∑
i
=
1
k
‍
∑
x
i
∈
s
i
‍
|
x
i
-
μ
i
|
2 arises, where k is the number of clusters and μ
i is the mean of cluster s
i.Let(7)
Var_
k
=
∑
i
=
1
k
‍
∑
x
i
∈
s
i
‍
|
x
i
-
μ
i
|
2
.Then it is enough to minimizeVar_1 without losing the generality since
(8)
min
⁡
s
i
∑
i
=
1
d
‍
∑
x
j
∈
s
i
‍
|
x
j
-
μ
i
|
2
≥
∑
l
=
1
d
‍
min
⁡
s
i
(
∑
i
=
1
k
‍
∑
x
j
∈
s
i
‍
|
x
j
l
-
μ
i
l
|
2
)
,
where d is the dimension of x
j
=
(
x
j
1
,
…
,
x
j
d
) and μ
j
=
(
μ
j
1
,
…
,
μ
j
d
).


### 2.2.2. Continuous Case

For the continuous casef
(
x
) is like the probability density function (PDF). Then for 2 classes, we need to minimize the following equation:
(9)
Var_2
=
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
.Ifd
≥
2 we can use (8)
(10)
min
⁡
Var_2
≥
∑
k
=
1
d
‍
min
⁡
(
∫
-
∞
M
‍
(
x
k
-
μ
1
k
)
2
f
(
x
k
)
d
x
+
∫
-
∞
M
‍
(
x
k
-
μ
2
k
)
2
f
(
x
k
)
d
x
)
.Letx
∈
R be any random variable with probability density function f
(
x
); we want to find μ
1 and μ
2 such that it minimizes (7) as follows:
(11)
(
μ
1
*
,
μ
2
*
)
=
arg
(
μ
1
,
μ
2
)
⁡
[
min
⁡
M
,
μ
1
,
μ
2
[
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
]
]
,
(12)
min
⁡
M
,
μ
1
,
μ
2
Var_2
=
min
⁡
M
[
min
⁡
μ
1
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
min
⁡
μ
2
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
]
,
(13)
min
⁡
M
,
μ
1
,
μ
2
Var_2
=
min
⁡
M
[
min
⁡
μ
1
∫
-
∞
M
‍
(
x
-
μ
1
(
M
)
)
2
f
(
x
)
d
x
+
min
⁡
μ
2
∫
M
∞
‍
(
x
-
μ
2
(
M
)
)
2
f
(
x
)
d
x
]
.We know that(14)
min
⁡
x
E
[
(
X
-
x
)
2
]
=
E
[
(
X
-
E
(
x
)
)
2
]
=
Var
⁡
(
x
)
.So we can findμ
1 and μ
2 in terms of M. Let us define Var
⁡
1 and Var
⁡
2 as follows:
(15)
Var
⁡
1
=
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
,
Var
⁡
2
=
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
.After simplification we get(16)
Var
⁡
1
=
∫
-
∞
M
‍
x
2
f
(
x
)
d
x
+
μ
1
2
∫
-
∞
M
‍
f
(
x
)
d
x
-
2
μ
1
∫
-
∞
M
‍
x
f
(
x
)
d
x
,
Var
⁡
2
=
∫
M
∞
‍
x
2
f
(
x
)
d
x
+
μ
2
2
∫
M
∞
‍
f
(
x
)
d
x
-
2
μ
2
∫
M
∞
‍
x
f
(
x
)
d
x
.To find the minimum,Var
⁡
1 and Var
⁡
2 are both partially differentiated with respect to μ
1 and μ
2, respectively. After simplification, we conclude that the minimization occurs when
(17)
μ
1
=
E
(
X
-
)
P
(
X
-
)
,
μ
2
=
E
(
X
+
)
P
(
X
+
)
,
where X
-
=
X
·
1
x
<
M and X
+
=
X
·
1
x
≥
M.To find the minimum ofVar_2, we need to find ∂
μ
/
∂
M
(18)
μ
1
=
(
∫
-
∞
M
‍
x
f
(
x
)
d
x
∫
-
∞
M
‍
f
(
x
)
d
x
)
⟹
∂
μ
1
∂
M
=
M
f
(
M
)
P
(
X
-
)
-
f
(
M
)
E
(
X
-
)
P
2
(
X
-
)
,
μ
2
=
(
∫
M
∞
‍
x
f
(
x
)
d
x
∫
M
∞
‍
f
(
x
)
d
x
)
⟹
∂
μ
2
∂
M
=
M
f
(
M
)
P
(
X
+
)
-
f
(
M
)
E
(
X
+
)
P
2
(
X
+
)
.After simplification we get(19)
∂
μ
1
∂
M
=
f
(
M
)
P
2
(
X
-
)
[
M
P
(
X
-
)
-
E
(
X
-
)
]
,
∂
μ
2
∂
M
=
f
(
M
)
P
2
(
X
+
)
[
M
P
(
X
+
)
-
E
(
X
+
)
]
.In order to find the minimum ofVar_2, we will find ∂
Var
⁡
1
/
∂
M and ∂
V
a
r
2
/
∂
M separately and add them as follows:
(20)
∂
∂
M
(
Var_2
)
=
∂
∂
M
(
Var
⁡
1
)
+
∂
∂
M
(
Var
⁡
2
)
∂
Var
⁡
1
∂
M
=
∂
∂
M
(
∫
-
∞
M
‍
x
2
f
(
x
)
d
x
+
μ
1
2
∫
-
∞
M
‍
f
(
x
)
d
x
-
2
μ
1
∫
-
∞
M
‍
x
f
(
x
)
d
x
)
=
M
2
f
(
M
)
+
2
μ
1
μ
1
′
P
(
X
-
)
+
μ
1
2
f
(
M
)
-
2
μ
1
′
E
(
X
-
)
-
2
μ
1
M
f
(
M
)
.After simplification we get(21)
∂
Var
⁡
1
∂
M
=
f
(
M
)
[
M
2
+
E
2
(
X
-
)
P
(
X
-
)
2
-
2
M
E
(
X
-
)
P
(
X
-
)
]
.Similarly(22)
∂
Var
⁡
2
∂
M
=
-
M
2
f
(
M
)
+
2
μ
2
μ
2
′
P
(
X
+
)
-
μ
2
2
f
(
M
)
-
2
μ
2
′
E
(
X
+
)
+
2
μ
2
M
f
(
M
)
.After simplification we get(23)
∂
Var
⁡
2
∂
M
=
f
(
M
)
[
-
M
2
+
E
2
(
X
+
)
P
(
X
+
)
2
-
2
M
E
(
X
+
)
P
(
X
+
)
]
.Adding both these differentials and after simplifying, it follows that(24)
∂
∂
M
(
Var_2
)
=
f
(
M
)
[
(
E
2
(
X
-
)
P
(
X
-
)
2
-
E
2
(
X
+
)
P
(
X
+
)
2
)
-
2
M
(
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
)
]
.Equating this with 0 gives(25)
∂
∂
M
(
Var_2
)
=
0
,
[
(
E
2
(
X
-
)
P
(
X
-
)
2
-
E
2
(
X
+
)
P
(
X
+
)
2
)
-
2
M
(
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
)
]
=
0
,
[
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
]
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
-
2
M
]
=
0
.But[
(
E
(
X
-
)
/
P
(
X
-
)
)
-
(
E
(
X
+
)
/
P
(
X
+
)
)
]
≠
0, as in that case μ
1
=
μ
2, which contradicts the initial assumption. Therefore
(26)
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
-
2
M
]
=
0
⟹
M
=
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
]
·
1
2
,
M
=
μ
1
+
μ
2
2
.Therefore, we have to find all values ofM that satisfy (26), which is easier than minimizing Var_2 directly.To find the minimum in terms ofM, it is enough to derive
(27)
∂
Var_2
(
μ
1
(
M
)
,
μ
2
(
M
)
)
∂
M
=
f
(
M
)
[
E
2
(
X
-
)
P
(
X
<
M
)
2
-
E
2
(
X
+
)
P
(
X
≥
M
)
2
-
2
M
[
E
(
X
-
)
P
(
X
≤
M
)
-
E
(
X
+
)
P
(
X
>
M
)
]
]
.Then we find(28)
M
=
[
E
(
X
+
)
P
(
X
+
)
+
E
(
X
-
)
P
(
X
-
)
]
·
1
2
.


### 2.2.3. Finding the Centers of Some Common Probability Density Functions

(
1
)
Uniform Distribution. The probability density function of the continuous uniform distribution is
(29)
f
(
x
)
=
{
1
b
-
a
,
a
<
x
<
b
0
otherwise
,
where a and b are the two parameters of the distribution. The expected values for X
·
1
x
≥
M and X
·
1
x
<
M are
(30)
E
(
X
+
)
=
1
2
(
b
2
-
M
2
b
-
a
)
,
E
(
X
-
)
=
1
2
(
M
2
-
a
2
b
-
a
)
.Putting all these in (28) and solving, we get
(31)
M
=
1
2
(
a
+
b
)
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
1
4
(
3
a
+
b
)
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
1
4
(
a
+
3
b
)
.(
2
)
Log-Normal Distribution. The probability density function of a log-normal distribution is
(32)
f
(
x
;
μ
,
σ
)
=
1
x
σ
2
π
e
-
(
ln
⁡
x
-
μ
)
2
/
2
σ
2
,
x
>
0
.The probabilities in left and right ofM are, respectively,
(33)
P
(
X
+
)
=
1
2
-
1
2
erf
⁡
(
1
2
2
(
ln
⁡
M
-
μ
)
σ
)
,
P
(
X
-
)
=
1
2
+
1
2
erf
⁡
(
1
2
2
(
ln
⁡
M
-
μ
)
σ
)
.The expected values forX
·
1
x
≥
M and X
·
1
x
<
M are
(34)
E
(
X
+
)
=
[
1
2
e
μ
+
(
1
/
2
)
σ
2
{
1
+
erf
⁡
(
1
2
(
μ
+
σ
2
-
ln
⁡
M
)
2
σ
)
}
]
,
E
(
X
-
)
=
[
1
2
e
μ
+
(
1
/
2
)
σ
2
{
1
+
erf
⁡
(
1
2
(
-
μ
-
σ
2
+
ln
⁡
M
)
2
σ
)
}
]
.Putting all these in (28) and assuming μ
=
0 and σ
=
1 and solving, we get
(35)
M
=
4.244942583
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
1.196827816
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
7.293057348
.(
3
)
Normal Distribution. The probability density function of a normal distribution is
(36)
f
(
x
,
μ
,
σ
)
=
1
σ
2
π
e
-
(
x
-
μ
)
2
/
2
σ
2
.The probabilities in left and right ofM are, respectively,
(37)
P
(
X
+
)
=
1
2
-
1
2
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
,
P
(
X
-
)
=
1
2
+
1
2
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
.The expected values forX
·
1
x
≥
M and X
·
1
x
<
M are
(38)
E
(
X
+
)
=
-
1
2
μ
{
-
1
+
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
}
,
E
(
X
-
)
=
1
2
μ
{
1
+
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
}
.Putting all these in (28) and solving we get M
=
μ. Assuming μ
=
0 and σ
=
1 and solving, we get
(39)
M
=
0
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
-
2
π
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
2
π
.


### 2.2.4. Discrete Case

LetX be a discrete random variable and assume we have a set of N observations from X. Let p
j be the mass density function for an observation x
j.Letx
i
-
1
<
M
i
-
1
<
x
i and x
i
<
M
i
<
x
i
+
1. We define μ
x
≤
x
i as the mean for all x
≤
x
i, and μ
x
≥
x
i as the mean for all x
≥
x
i. Define Var_2
(
M
i
-
1
) as the variance of two centers forced by M
i
-
1 as a separator
(40)
Var
⁡
(
M
i
-
1
)
=
∑
x
j
≤
x
i
-
1
‍
(
x
j
-
μ
x
j
≤
x
i
-
1
)
2
p
j
+
∑
x
j
≥
x
i
‍
(
x
j
-
μ
x
j
≥
x
i
)
2
p
j
,
(41)
Var
⁡
(
M
i
)
=
∑
x
j
≤
x
i
‍
(
x
j
-
μ
x
j
≤
x
i
)
2
p
j
+
∑
x
j
≥
x
i
+
1
‍
(
x
j
-
μ
x
j
≥
x
i
+
1
)
2
p
j
.The first part ofVar
⁡
(
M
i
-
1
) simplifies as follows:
(42)
∑
x
j
≤
x
i
-
1
‍
(
x
j
-
μ
x
j
≤
x
i
-
1
)
2
p
j
=
∑
x
j
≤
x
i
-
1
‍
x
j
2
p
j
-
2
μ
x
j
≤
x
i
-
1
×
∑
x
j
≤
x
i
-
1
‍
x
j
p
j
+
μ
x
j
≤
x
i
-
1
2
∑
x
j
≤
x
i
-
1
‍
p
j
=
E
(
X
2
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
.Similarly,(43)
∑
x
j
≥
x
i
‍
(
x
j
-
μ
x
j
≥
x
i
)
2
p
j
=
E
(
X
2
·
1
x
j
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
,
Var_2
(
M
i
-
1
)
=
E
(
X
2
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
2
·
1
x
j
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
,
Var_2
(
M
i
)
=
E
(
X
2
·
1
x
j
≤
x
i
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
+
E
(
X
2
·
1
x
j
≥
x
i
+
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
.
Define Δ
Var_2
=
Var_2
(
M
i
)
-
Var_2
(
M
i
-
1
)
(44)
Δ
Var_2
=
[
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
[
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
]
.Now, in order to simplify the calculation, we rearrangeΔ
Var_2 as follows:
(45)
Δ
Var_2
=
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
]
+
[
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
.The first part is simplified as follows:(46)
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
=
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
+
x
i
p
i
]
2
P
(
X
≤
x
i
-
1
)
+
p
i
=
1
P
(
X
≤
x
i
-
1
)
[
P
(
X
≤
x
i
-
1
)
+
p
i
]
×
{
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
)
+
p
i
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
2
x
i
p
i
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
x
i
2
p
i
2
P
(
X
≤
x
i
-
1
)
}
=
1
P
(
X
≤
x
i
-
1
)
[
P
(
X
≤
x
i
-
1
)
+
p
i
]
×
{
p
i
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
-
2
x
i
p
i
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
x
i
2
p
i
2
P
(
X
≤
x
i
-
1
)
}
.Similarly, simplifying the second part yields(47)
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
=
[
E
(
X
·
1
x
j
≥
x
i
+
1
)
+
x
i
p
i
]
2
P
(
X
≥
x
i
+
1
)
+
p
i
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
=
1
P
(
X
≥
x
i
+
1
)
[
P
(
X
≥
x
i
+
1
)
+
p
i
]
×
{
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
-
p
i
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
2
x
i
p
i
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
x
i
2
p
i
2
P
(
X
≥
x
i
+
1
)
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
}
=
1
P
(
X
≥
x
i
+
1
)
[
P
(
X
≥
x
i
+
1
)
+
p
i
]
×
{
-
p
i
+
1
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
+
2
x
i
p
i
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
x
i
2
p
i
+
1
2
P
(
X
≥
x
i
+
1
)
-
p
i
+
1
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
}
.In general, these types of optimization arise when the dataset is very big.If we considerN to be very large then p
k
≪
max
⁡
(
∑
i
≥
k
+
1
‍
P
i
;
∑
i
≤
k
-
1
‍
P
i
), P
(
X
≤
x
i
-
1
)
+
p
i
≈
P
(
X
≤
x
i
), and P
(
X
≥
x
i
+
1
)
+
p
i
≈
P
(
X
≥
x
i
+
1
):
(48)
Δ
Var_2
≈
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
p
i
P
2
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
p
i
P
2
(
X
≥
x
i
+
1
)
]
-
2
p
i
x
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
p
i
2
x
i
2
[
1
P
(
X
≤
x
i
-
1
)
-
1
P
(
X
≥
x
i
+
1
)
]
and, as lim
⁡
N
→
∞
{
p
i
2
x
i
2
[
(
1
/
P
(
X
≤
x
i
-
1
)
)
-
(
1
/
P
(
X
≥
x
i
+
1
)
)
]
}
=
0, we can further simplify Δ
Var_2 as follows:
(49)
Δ
Var_2
≈
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
p
i
P
2
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
p
i
P
2
(
X
≥
x
i
+
1
)
]
-
2
p
i
x
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
=
p
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
×
[
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
)
-
2
x
i
]
.To find the minimum ofΔ
Var_2, let us locate when Δ
Var_2
=
0:
(50)
⟹
p
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
×
[
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
)
-
2
x
i
]
=
0
.Butp
i
≠
0, and so is [
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
/
P
(
X
≤
x
i
-
1
)
)
-
(
E
(
X
·
1
x
j
≥
x
i
+
1
)
/
P
(
X
≥
x
i
+
1
)
)
]
<
0, since x
i is an increasing sequence, which implies that μ
i
-
1
<
μ
i where
(51)
μ
i
-
1
=
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
,
μ
i
=
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
.Therefore,(52)
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
2
x
i
=
0
.Let us define(53)
F
(
x
i
)
=
2
x
i
-
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
F
(
x
i
)
=
0
⟹
x
i
=
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
·
1
2
⟹
x
i
=
μ
i
-
1
+
μ
i
2
,
where μ
i
-
1 and μ
i are the means of the first and second parts, respectively.To find the minimum for the total variation, it is enough to find the good separatorM between x
i and x
i
-
1 such that F
(
x
i
)
·
F
(
x
i
+
1
)
<
0 and F
(
x
i
)
<
0 (see Figure 3).Figure 3
Discrete random variables.After finding few local minima, we can get the global minimum by finding the smallest among them.


## 2.2.1. Optimized K-Means Algorithm

The proposed K-means modified algorithm is intended to solve the limitations of the standard version by using differential equations to determine the optimum separation point. This algorithm finds the optimal centers for each cluster which corresponds to the global minimum of the k-means cluster. As an example, say we would like to cluster an image into two subclasses: face and nonface. So we look for a separator which will separate the image into two different clusters in two dimensions. Then the means of the clusters can be found easily. Figure2 shows the separation points and the mean of each class. To find these points we proposed the following modification in continuous and discrete cases.Figure 2
The separation points in two dimensions.Let{
x
1
,
x
2
,
x
3
,
…
,
x
N
} be a dataset and let x
i be an observation. Let f
(
x
) be the probability mass function (PMF) of x
(6)
f
(
x
)
=
∑
-
∞
∞
‍
1
N
δ
(
x
-
x
i
)
,
x
i
∈
R
,
where δ
(
x
) is the Dirac function.When the size of the data is very large, the question of the cost of the minimization of∑
i
=
1
k
‍
∑
x
i
∈
s
i
‍
|
x
i
-
μ
i
|
2 arises, where k is the number of clusters and μ
i is the mean of cluster s
i.Let(7)
Var_
k
=
∑
i
=
1
k
‍
∑
x
i
∈
s
i
‍
|
x
i
-
μ
i
|
2
.Then it is enough to minimizeVar_1 without losing the generality since
(8)
min
⁡
s
i
∑
i
=
1
d
‍
∑
x
j
∈
s
i
‍
|
x
j
-
μ
i
|
2
≥
∑
l
=
1
d
‍
min
⁡
s
i
(
∑
i
=
1
k
‍
∑
x
j
∈
s
i
‍
|
x
j
l
-
μ
i
l
|
2
)
,
where d is the dimension of x
j
=
(
x
j
1
,
…
,
x
j
d
) and μ
j
=
(
μ
j
1
,
…
,
μ
j
d
).


## 2.2.2. Continuous Case

For the continuous casef
(
x
) is like the probability density function (PDF). Then for 2 classes, we need to minimize the following equation:
(9)
Var_2
=
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
.Ifd
≥
2 we can use (8)
(10)
min
⁡
Var_2
≥
∑
k
=
1
d
‍
min
⁡
(
∫
-
∞
M
‍
(
x
k
-
μ
1
k
)
2
f
(
x
k
)
d
x
+
∫
-
∞
M
‍
(
x
k
-
μ
2
k
)
2
f
(
x
k
)
d
x
)
.Letx
∈
R be any random variable with probability density function f
(
x
); we want to find μ
1 and μ
2 such that it minimizes (7) as follows:
(11)
(
μ
1
*
,
μ
2
*
)
=
arg
(
μ
1
,
μ
2
)
⁡
[
min
⁡
M
,
μ
1
,
μ
2
[
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
]
]
,
(12)
min
⁡
M
,
μ
1
,
μ
2
Var_2
=
min
⁡
M
[
min
⁡
μ
1
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
+
min
⁡
μ
2
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
]
,
(13)
min
⁡
M
,
μ
1
,
μ
2
Var_2
=
min
⁡
M
[
min
⁡
μ
1
∫
-
∞
M
‍
(
x
-
μ
1
(
M
)
)
2
f
(
x
)
d
x
+
min
⁡
μ
2
∫
M
∞
‍
(
x
-
μ
2
(
M
)
)
2
f
(
x
)
d
x
]
.We know that(14)
min
⁡
x
E
[
(
X
-
x
)
2
]
=
E
[
(
X
-
E
(
x
)
)
2
]
=
Var
⁡
(
x
)
.So we can findμ
1 and μ
2 in terms of M. Let us define Var
⁡
1 and Var
⁡
2 as follows:
(15)
Var
⁡
1
=
∫
-
∞
M
‍
(
x
-
μ
1
)
2
f
(
x
)
d
x
,
Var
⁡
2
=
∫
M
∞
‍
(
x
-
μ
2
)
2
f
(
x
)
d
x
.After simplification we get(16)
Var
⁡
1
=
∫
-
∞
M
‍
x
2
f
(
x
)
d
x
+
μ
1
2
∫
-
∞
M
‍
f
(
x
)
d
x
-
2
μ
1
∫
-
∞
M
‍
x
f
(
x
)
d
x
,
Var
⁡
2
=
∫
M
∞
‍
x
2
f
(
x
)
d
x
+
μ
2
2
∫
M
∞
‍
f
(
x
)
d
x
-
2
μ
2
∫
M
∞
‍
x
f
(
x
)
d
x
.To find the minimum,Var
⁡
1 and Var
⁡
2 are both partially differentiated with respect to μ
1 and μ
2, respectively. After simplification, we conclude that the minimization occurs when
(17)
μ
1
=
E
(
X
-
)
P
(
X
-
)
,
μ
2
=
E
(
X
+
)
P
(
X
+
)
,
where X
-
=
X
·
1
x
<
M and X
+
=
X
·
1
x
≥
M.To find the minimum ofVar_2, we need to find ∂
μ
/
∂
M
(18)
μ
1
=
(
∫
-
∞
M
‍
x
f
(
x
)
d
x
∫
-
∞
M
‍
f
(
x
)
d
x
)
⟹
∂
μ
1
∂
M
=
M
f
(
M
)
P
(
X
-
)
-
f
(
M
)
E
(
X
-
)
P
2
(
X
-
)
,
μ
2
=
(
∫
M
∞
‍
x
f
(
x
)
d
x
∫
M
∞
‍
f
(
x
)
d
x
)
⟹
∂
μ
2
∂
M
=
M
f
(
M
)
P
(
X
+
)
-
f
(
M
)
E
(
X
+
)
P
2
(
X
+
)
.After simplification we get(19)
∂
μ
1
∂
M
=
f
(
M
)
P
2
(
X
-
)
[
M
P
(
X
-
)
-
E
(
X
-
)
]
,
∂
μ
2
∂
M
=
f
(
M
)
P
2
(
X
+
)
[
M
P
(
X
+
)
-
E
(
X
+
)
]
.In order to find the minimum ofVar_2, we will find ∂
Var
⁡
1
/
∂
M and ∂
V
a
r
2
/
∂
M separately and add them as follows:
(20)
∂
∂
M
(
Var_2
)
=
∂
∂
M
(
Var
⁡
1
)
+
∂
∂
M
(
Var
⁡
2
)
∂
Var
⁡
1
∂
M
=
∂
∂
M
(
∫
-
∞
M
‍
x
2
f
(
x
)
d
x
+
μ
1
2
∫
-
∞
M
‍
f
(
x
)
d
x
-
2
μ
1
∫
-
∞
M
‍
x
f
(
x
)
d
x
)
=
M
2
f
(
M
)
+
2
μ
1
μ
1
′
P
(
X
-
)
+
μ
1
2
f
(
M
)
-
2
μ
1
′
E
(
X
-
)
-
2
μ
1
M
f
(
M
)
.After simplification we get(21)
∂
Var
⁡
1
∂
M
=
f
(
M
)
[
M
2
+
E
2
(
X
-
)
P
(
X
-
)
2
-
2
M
E
(
X
-
)
P
(
X
-
)
]
.Similarly(22)
∂
Var
⁡
2
∂
M
=
-
M
2
f
(
M
)
+
2
μ
2
μ
2
′
P
(
X
+
)
-
μ
2
2
f
(
M
)
-
2
μ
2
′
E
(
X
+
)
+
2
μ
2
M
f
(
M
)
.After simplification we get(23)
∂
Var
⁡
2
∂
M
=
f
(
M
)
[
-
M
2
+
E
2
(
X
+
)
P
(
X
+
)
2
-
2
M
E
(
X
+
)
P
(
X
+
)
]
.Adding both these differentials and after simplifying, it follows that(24)
∂
∂
M
(
Var_2
)
=
f
(
M
)
[
(
E
2
(
X
-
)
P
(
X
-
)
2
-
E
2
(
X
+
)
P
(
X
+
)
2
)
-
2
M
(
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
)
]
.Equating this with 0 gives(25)
∂
∂
M
(
Var_2
)
=
0
,
[
(
E
2
(
X
-
)
P
(
X
-
)
2
-
E
2
(
X
+
)
P
(
X
+
)
2
)
-
2
M
(
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
)
]
=
0
,
[
E
(
X
-
)
P
(
X
-
)
-
E
(
X
+
)
P
(
X
+
)
]
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
-
2
M
]
=
0
.But[
(
E
(
X
-
)
/
P
(
X
-
)
)
-
(
E
(
X
+
)
/
P
(
X
+
)
)
]
≠
0, as in that case μ
1
=
μ
2, which contradicts the initial assumption. Therefore
(26)
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
-
2
M
]
=
0
⟹
M
=
[
E
(
X
-
)
P
(
X
-
)
+
E
(
X
+
)
P
(
X
+
)
]
·
1
2
,
M
=
μ
1
+
μ
2
2
.Therefore, we have to find all values ofM that satisfy (26), which is easier than minimizing Var_2 directly.To find the minimum in terms ofM, it is enough to derive
(27)
∂
Var_2
(
μ
1
(
M
)
,
μ
2
(
M
)
)
∂
M
=
f
(
M
)
[
E
2
(
X
-
)
P
(
X
<
M
)
2
-
E
2
(
X
+
)
P
(
X
≥
M
)
2
-
2
M
[
E
(
X
-
)
P
(
X
≤
M
)
-
E
(
X
+
)
P
(
X
>
M
)
]
]
.Then we find(28)
M
=
[
E
(
X
+
)
P
(
X
+
)
+
E
(
X
-
)
P
(
X
-
)
]
·
1
2
.


## 2.2.3. Finding the Centers of Some Common Probability Density Functions

(
1
)
Uniform Distribution. The probability density function of the continuous uniform distribution is
(29)
f
(
x
)
=
{
1
b
-
a
,
a
<
x
<
b
0
otherwise
,
where a and b are the two parameters of the distribution. The expected values for X
·
1
x
≥
M and X
·
1
x
<
M are
(30)
E
(
X
+
)
=
1
2
(
b
2
-
M
2
b
-
a
)
,
E
(
X
-
)
=
1
2
(
M
2
-
a
2
b
-
a
)
.Putting all these in (28) and solving, we get
(31)
M
=
1
2
(
a
+
b
)
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
1
4
(
3
a
+
b
)
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
1
4
(
a
+
3
b
)
.(
2
)
Log-Normal Distribution. The probability density function of a log-normal distribution is
(32)
f
(
x
;
μ
,
σ
)
=
1
x
σ
2
π
e
-
(
ln
⁡
x
-
μ
)
2
/
2
σ
2
,
x
>
0
.The probabilities in left and right ofM are, respectively,
(33)
P
(
X
+
)
=
1
2
-
1
2
erf
⁡
(
1
2
2
(
ln
⁡
M
-
μ
)
σ
)
,
P
(
X
-
)
=
1
2
+
1
2
erf
⁡
(
1
2
2
(
ln
⁡
M
-
μ
)
σ
)
.The expected values forX
·
1
x
≥
M and X
·
1
x
<
M are
(34)
E
(
X
+
)
=
[
1
2
e
μ
+
(
1
/
2
)
σ
2
{
1
+
erf
⁡
(
1
2
(
μ
+
σ
2
-
ln
⁡
M
)
2
σ
)
}
]
,
E
(
X
-
)
=
[
1
2
e
μ
+
(
1
/
2
)
σ
2
{
1
+
erf
⁡
(
1
2
(
-
μ
-
σ
2
+
ln
⁡
M
)
2
σ
)
}
]
.Putting all these in (28) and assuming μ
=
0 and σ
=
1 and solving, we get
(35)
M
=
4.244942583
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
1.196827816
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
7.293057348
.(
3
)
Normal Distribution. The probability density function of a normal distribution is
(36)
f
(
x
,
μ
,
σ
)
=
1
σ
2
π
e
-
(
x
-
μ
)
2
/
2
σ
2
.The probabilities in left and right ofM are, respectively,
(37)
P
(
X
+
)
=
1
2
-
1
2
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
,
P
(
X
-
)
=
1
2
+
1
2
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
.The expected values forX
·
1
x
≥
M and X
·
1
x
<
M are
(38)
E
(
X
+
)
=
-
1
2
μ
{
-
1
+
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
}
,
E
(
X
-
)
=
1
2
μ
{
1
+
erf
⁡
(
1
2
2
(
M
-
μ
)
σ
)
}
.Putting all these in (28) and solving we get M
=
μ. Assuming μ
=
0 and σ
=
1 and solving, we get
(39)
M
=
0
,
μ
1
(
M
)
=
E
(
X
-
)
P
(
X
-
)
=
-
2
π
,
μ
2
(
M
)
=
E
(
X
+
)
P
(
X
+
)
=
2
π
.


## 2.2.4. Discrete Case

LetX be a discrete random variable and assume we have a set of N observations from X. Let p
j be the mass density function for an observation x
j.Letx
i
-
1
<
M
i
-
1
<
x
i and x
i
<
M
i
<
x
i
+
1. We define μ
x
≤
x
i as the mean for all x
≤
x
i, and μ
x
≥
x
i as the mean for all x
≥
x
i. Define Var_2
(
M
i
-
1
) as the variance of two centers forced by M
i
-
1 as a separator
(40)
Var
⁡
(
M
i
-
1
)
=
∑
x
j
≤
x
i
-
1
‍
(
x
j
-
μ
x
j
≤
x
i
-
1
)
2
p
j
+
∑
x
j
≥
x
i
‍
(
x
j
-
μ
x
j
≥
x
i
)
2
p
j
,
(41)
Var
⁡
(
M
i
)
=
∑
x
j
≤
x
i
‍
(
x
j
-
μ
x
j
≤
x
i
)
2
p
j
+
∑
x
j
≥
x
i
+
1
‍
(
x
j
-
μ
x
j
≥
x
i
+
1
)
2
p
j
.The first part ofVar
⁡
(
M
i
-
1
) simplifies as follows:
(42)
∑
x
j
≤
x
i
-
1
‍
(
x
j
-
μ
x
j
≤
x
i
-
1
)
2
p
j
=
∑
x
j
≤
x
i
-
1
‍
x
j
2
p
j
-
2
μ
x
j
≤
x
i
-
1
×
∑
x
j
≤
x
i
-
1
‍
x
j
p
j
+
μ
x
j
≤
x
i
-
1
2
∑
x
j
≤
x
i
-
1
‍
p
j
=
E
(
X
2
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
.Similarly,(43)
∑
x
j
≥
x
i
‍
(
x
j
-
μ
x
j
≥
x
i
)
2
p
j
=
E
(
X
2
·
1
x
j
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
,
Var_2
(
M
i
-
1
)
=
E
(
X
2
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
2
·
1
x
j
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
,
Var_2
(
M
i
)
=
E
(
X
2
·
1
x
j
≤
x
i
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
+
E
(
X
2
·
1
x
j
≥
x
i
+
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
.
Define Δ
Var_2
=
Var_2
(
M
i
)
-
Var_2
(
M
i
-
1
)
(44)
Δ
Var_2
=
[
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
[
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
]
.Now, in order to simplify the calculation, we rearrangeΔ
Var_2 as follows:
(45)
Δ
Var_2
=
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
]
+
[
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
.The first part is simplified as follows:(46)
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
)
P
(
X
≤
x
i
)
=
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
+
x
i
p
i
]
2
P
(
X
≤
x
i
-
1
)
+
p
i
=
1
P
(
X
≤
x
i
-
1
)
[
P
(
X
≤
x
i
-
1
)
+
p
i
]
×
{
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
)
+
p
i
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
2
x
i
p
i
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
x
i
2
p
i
2
P
(
X
≤
x
i
-
1
)
}
=
1
P
(
X
≤
x
i
-
1
)
[
P
(
X
≤
x
i
-
1
)
+
p
i
]
×
{
p
i
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
-
2
x
i
p
i
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
x
i
2
p
i
2
P
(
X
≤
x
i
-
1
)
}
.Similarly, simplifying the second part yields(47)
E
2
(
X
·
1
x
j
≥
x
i
)
P
(
X
≥
x
i
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
=
[
E
(
X
·
1
x
j
≥
x
i
+
1
)
+
x
i
p
i
]
2
P
(
X
≥
x
i
+
1
)
+
p
i
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
=
1
P
(
X
≥
x
i
+
1
)
[
P
(
X
≥
x
i
+
1
)
+
p
i
]
×
{
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
-
p
i
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
2
x
i
p
i
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
x
i
2
p
i
2
P
(
X
≥
x
i
+
1
)
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
}
=
1
P
(
X
≥
x
i
+
1
)
[
P
(
X
≥
x
i
+
1
)
+
p
i
]
×
{
-
p
i
+
1
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
+
2
x
i
p
i
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
+
x
i
2
p
i
+
1
2
P
(
X
≥
x
i
+
1
)
-
p
i
+
1
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
}
.In general, these types of optimization arise when the dataset is very big.If we considerN to be very large then p
k
≪
max
⁡
(
∑
i
≥
k
+
1
‍
P
i
;
∑
i
≤
k
-
1
‍
P
i
), P
(
X
≤
x
i
-
1
)
+
p
i
≈
P
(
X
≤
x
i
), and P
(
X
≥
x
i
+
1
)
+
p
i
≈
P
(
X
≥
x
i
+
1
):
(48)
Δ
Var_2
≈
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
p
i
P
2
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
p
i
P
2
(
X
≥
x
i
+
1
)
]
-
2
p
i
x
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
p
i
2
x
i
2
[
1
P
(
X
≤
x
i
-
1
)
-
1
P
(
X
≥
x
i
+
1
)
]
and, as lim
⁡
N
→
∞
{
p
i
2
x
i
2
[
(
1
/
P
(
X
≤
x
i
-
1
)
)
-
(
1
/
P
(
X
≥
x
i
+
1
)
)
]
}
=
0, we can further simplify Δ
Var_2 as follows:
(49)
Δ
Var_2
≈
[
E
2
(
X
·
1
x
j
≤
x
i
-
1
)
p
i
P
2
(
X
≤
x
i
-
1
)
-
E
2
(
X
·
1
x
j
≥
x
i
+
1
)
p
i
P
2
(
X
≥
x
i
+
1
)
]
-
2
p
i
x
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
=
p
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
×
[
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
)
-
2
x
i
]
.To find the minimum ofΔ
Var_2, let us locate when Δ
Var_2
=
0:
(50)
⟹
p
i
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
-
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
×
[
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
)
-
2
x
i
]
=
0
.Butp
i
≠
0, and so is [
(
E
(
X
·
1
x
j
≤
x
i
-
1
)
/
P
(
X
≤
x
i
-
1
)
)
-
(
E
(
X
·
1
x
j
≥
x
i
+
1
)
/
P
(
X
≥
x
i
+
1
)
)
]
<
0, since x
i is an increasing sequence, which implies that μ
i
-
1
<
μ
i where
(51)
μ
i
-
1
=
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
,
μ
i
=
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
.Therefore,(52)
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
-
2
x
i
=
0
.Let us define(53)
F
(
x
i
)
=
2
x
i
-
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
F
(
x
i
)
=
0
⟹
x
i
=
[
E
(
X
·
1
x
j
≤
x
i
-
1
)
P
(
X
≤
x
i
-
1
)
+
E
(
X
·
1
x
j
≥
x
i
+
1
)
P
(
X
≥
x
i
+
1
)
]
·
1
2
⟹
x
i
=
μ
i
-
1
+
μ
i
2
,
where μ
i
-
1 and μ
i are the means of the first and second parts, respectively.To find the minimum for the total variation, it is enough to find the good separatorM between x
i and x
i
-
1 such that F
(
x
i
)
·
F
(
x
i
+
1
)
<
0 and F
(
x
i
)
<
0 (see Figure 3).Figure 3
Discrete random variables.After finding few local minima, we can get the global minimum by finding the smallest among them.


## 3. Face Localization Using Proposed Method

Now, it becomes clear the superiority of the K-means modified algorithm over the standard algorithm in terms of finding a global minimum and we can explain our solution based on the modified algorithm. Suppose we have an image with pixels values from 0 to 255. First, we reshaped the input image into a vector of pixels; Figure4 shows an input image before the operation of the K-means modified algorithm. The image contains a face part, a normal background, and a background with illumination. Then we split the image pixels into two classes: one from 0 to x
1 and another from x
1 to 255 by applying a certain threshold determined by the algorithm.Figure 4
The original image.Let[
0
,
…
,
x
1
] be the class 1-which represents the nonface part and let [
x
1
,
…
,
255
] be the class 2-which represents the face with the shade.All pixels with values less than the threshold belong to the first class and will be assigned to 0 values. This will keep the pixels with values above the threshold and these pixels belong to both the face part and the illuminated part of the background, as shown in Figure5.Figure 5
The face with shade only.In order to remove the illuminated background part, another threshold will be applied to the pixels using the algorithm meant to separate between the face part and the illuminated background part. New classes will be obtained, fromx
1 to x
2 and from x
2 to 255. Then, the pixels with values less than the threshold will be assigned a value of 0 and the nonzero pixels left represent the face part. Let [
x
1
,
…
,
x
2
] be the class 1 which represents the illuminated background part and let [
x
2
,
…
,
255
] be the class 2 which represents the face part.The result of the removal of the illuminated back-ground part is shown in Figure6.Figure 6
The face without illumination.One can see that some face pixels were assigned to 0 due to the effects of the illumination. Therefore, there is a need to return to the original image in Figure4 to fully extract the image window corresponding to the face. A filtering process is then performed to reduce the noise and the result is shown in Figure 7.Figure 7
The filtered image.Figure8 shows the final result which is a window containing the face extracted from the original image.Figure 8
The corresponding window in the original image.


## 4. Results

In this Section, the experimental results of the proposed algorithm will be analyzed and discussed. First, it starts with presenting a description of the datasets used in this work, followed by the results of the previous part of the work. Three popular databases were used to test the proposed face localization method. The MIT-CBCL database has image size115
×
115, and 300 images were selected from this dataset based on different conditions such as light variation, pose, and background variations. The second dataset is the BioID dataset with image size 384
×
286 and 300 images from the database were selected based on different conditions such as indoor/outdoor. The last dataset is the Caltech dataset with image size 896
×
592 and the 300 images from the database were selected based on different conditions such as indoor/outdoor.


### 4.1. MIT-CBCL Dataset

This dataset was established by the Center for Biological and Computational Learning (CBCL) at Massachusetts Institute of Technology (MIT) [41]. It has 10 persons with 200 images per person and an image size of 115
×
115 pixels. The images of each person are with different light conditions, clutter background, and scale and different poses. Few examples of these images are shown in Figure 9.Figure 9
Samples from MIT-CBCL dataset.


### 4.2. BioID Dataset

This dataset [42] consists of 1521 gray level images with a resolution of 384
×
286 pixels. Each one shows the frontal view of a face of one out of 23 different test persons. The images were taken under different lighting conditions, with a complex background, and contain tilted and rotated faces. This dataset is considered as the most difficult one for eye detection (see Figure 10).Figure 10
Sample of faces from BioID dataset for localization.


### 4.3. Caltech Dataset

This database [43] has been collected by Markus Weber at California Institute of Technology. The database contains 450 face images of 27 subjects under different lighting conditions, different face expressions, and complex backgrounds (see Figure 11).Figure 11
Sample of faces from Caltech dataset for localization.


### 4.4. Proposed Method Result

The proposed method is evaluated on these datasets where the result was100
% with a localization time of 4 s. Table 1 shows the proposed method results with the MIT-CBCL, Caltech, and BioID databases. In this test, we have focused on the use of images of faces from different angles from 30 degrees left to 30 degrees right and the variation in the background as well as the cases of indoor and outdoor images. Two sets of images are selected from the MIT-CBCL dataset and two other sets from the other databases to evaluate the proposed method. The first set from MIT-CBCL contains the images with different poses and the second set contains images with different backgrounds. Each set contains 150 images. The BioID and Caltech sets contain images with indoor and outdoor settings. The results showed the efficiency of the proposed K-mean modified algorithm to locate the faces from the image with high background and poses changes. In addition, another parameter was considered in this test which is the localization time. From the results, the proposed method can locate the face in a very short time because it has less complexity due to the use of differential equations.Table 1
Face localization usingK-mean modified algorithm on MIT-CBCL, BioID, and Caltech; the sets contain faces with different poses and faces with clutter background as well as indoors/outdoors.
Dataset
Image
Accuracy (%)
MIT-Set 1
150
100
MIT-Set 2
150
100
BioID
300
100
Caltech
300
100Figure12 shows some examples of face localization on MIT-CBCL.Figure 12
Examples of face localization using k-mean modified algorithm on MIT-CBCL dataset.In Figure13(a), an original image from BioID dataset is shown. The localization position is shown in Figure 13(b). But there is need to do a filtering operation in order to remove the unwanted parts and this is shown in Figure 13(c) and then the area of ROI will be determined.Example of face localization using k-mean modified algorithm on BioID dataset: (a) the original image, (b) face position, and (c) filtered image.
(a)
(b)
(c)Figure14 shows the face position on an image from the Caltech dataset.Example of face localization using k-mean modified algorithm on Caltech dataset: (a) the original image, (b) face position, and (c) filtered image.
(a)
(b)
(c)


## 4.1. MIT-CBCL Dataset

This dataset was established by the Center for Biological and Computational Learning (CBCL) at Massachusetts Institute of Technology (MIT) [41]. It has 10 persons with 200 images per person and an image size of 115
×
115 pixels. The images of each person are with different light conditions, clutter background, and scale and different poses. Few examples of these images are shown in Figure 9.Figure 9
Samples from MIT-CBCL dataset.


## 4.2. BioID Dataset

This dataset [42] consists of 1521 gray level images with a resolution of 384
×
286 pixels. Each one shows the frontal view of a face of one out of 23 different test persons. The images were taken under different lighting conditions, with a complex background, and contain tilted and rotated faces. This dataset is considered as the most difficult one for eye detection (see Figure 10).Figure 10
Sample of faces from BioID dataset for localization.


## 4.3. Caltech Dataset

This database [43] has been collected by Markus Weber at California Institute of Technology. The database contains 450 face images of 27 subjects under different lighting conditions, different face expressions, and complex backgrounds (see Figure 11).Figure 11
Sample of faces from Caltech dataset for localization.


## 4.4. Proposed Method Result

The proposed method is evaluated on these datasets where the result was100
% with a localization time of 4 s. Table 1 shows the proposed method results with the MIT-CBCL, Caltech, and BioID databases. In this test, we have focused on the use of images of faces from different angles from 30 degrees left to 30 degrees right and the variation in the background as well as the cases of indoor and outdoor images. Two sets of images are selected from the MIT-CBCL dataset and two other sets from the other databases to evaluate the proposed method. The first set from MIT-CBCL contains the images with different poses and the second set contains images with different backgrounds. Each set contains 150 images. The BioID and Caltech sets contain images with indoor and outdoor settings. The results showed the efficiency of the proposed K-mean modified algorithm to locate the faces from the image with high background and poses changes. In addition, another parameter was considered in this test which is the localization time. From the results, the proposed method can locate the face in a very short time because it has less complexity due to the use of differential equations.Table 1
Face localization usingK-mean modified algorithm on MIT-CBCL, BioID, and Caltech; the sets contain faces with different poses and faces with clutter background as well as indoors/outdoors.
Dataset
Image
Accuracy (%)
MIT-Set 1
150
100
MIT-Set 2
150
100
BioID
300
100
Caltech
300
100Figure12 shows some examples of face localization on MIT-CBCL.Figure 12
Examples of face localization using k-mean modified algorithm on MIT-CBCL dataset.In Figure13(a), an original image from BioID dataset is shown. The localization position is shown in Figure 13(b). But there is need to do a filtering operation in order to remove the unwanted parts and this is shown in Figure 13(c) and then the area of ROI will be determined.Example of face localization using k-mean modified algorithm on BioID dataset: (a) the original image, (b) face position, and (c) filtered image.
(a)
(b)
(c)Figure14 shows the face position on an image from the Caltech dataset.Example of face localization using k-mean modified algorithm on Caltech dataset: (a) the original image, (b) face position, and (c) filtered image.
(a)
(b)
(c)


## 5. Conclusion

In this paper, we focus on developing the RIO detection by suggesting a new method for face localization. In this method of localization, the input image is reshaped into a vector and then the optimized k-means algorithm is applied twice to cluster the image pixels into classes. At the end, the block window corresponding to the class of pixels containing the face only is extracted from the input image. Three sets of images taken from MIT, BioID, and Caltech datasets and differing in illumination, and backgrounds, as well as in-door/outdoor settings were used to evaluate the proposed method. The results show that the proposed method achieved significant localization accuracy. Our future research direction includes finding the optimal centers for higher dimension data, which we believe is just generalization of our approach to higher dimension (8) and to higher subclusters.



---

*Source: 506480-2014-09-07.xml*