A Batch Rival Penalized Expectation-Maximization Algorithm for Gaussian Mixture Clustering with Automatic Model Selection

Within the learning framework of maximum weighted likelihood (MWL) proposed by Cheung, 2004 and 2005, this paper will develop a batch Rival Penalized Expectation-Maximization (RPEM) algorithm for density mixture clustering provided that all observations are available before the learning process. Compared to the adaptive RPEM algorithm in Cheung, 2004 and 2005, this batch RPEM need not assign the learning rate analogous to the Expectation-Maximization (EM) algorithm (Dempster et al., 1977), but still preserves the capability of automatic model selection. Further, the convergence speed of this batch RPEM is faster than the EM and the adaptive RPEM in general. The experiments show the superior performance of the proposed algorithm on the synthetic data and color image segmentation.


Introduction
As a typical statistical technique, clustering analysis has been widely applied to a variety of scientific areas such as data mining [1], vector quantization [2,3], image processing [4][5][6][7], and so forth. In particular, clustering analysis provides a useful tool to solve the several computer vision problems, for example, multithresholding of gray level images, analysis of the Hough space, and range image segmentation, formulated in the feature space paradigm [8]. In general, one kind of clustering analysis can be formulated as a density mixture modeling, in which each mixture component represents the density distribution of a data cluster. Subsequently, the task of clustering analysis is to identify the dense regions of the input (also called observation interchangeably) densities in a mixture. Such a clustering is therefore called a density mixture clustering.
In general, the Expectation-Maximum (EM) algorithm [9,10] has provided a general solution for the parameter estimation of a density mixture model. Nevertheless, it needs to preassign an appropriate number of density components, that is, the number of clusters. Roughly, the mixture may overfit the data if too many components are utilized, whereas a mixture with too few components may not be flexible enough to approximate the true underlying model. Subsequently, the EM almost always leads to a poor estimate result if the number of components is misspecified. Unfortunately, from the practical viewpoint, it is hard or even impossible to know the exact cluster number in advance. In the literature, one promising way is to develop a clustering algorithm that is able to perform a correct clustering without preassigning the exact number of clusters. Such algorithms include the RPCL algorithm [11] and its improved version, namely, RPCCL [12]. More recently, Cheung [13,14] has proposed a general learning framework, namely, Maximum Weighted Likelihood (MWL), through which an adaptive Rival Penalized EM (RPEM) algorithm has been proposed for density mixture clustering. The RPEM learns the density parameters by making mixture component compete each other at each time step. Not only are the associated parameters of the winning density component updated to adapt to an input, but also all rivals' parameters are penalized with the strength proportional to the corresponding posterior density probabilities. Therefore, this intrinsic rival penalization mechanism enables the RPEM to automatically select an appropriate number of densities by gradually fading out the redundant densities from a density mixture. Furthermore, a simplified version of RPEM has included RPCL and RPCCL as its special cases with some new extensions.
In the papers [13,14], the RPEM algorithm learns the parameters via a stochastic gradient ascending method; that is, we update the parameters immediately and adaptively once the current observation is available. In general, the adaptiveness of the RPEM makes it more applicable to the environment changed over time. Nevertheless, the convergence speed of the RPEM relies on the value of learning rate. Often, by a rule of thumb, we arbitrarily set the learning rate at a small positive constant. If the value of learning rate is assigned too small, the algorithm will converge at a very slow speed. On the contrary, if it is too large, the algorithm may even oscillate. In general, it is a nontrivial task to assign an appropriate value to the learning rate, although we can pay extra efforts to make the learning rate dynamically change over time, for example, see [15].
In this paper, we further study the MWL learning framework and develop a batch RPEM algorithm accordingly provided that all observations are available before the learning process. Compared to the adaptive RPEM, this batch one need not assign the learning rate analogous to the EM, but still preserves the capability of automatic model selection. Further, the convergence speed of this batch RPEM is faster than the EM and the adaptive RPEM in general. The experiments have shown the superior performance of the proposed algorithm on the synthetic data and color image segmentation.
The remainder of this paper is organized as follows. Section 2 reviews the MWL learning framework. In Section 3, we present the batch RPEM algorithm in detail, in which the weights involve a coefficient ε. We will therefore further explore the assignment of ε in Section 4. Section 5 shows the detailed experiment results. Finally, we draw a conclusion in Section 6.

Overview of Maximum Weighted Likelihood (MWL) Learning Framework
Suppose that an input x ∈ d comes from the following density mixture model: where Θ is the parameter set of {α j , θ j } k j=1 . Furthermore, k is the number of components, α j is the mixture proportion of the jth component, and p(x | θ j ) is a multivariate probability density function (pdf) of x parameterized by θ j . As long as we know the value of Θ, an input x can be classified into a certain cluster via its posterior probability: using the winner-take-all rule, that is, assigning an input x to Cluster c if c = arg max j h( j | x, Θ) or taking its soft version which assigns x to Cluster j with the probability h( j | x, Θ). Therefore, how to estimate the parameter set Θ, particularly without knowing the correct value of k in advance, is a key issue in density mixture clustering.
In the MWL learning framework [13,14], the parameter set Θ is learned via maximizing the following Weighted Likelihood (WL) cost function: with where g( j | x, Θ)'s are the designable weights satisfying the two conditions: Suppose that a set of N i.i.d. observations, denoted as X = {x 1 , x 2 , . . . , x N }, comes from the density mixture model in (1). The empirical WL function of (3), written as Q(Θ; X), can be given as with Moreover, the weights g( j | x t , Θ)'s have been generally designed as [13,14] g j | x t , where ε t is a coefficient varying with the time step t in general. Please note that g( j | x t , Θ)'s in (7) can be negative as well as positive. For simplicity, we hereinafter set ε t as a constant, denoted as ε. Moreover, I( j | x t , Θ) is an indicator function with Subsequently, under a specific weight design, the papers [13,14] have presented the adaptive RPEM to learn Θ via maximizing the WL function of (5) using a stochastic gradient ascent method, which is able to fade out the redundant densities gradually from a density mixture. Consequently, it can automatically select an appropriate number of density components in density mixture clustering.
Computational and Mathematical Methods in Medicine 3 Interested readers may refer to the paper [14] for more details. We summarize the main steps of the adaptive RPEM in Algorithm 1. Although the experiments have shown the superior performance of the adaptive RPEM on automatic model selection, its convergence speed, however, relies on the value of learning rate. Under the circumstances, we will present a batch version without the learning rate in the next section.

Batch RPEM Algorithm
To estimate the parameter set within the MWL framework, we have to maximize the empirical WL function Q(Θ; X) in (5). In general, we update the parameters via maximizing the first term of (5), that is, ω(Θ; X), by fixing the second term ν(Θ; X). Subsequently, we need to solve the following nonlinear optimization problem: subject to the constraints as shown in (1). We adopt the Lagrange method analogous to the EM by introducing a Lagrange multiplier λ into the Lagrange function. Subsequently, we have In this paper, we concentrate on the Gaussian mixture model only, that is, each component p(x | θ j ) in (1) is a Gaussian density. We then have where θ j = (μ j , C j ), μ j and C j are the mean (also called seed points interchangeably) and the covariance of the jth density, respectively. Through optimizing (10), we then finally obtain the batch RPEM algorithm as shown in Algorithm 2. If a covariance matrix C (n+1) j is singular, then it indicates that the corresponding jth density component is degenerated and can be simply discarded without being learned any more in the subsequent iterations. In this case, we have to normalize those remaining α (n+1) r 's (r / = j) so that their sum is always kept to be 1.
In the above batch RPEM, its capability of automatic model selection is controlled by the weight functions g( j | x t , θ)'s, which further rely on the parameter ε as shown in (7). Subsequently, a new question is arisen: how to assign an appropriate value of ε? The next section will answer this question.

How to Assign Parameter ε?
To deal with how to assign an appropriate value of ε, we rewrite (7) as the following form: can be regarded as the award received by the winning density component (i.e., the cth density with I(c | x t , Θ) = 1), and meanwhile is the penalty of the rival components (i.e., those densities with I( j | x t , Θ) = 0). In general, it is expected that the award is positive and the penalty is negative. That is, ε should be greater than −1.
Otherwise, as ε < −1, we will meet an awkward situation; that is, the amount of award is negative and the penalty one becomes positive. This implies that we will penalize the winner and award the rivals, which evidently violates our expectations. Furthermore, as ε = −1, both of the award and penalty amount become zero. In this special case, the batch RPEM is actually degenerated into the EM without the property of automatic model selection. As a result, ε is required to be greater than −1. In addition, ε in the batch RPEM should be a negative value. Otherwise, the weights of the rival components g( j | x t , Θ) = −εh( j | x t )'s become negative, resulting in some α j 's to be negative finally. Hence, an appropriate selection of ε in the batch RPEM would be a negative value and greater than −1. That is, ε should be fallen into the range of (−1, 0). Furthermore, our empirical studies have found that a smaller ε will lead the batch RPEM algorithm to a more robust performance, especially when the value of k is large and the data are overlapped considerably. In other words, the algorithm has a poor capability of automatic model selection if ε is close to zero. To illustrate this scenario, we have utilized two synthetic data sets that are generated from the two bivariate three-Gaussian mixtures individually as shown in Figures 1(a) and 1(b), where each data set consists of 1, 000 observations with the true mixture proportions: α * 1 = 0.4, α * 2 = 0.3, and α * 3 = 0.3. Also, the true μ * j 's and C * j 's of data set 1 in Figure 1 4

Computational and Mathematical Methods in Medicine
Initialization: Given a specific k (k ≥ k * , k * is the true number of clusters), initialize the parameter Θ.
Step 2: Given h( j | x t , Θ (n) )'s and g( j | x t , Θ (n) )'s, we update Θ by , where η is a small positive learning rate.
Step 3: Let n = n + 1, and go to Step 1 for the next iteration until Θ is converged.
Initialization: Given a specific k (k ≥ k * , k * is the true number of clusters), initialize the parameter Θ.
Step 2: Fixing h( j | x t , Θ (n) )'s and g( j | x t , Θ (n) )'s, we update Θ by Step 3: Let n = n + 1, and go to Step 1 for the next iteration until Θ is converged.
while the true parameters of data set 2 in Figure 1 It can be seen that the clusters in data set 1 are well separated, whereas the clusters in data set 2 are overlapped considerably.
For each data set, we conducted the three experiments by setting k = 3, k = 8, and k = 20, respectively. Also, all α j 's and C j 's were initialized at 1/k and the identity matrix, respectively. During the learning process, we discarded those densities whose covariance matrices C j 's were singular. Table 1 shows the performance of the batch RPEM over the parameter ε. We found that, as k = 3 and k = 8, all ε's we have tried from −0.9 to −0.1 lead to the good performance of the algorithm when using the data set 1. For example, as k = 8 and ε = −0.8, we randomly initialized the eight seed points in the input space as shown in Figure 2(a). After all the parameters were converged, 2 out of 8 density components had been discarded and the mixture proportions of the remaining components were converged to α 1 = 0.2960, α 2 = 0.0036, α 3 = 0.2900, α 4 = 0.0058, α 5 = 0.0136, and α 6 = 0.3910. It can be seen that the three principal mixing proportions, α 1 , α 3 , and α 6 , have well estimated the true ones, while the other proportions were inclined to zero. The corresponding three μ j 's and C j 's were μ 1 As shown in Figure 2(b), the three μ j 's have successfully stabilized at the corresponding cluster centers, meanwhile the other three redundant seed points have been pushed away and stably located at the boundary of the clusters. That is, the redundant densities have been fade out through the learning, thus the batch RPEM can select the model automatically as well as the adaptive version. Nevertheless, when k is set at a large value, for example, say k = 20, it is found that the proposed algorithm could not fade out the redundant density components from a mixture if ε is close to 0. Instead, we should set ε at a value close to −1. For example, as k * = 3, k = 20, and ε = −0.9, we ran the proposed algorithm. It was found that  In other words, the algorithm has a poor performance as ε get close to zero. Hence, if k is large, it would be better to choose a relative smaller value of ε between −1 and 0.
In addition, we also investigated the assignment of ε on data set 2, where the data are considerably overlapped. We take k * = 3, k = 20, and ε = −0.9 for instance. The converged positions of the seed points are shown in Figure 4(a), where the learned positions converged to the cluster centers while driving the redundant seed points to the boundaries of the clusters. That is, the proposed batch algorithm can work quite well as ε = −0.9. Also, we let k * = 3, k = 20, and ε = −0.2 to run the algorithm again for comparison. As shown in Figure 4(b), the converged positions of the seed points have a bias from the cluster centers. This implies that the values of ε that are close to zero cannot work well in this case. More examples can be found in Table 1. It can be seen that the feasible region of ε is shrunk over the overlap level of the data. For example, the appropriate values of ε are in the range of [−0.9, −0.6] only when using the date set 2 with k * = 3 and k = 3 or 8, respectively. In contrast, ε is feasible in the full range of [−0.9, −0.1] where we have tried so far as data set 1 is used. Hence, by a rule of thumb, we should choose an appropriate value of ε close to −1. Nevertheless, we have also noted that it is not a good choice if ε is too close to −1. In fact, the proposed algorithm will gradually degenerate to the EM as ε tends to −1; that is, the capability of the proposed algorithm on model selection will be reduced as ε tends to −1. Hence, to sum up, empirical studies have found that [−0.9, −0.8] is an appropriate feasible region of ε. In the next section, we therefore arbitrarily set ε at −0.8.

Experimental Results
To evaluate the performance of the batch RPEM algorithm, we have conducted the following three experiments.

Experiment 1: Batch RPEM on Synthetic Data with
This experiment was to evaluate the convergence speed of the batch RPEM. We utilized 1, 000 data points from a mixture of three bivariate Gaussian densities with the true parameters as follows:  where "G" stands for a good model selection capability of the algorithm, while "P" represents a poor model selection capability.
We let k = 3, which is equal to the true mixture number k * = 3. The three seed points were randomly allocated in the observation space as shown in Figure 5(a), where the data are considerably overlapped. Moreover, all α j 's and C j 's were initialized at 1/k and the identity matrix, respectively. Figure 5(d) shows the positions of the three converged seed points, which are all stably located at the corresponding cluster centers. For comparison, we also implemented the EM under the same experimental environment. Figure 5(b) shows that the EM had successfully located the three seed points as well as the batch RPEM. Nevertheless, as shown in Figures 6(c) and 7, the batch RPEM converges at 20 epochs, while the EM needs 60 epochs as shown in Figure 6(a). That is, the convergence speed of the batch RPEM is significantly faster than the EM. This indicates that the intrinsic rival-penalization scheme of the batch RPEM, analogous to the RPCL [11], RPCCL [12], and the adaptive RPEM [14], is able to drive away the rival seed points so that they can be more quickly towards the other cluster centers. As a result, the batch RPEM converges much faster than the EM. Moreover, we also compared it with the adaptive RPEM, in which we set the learning rate η = 0.001. Figure 5(c) shows the convergent results of the seed points. It can be seen that the adaptive RPEM works quite well in this case, but it needs around 70 epochs as shown in Figure 6     appropriate learning rate is adopted, which, however, is not a trivial task.

Experiment 2:
Batch RPEM on Synthetic Data with K > K * . This experiment will investigate the performance of batch RPEM performance as k > k * . We generated 1, 000 observations from a mixture of five bivariate Gaussian density distributions with the mixing proportions:  the red-green-blue (RGB) color space model that represents each pixel in an image by a three-color vector. We conducted color image segmentation on a 122 × 152 hand image and a 96 × 128 house image as shown in Figures 9(a) and 10, respectively. For the former, we initially assigned 10 seed points randomly. After the convergence of the algorithms' performance, a snapshot of their segmentation results is shown in Figures 9(b) and 9(c). It can be seen that the blue tiny swim ring-shaped region after segmentation process by the batch RPEM is smoother than the EM. Further, the tiny nail regions have been partitioned by the batch RPEM but the EM is not. In other words, the batch RPEM algorithm performs better than the EM algorithm. For the house image, we initially assigned the seed points to be 80. A snapshot of the converged segmentation results of the EM and the batch RPEM is shown in Figure 11. It can be seen that the texture on the red wall and the green lawn has no longer maintained after the segmentation process both by the EM and the RPEM. However, the small white regions of windows on red wall were disappeared by the EM as well as the triangle shadow area on the wall. In contrast, the batch RPEM algorithm partitioned these regions well as shown in Figure 11(b). Actually, the batch RPEM has drove out the redundant seed points far away and maintained some principal components correctly, which therefore leads to a better performance in color image segmentation.

Conclusion
In this paper, we have developed a batch RPEM algorithm based on MWL learning framework for Gaussian mixture clustering. Compared to the adaptive RPEM, this new one need not select the value of learning rate. As a result, it can learn faster in general and still preserve the capability of automatic model selection analogous to the adaptive one. We have evaluated the proposed batch RPEM algorithm on both synthetic data and color image segmentation. The numerical results have shown the efficacy of the proposed algorithm.