J J Verbeek
email: jverbeek@science.uva.nl

N Vlassis
email: vlassis@science.uva.nl

B Kröse
email: krose@science.uva.nl

Efficient Greedy Learning of Gaussian Mixture Models

Keywords: unsupervised learning, finite mixture models, Gaussian mixtures, EM algorithm, greedy learning

This paper concerns the greedy learning of Gaussian mixtures. In the greedy approach, mixture components are inserted into the mixture one after the other. We propose a heuristic for searching for the optimal component to insert. In a randomized manner a set of candidate new components is generated. For each of these candidates we find the locally optimal new component. The best local optimum is then inserted into the existing mixture. The resulting algorithm resolves the sensitivity to initialization of state-of-the-art methods, like EM, and has running time linear in the number of data points and quadratic in the (final) number of mixture components. Due to its greedy nature the algorithm can be particularly useful when the optimal number of mixture components is unknown. Experimental results comparing the proposed algorithm to other methods on density estimation and texture segmentation are provided.

Introduction

This paper concerns the learning (fitting the parameters of) a mixture of Gaussian distributions [START_REF] Mclachlan | Finite Mixture Models[END_REF]. Mixture models form an expressive class of models for density estimation. Applications in a wide range of fields have emerged in the past decades. They are used for density estimation in 'unsupervised' problems, for clustering purposes, for estimating class-conditional densities in supervised learning settings, in situations where the data is partially supervised and in situations where some observations have 'missing values'.

The most popular algorithm to learn mixture models is the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. For a given finite data set X n of n observations and an initial mixture f 0 , the algorithm provides a means to generate a sequence of mixture models {f i } with non-decreasing log-likelihood on X n . The EM algorithm is known to converge to a locally optimal solution. However, convergence to a globally optimal solution is not guaranteed. The log-likelihood of the given data set under the found mixture distribution is highly dependent on the initial mixture f 0 .

The standard procedure to overcome the high dependence on initialization is to start the EM algorithm for several random initializations and use the mixture yielding maximum likelihood on the data. In this paper we present a method to learn finite Gaussian mixtures, based on the approximation result of [START_REF] Li | Mixture density estimation[END_REF], which states that we can 'quickly' approximate any target density with a finite mixture model, as compared to the best we can do with any (possibly non-finite) mixture from the same component class. This result also holds if we learn the mixture models in a greedy manner: start with one component and optimally add new components one after the other. However locating the optimal new component boils down to finding the global maximum of a log-likelihood surface. In [START_REF] Li | Estimation of Mixture Models[END_REF] it is proposed to grid the parameter space to locate the global maximum, but this is unfeasible when learning mixtures in high dimensional spaces.

Here, we propose a search heuristic to locate the global maximum. Our search procedure selects the best member from a set of candidate new components. The set of candidate components is 'dynamic' in the sense that it depends on the mixture learned so far. Our search procedure for a new component has running time O(n). For learning a k component mixture this results in a run time of O(nk 2 ) if the complete mixture is updated with EM after each insertion. If the mixture is not updated between the component insertions (the approximation result still holds in this case) the run time would be O(nk).

The paper is organized as follows: In Section 2 we recapitulate the definition and EM-learning of Gaussian mixtures. Section 3 discusses greedy learning of Gaussian mixtures. Our component insertion procedure is discussed in Section 4, which forms the core of the paper. Then, in Section 5, we present experimental results on two tasks: modeling of artificially generated data drawn from several types of Gaussian mixtures and texture segmentation. The experiments compare the performance of the new algorithm with the performance of several other methods to learn Gaussian mixtures. Section 6 ends the paper with conclusions and a discussion.

Gaussian mixtures and the EM algorithm

A Gaussian mixture (GM) is defined as a convex combination of Gaussian densities. A Gaussian density in a d-dimensional space, characterized by its mean m ∈ IR d and d × d covariance matrix C, is defined as:

φ(x; θ) = (2π) -d/2 det(C) -1/2 exp (-(x -m) ⊤ C -1 (x -m)/2), (1)
where θ denotes the parameters m and C. A k component GM is then defined as:

f k (x) = k j=1 π j φ(x; θ j ), with k j=1
π j = 1 and for j ∈ {1, . . . , k} : π j ≥ 0.

(2)

The π j are called the mixing weights and φ(x; θ j ) the components of the mixture. As compared to non-parametric density estimation, mixture density estimation (or: semi-parametric density estimation) allows for efficient evaluation of the density, since only relatively few density functions have to be evaluated. Furthermore, by using few mixture components a smoothness constraint on the resulting density is enforced, allowing for robust density estimates.

The well known Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] enables us to update the parameters of a given k-component mixture with respect to a data set X n = {x 1 , . . . , x n } with all x i ∈ IR d , such that the likelihood of X n is never smaller under the new mixture. The updates for a GM can be accomplished by iterative application of the following equations for all components j ∈ {1, . . . , k}:

P (j | x i ) := π j φ(x i ; θ j )/f k (x i ),
(3)

π j := n i=1 P (j | x i )/n, (4)
m j := n i=1 P (j | x i )x i /(nπ j ), (5)
C j := n i=1 P (j | x i )(x i -m j )(x i -m j ) ⊤ /(nπ j ). ( 6
)
As already mentioned in the previous section, the EM algorithm is not guaranteed to lead us to the best solution, where 'best' means the solution yielding maximal likelihood on X n among all maxima of the likelihood. (We implicitly assume throughout that the likelihood is bounded, by restricting the parameter space, and hence the maximum likelihood estimator is known to exist [START_REF] Lindsay | The geometry of mixture likelihoods: a general theory[END_REF].) The good thing is that if we are 'close' to the global optimum of the parameter space, then it is very likely that by using EM we obtain the globally optimal solution.

In this section we discuss the greedy approach to learning GMs. The basic idea is simple: Instead of starting with a (random) configuration of all components and improve upon this configuration with EM, we build the mixture component-wise. We start with the optimal one-component mixture, whose parameters are trivially computed. Then we start repeating two steps until a stopping criterion is met. The steps are: (i) insert a new component and (ii) apply EM until convergence.

The stopping criterion can implement the choice for a pre-specified number of components or it can be any model complexity selection criterion (like Minimum Description Length, Akaike Information Criterion, Cross Validation, etc.).

Motivation

The intuitive motivation for the greedy approach is that the optimal component insertion problem involves a factor k fewer parameters than the original problem of finding a near optimal start configuration, we therefore expect it to be an easier problem. Of course, the intuitive assumption is that we can find the optimal (with respect to log-likelihood) (k + 1)-component mixture by local search if we start the local search from the mixture obtained by inserting optimally a new component in the optimal k-component mixture. We recently applied a similar strategy to the k-means clustering problem [START_REF] Likas | The global k-means clustering algorithm[END_REF], for which encouraging results were obtained. Two recent theoretical results provide complementary motivation. The Kullback-Leibler (KL) divergence between two densities f and g is defined as:

D(f g) = Ω f (x) log [f (x)/g(x)]dx, ( 7
)
where Ω is the domain of the densities f and g, see [START_REF] Cover | Elements of Information Theory[END_REF] for details. In [START_REF] Li | Mixture density estimation[END_REF] it is shown that for an arbitrary probability density function f there exists a sequence {f i } of finite mixtures such that f k (x) = k i=1 π i φ(x; θ i ) achieves Kullback-Leibler (KL) divergence D(f f k ) ≤ D(f g P ) + c/k for every g P = φ(x; θ)P (dθ). Hence, the difference in KL divergence achievable by k-component mixtures and the KL divergence achievable by any (possibly non-finite) mixture from the same family of components tends to zero with a rate of c/k, where c is a constant not dependent on k but only on the component family. Furthermore, it is shown in [START_REF] Li | Mixture density estimation[END_REF] that this bound is achievable by employing the greedy scheme discussed in Section 3.2. This tells us that we can 'quickly' approximate any density by the greedy procedure. Therefore, we might expect the results of the greedy procedure-as compared to methods that initialize all-at-once-to differ more when fitting mixtures with many components.

The sequence of mixtures generated by the greedy learning method can conveniently be used to guide a model selection process in the case of an unknown number of components. Recently a result of 'almost' concavity of the log-likelihood of a data set under the maximum-likelihood kcomponent mixture, as function of the number of components k was presented [START_REF] Cadez | On model selection and concavity for finite mixture models[END_REF]. The result states that the first order Taylor approximation of the log-likelihood of the maximum likelihood models as function of k is concave under very general conditions. Hence, if we use a penalized log-likelihood model selection criterion based on a penalty term which is concave or linear in k then the penalized log-likelihood is almost concave. This implies that if the 'almost' concave turns out to be concave then there is only one peak in the penalized log-likelihood and hence this peak can be relatively easily identified, e.g. with a greedy construction method.

A General scheme for greedy learning of mixtures

Let L(X n , f k ) = n i=1 log f k (x i ) (we will just write L k if no confusion arises) denote the loglikelihood of the data set X n under the k-component mixture f k . The greedy learning procedure outlined above can be summarized as follows:

1. Compute the optimal (yielding maximal log-likelihood) one-component mixture f 1 . Set k:=1.

2. Find the optimal new component φ(x; θ * ) and corresponding mixing weight α * :

{θ * , α * } = arg max {θ,α} n i=1 log [(1 -α)f k (x i ) + αφ(x i ; θ)] (8)
while keeping f k fixed.

3. Set f k+1 (x) := (1 -α * )f k (x) + α * φ(x; θ * ) and k := k + 1; 4. Update f k using EM until convergence. [optional]
5. If a stopping criterion is met then quit, else go to step 2.

Clearly, the crucial step is the component insertion (step 2), which we will be the topic of the next section. Let us first make a few remarks before we turn to the next section: (i) The stopping criterion in step 5 can be used to force the algorithm to find a mixture of a pre-specified number of components. Of course, step 5 may also implement any kind of model selection criterion. (ii) Note that step 4 may also be implemented by other algorithms than EM. Furthermore, step 4 is not needed to obtain the approximation result of [START_REF] Li | Mixture density estimation[END_REF]. (iii) One may note similarity to the Vertex Direction Method (VDM) [START_REF] Böhning | A review of reliable maximum likelihood algorithms for semiparametric mixture models[END_REF][START_REF] Lindsay | The geometry of mixture likelihoods: a general theory[END_REF][START_REF] Wu | Some algorithmic aspects of the theory of optimal designs[END_REF] to learn Gaussian Mixtures. Indeed, VDM may be regarded as a specific choice for implementing the search needed for step 2, as we explain in the discussion in Section 6.

Efficient search for new components

Suppose we have obtained a k-component mixture f k , the quest is to find the component characterized by equation ( 8). It is easily shown that if we fix f k and θ then L k+1 is concave as function of α only, allowing efficient optimization. For example, in [START_REF] Vlassis | A greedy EM algorithm for Gaussian mixture learning[END_REF] a second order Taylor approximation was used. However, L k+1 as function of θ can have multiple maxima. Hence, we have to perform a global search among the new components in order to identify the optimum since no constructive method to find the global maximum is known.

In our new algorithm, the global search for the optimal new component is achieved by starting km 'partial' EM searches from initial parameter pairs {(θ i , α i )}, 1 ≤ i ≤ km. By 'partial' EM search we mean that we fix f k and optimize L k+1 over θ and α only. The use of partial EM searches is dictated by the general scheme above, for which the results of [START_REF] Li | Mixture density estimation[END_REF] hold. One may wonder why we would use such partial searches, since we might as well optimize over all parameters of the resulting (k + 1)-component mixture. The answer is that if we would update all parameters, then the number of computations needed in each of the km individual searches would be O(nk) instead of O(n).

Each partial search starts with a different initial configuration. After these multiple partial searches we end up with km 'candidate' new components together with their mixing weights. We pick that candidate component φ(x; θ) that maximizes the likelihood when mixed into the previous mixture by a factor α as in (8). Then, in step 3 of the general algorithm, in place of the (unknown) optimal parameter pair (θ * , α * ) we use ( θ, α).

Note that a simple uniform quantization of the parameter space, as implemented for VDM for mixtures of univariate Gaussian distributions in [START_REF] Dersimonian | Maximum likelihood estimation of a mixing distribution[END_REF], is not feasible in general since the number of quantization levels grows exponentially with the dimension of the parameter space. In the following section we discuss the step 2 as proposed in [START_REF] Vlassis | A greedy EM algorithm for Gaussian mixture learning[END_REF] and its drawbacks. In Section 4.2 we present our improved search procedure.

Every data point generates a candidate

In [START_REF] Vlassis | A greedy EM algorithm for Gaussian mixture learning[END_REF] (we will from now on refer to this method as VL) it is proposed to use n candidate components. Every data point is the mean of a corresponding candidate. All candidates have the same covariance matrix σ 2 I, where σ is taken at a theoretically justifiable value. For each candidate component the mixing weight α is set to the mixing weight maximizing the second order Taylor approximation of L k+1 around α = 1/2. The candidate yielding highest log-likelihood when inserted in the existing mixture f k is selected. The selected component is then updated using partial EM steps before it is actually inserted into f k to give f k+1 . Note that for every component insertion that is performed, all point-candidates are considered. There are two main drawbacks to the aforementioned method:

1. Using n candidates in each step results in a time complexity O(n 2 ) for the search which is unacceptable in most applications. O(n 2 ) computations are needed since the likelihood of every data point under every candidate component has to be evaluated.

2. By using fixed small (referring to the size of the determinant of the covariance matrix) candidate components, the method keeps inserting small components in high density areas while in low density areas the density is modeled poorly. We observed this experimentally. Larger components that give greater improvement of the mixture are not among the candidates, nor are they found by the EM procedure following the component insertion.

Working on improvements of the above algorithm, we noticed that both better and faster performance could be achieved by using a smaller set of candidates that are tailored to the existing mixture. In the next section we discuss how we can generate the candidates in a data-and mixture-dependent manner.

Generating few good candidates

Based on our experience with the above method, we propose a new search heuristic for finding the optimal new mixture component. Two observations motivate the new search method:

1. The size (i.e. the determinant of the covariance matrix) of the new component should in general be smaller than the size of the existing components of the mixture.

2. As the existing mixture f k contains more components, the search for the optimal new component should become more thorough.

In our search strategy we account for both observations by (i) setting the size of the initial candidates to a value related to (and in general smaller than) the size of the components in the existing mixture and (ii) increasing the number of candidates linearly with k.

For each insertion problem, our method constructs a fixed number of candidates per existing mixture component. Based on the posterior distributions, we partition the data set X n in k disjoint subsets A i = {x ∈ X n : i = arg max j {P (j | x)}}. If one is using EM for step 4 (of the general scheme given above), then the posteriors P (i | x) are available directly since we have already computed them for the EM updates of the k-component mixture. For each set A i (i = 1, . . . , k) we construct m candidate components. In our experiments we used m = 10 candidates but more can be used, trading run-time for better performance. In fact, like in [START_REF] Smola | Sparse greedy matrix approximation for machine learning[END_REF], we can compute confidence bounds of the type: "With probability at least 1 -δ the best split among m uniformly randomly selected splits is among the best ǫ fraction of all splits, if m ≥ ⌈log δ/ log (1 -ǫ)⌉."

To generate new candidates from A i , we pick uniformly random two data points x l and x r in A i . Then, we partition A i into two disjoint subsets A il and A ir . For elements of A il the point x l is closer than x r and vice versa for A ir . The mean and covariance of the sets A il and A ir are used as parameters for two candidate components. The initial mixing weights for candidates generated from A i are set to π i /2. To obtain the next two candidates we draw new x l and x r , until the desired number of candidates is reached.

The initial candidates can be replaced easily by candidates that yield higher likelihood when mixed into the existing mixture. To obtain these better candidates we perform the partial EM searches, starting from the initial candidates. Each iteration of the km partial updates takes O(mnk) computations, since we have to evaluate the likelihood of each datum under each of the mk candidate components. In the following section we discuss how we can reduce the amount of computations needed by a factor k resulting in only O(mn) computations to perform one iteration of the partial updates.

We stop the partial updates if the change in log-likelihood of the resulting (k + 1)-component mixtures drops below some threshold or if some maximal number of iterations is reached. After these partial updates we set the new component φ(x; θ k+1 ) as the candidate that maximizes the log-likelihood when mixed into the existing mixture with its corresponding mixing weight.

An example is given in Figure 1, which depicts the evolution of a solution for artificially generated data. The mixtures f 1 , . . . , f 4 are depicted, each component is shown as an ellipse which has the eigenvectors of the covariance matrix as axes and radii of twice the square root of the corresponding eigenvalue.

Speeding up the partial EM searches

In order to achieve O(mn) time complexity for the partial EM searches initiated at the km initial candidates, we use a lower bound on the true log-likelihood instead of the true log-likelihood itself. The lower-bound is obtained by assuming that the support (non-zero density points in X n ) for a candidate component φ(x, m k+1 , C k+1 ) originating from existing component i is fully contained in A i . The lower-bound, termed negative free-energy [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF], is detailed in Appendix A. Optimizing the lower-bound allows us to base the partial updates for a candidate from A i purely on the data in A i . The update equations are:

p(k + 1 | x j ) = αφ(x j ; m k+1 , C k+1 ) (1 -α)f k (x j ) + αφ(x j ; m k+1 , C k+1 ) (9) m k+1 = j∈Ai p(k + 1 | x j )x j j∈Ai p(k + 1 | x j ) (10)
C k+1 = j∈Ai p(k + 1 | x j )(x j -m k+1 )(x j -m k+1 ) ⊤ j∈Ai p(k + 1 | x j ) (11) α = j∈Ai p(k + 1 | x j ) N (12)

Time Complexity Analysis

As mentioned in the introduction, the total run time of the algorithm is O(k 2 n) (or O(kmn) if k < m). This is due to the updates of the mixtures f i , which cost O(ni) computations each if we use EM for these updates. Therefore, summing over all mixtures we get O(

k i=1 ni) ∈ O(nk 2
). This is a factor k slower than the standard EM procedure. The run times of our experiments confirm this analysis, the new method performed on average about k/2 times slower than standard EM. However, if we would use EM to learn a sequence of mixtures consisting of 1, . . . , k components then the run time would also be O(nk 2 ).

Note that if we do not use the aforementioned speed-up, the run time would become O(k 2 mn) since in that case the component allocation step for the i + 1-st mixture component would cost O(imn). This is a factor m more than the preceding EM steps, hence the total run-time goes up by a factor m.

Experimental Results

In this section we present results of two experiments. The experiments compare the results as obtained with k-means initialized EM, the method of VL discussed in Section 4.1 , Split and Merge EM (SMEM) [START_REF] Ueda | SMEM algorithm for mixture models[END_REF] and the new greedy approach presented here. The SMEM algorithm checks after convergence of EM if the log-likelihood can be improved by splitting one mixture component in two and merging two other mixture components. To initialize SMEM we also used the k-means initialization. The k-means initialization makes use of the results of applying the k-means or Generalized Lloyd Algorithm to the data [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]. The means are initialized as the centroids of the Voronoi Regions (VR) found, and the covariance matrices as the covariance matrix corresponding to the data in the VRs. The mixing weights were initialized as the proportion of the data present in each VR. The centroids of the k-means algorithm were initialized as a uniformly random selected subset of the data.

Artificial Data

In this experiment we generated data sets of 400 points in IR d with d ∈ {2, 5}. The data was drawn from a Gaussian mixture of k ∈ {4, 6, 8, 10} components. The separation of the components was chosen c ∈ {1, 2, 3, 4}. A separation of c means [START_REF] Dasgupta | Learning mixtures of Gaussians[END_REF]:

∀ i =j : µ i -µ j 2 ≥ c • max {i,j} {trace(C i ), trace(C j )}. ( 13
)
For each mixture configuration we generated 50 data sets. We allowed a maximum eccentricity (i.e. largest singular value of the covariance matrix over the smallest) of 15 for each component. Also, for each generated mixture, we generated a test set of 200 points not presented to the mixture learning algorithm. In the comparison below, we evaluate the difference in log-likelihood of the test sets under the mixtures provided by the different methods.

In Figure 2 we give show experimental results comparing our greedy approach to VL, SMEM and several runs of normal EM initialized with k-means. The average of the difference in loglikelihood on the test set is given for different characteristics of the generating mixture. The results show that the greedy algorithm performs comparable to SMEM and generally outperforms VL. Both greedy methods and SMEM outperform the k-means initialization by far.

Texture segmentation

In this experiment the task is to cluster a set of 16 × 16 = 256 pixel image patches. The patches are extracted from 256 × 256 Brodatz texture images, four of these textures are shown in Figure 3. Since different patches from the same texture differ, roughly speaking, from each other only in limited number of degrees of freedom (translation, rotation, scaling and lightness), patches from the same texture are assumed to form clusters in the 256 dimensional patch-space. Figure 2: Detailed exposition of log-likelihood differences. The upper three tables give the loglikelihood of a method subtracted from the log-likelihood obtained using our greedy approach.

Vlassis & Likas

The bottom table gives averages of the log-likelihood under the generating mixture minus the log-likelihood given by our method. We conducted experiments where the number of textures from which patches were extracted ranged from k = {2, . . . , 6}. For each value of k, we created 100 data sets of 500k patches each by picking up 500 patches at random from each of k randomly selected textures. The experiment comprised learning a Gaussian mixture model with as many components as different textures involved. We compared our greedy approach to the same methods as in the previous experiment.

To speed-up the experiment, we projected the data into a linear subspace with PCA. We used a 50 dimensional subspace, or to a lower than 50 dimensional space if that could capture at least 80% of the data variance. Once the mixture model was learned, we segmented the patches by assigning each patch to the mixture component that has the highest posterior probability for this patch. To evaluate the quality of the different segmentations discovered by the different learning methods, we considered how informative the found segmentations are. We measure how predictable the texture label of a patch is, given the cluster of the class. This can be done using the conditional entropy [START_REF] Cover | Elements of Information Theory[END_REF] of the texture label given the cluster label. In Figure 4 the average conditional entropies (over 50 experiments) are given. The results of the greedy method and SMEM are roughly comparable, although the greedy approach was much faster. Both provide generally more informative segmentations than using the k-means initialization. Note that VL fails to produce informative clusters in this experiment. This is probably due to the high dimensionality of the data that renders spherical candidate components rather unfit.

Discussion and conclusions

Discussion Both VDM [START_REF] Böhning | A review of reliable maximum likelihood algorithms for semiparametric mixture models[END_REF][START_REF] Lindsay | The geometry of mixture likelihoods: a general theory[END_REF][START_REF] Wu | Some algorithmic aspects of the theory of optimal designs[END_REF] and the proposed method are instantiations of the more general greedy scheme given in section 3.2 (although VDM skips step 4 of the scheme). VDM makes use of the directional derivative D f k (φ) of the log-likelihood for the current mixture f k , where

D f k (φ) = lim α→0 [L(X n , (1 -α)f k + αφ) -L(X n , f k )]/α.
VDM proceeds by picking the φ * that maximizes D f k (φ) and inserting this in the mixture with a factor α * such that α * = arg max α {L(X n , (1 -α)f k + αφ * }. Using the directional derivative at the current mixture can be seen as an instantiation of step 2 of the general scheme. The optimization over both φ and α is replaced by (i) an optimization over D f k (φ) (typically implemented by griding the parameter space) and then (ii) an optimization over α. Note that by moving in the direction of maximum D f k (φ) does not guarantee that we move in the direction of maximum improvement of log-likelihood if we optimize over α subsequently.

Recently, several novel methods to learn mixtures (of Gaussian densities) were proposed among which we mention [START_REF] Dasgupta | Learning mixtures of Gaussians[END_REF][START_REF] Figueiredo | Unsupervised learning of finite mixture models[END_REF][START_REF] Ueda | SMEM algorithm for mixture models[END_REF]. The first tries to overcome the difficulties of learning mixtures of Gaussian densities in high dimensional spaces. By projecting the data to a lower dimensional subspace (which is chosen uniformly randomly), finding the means in that space and then projecting them back, the problems of high dimensionality are reduced. The last two methods try to overcome the dependence of EM on the initial configuration as does our method. In [START_REF] Ueda | SMEM algorithm for mixture models[END_REF] split and merge operations are applied to local optima solutions found by applying EM. The split and merge operations constitute jumps in the parameter space that allow the algorithm to jump from a local optimum to a better region in the parameter space. By then reapplying EM a better (hopefully global) optimum is found.

One may note that the latter method and the greedy approach are complementary, i.e. they can be combined. The split-and-merge operations can be incorporated in the greedy algorithm by checking after each component insertion whether a component i can be removed such that the resulting log-likelihood is greater than the likelihood before insertion. If so, the component i is removed and the algorithm continues. If not, the algorithm just proceeds as usual. However, in experimental studies on artificially generated data we found that this combination hardly gave improvement over the individual methods.

An important benefit of our new method over [START_REF] Ueda | SMEM algorithm for mixture models[END_REF] is that the new algorithm produces a sequence of mixtures that can be used to perform model complexity selection as the mixtures are learned. For example a kurtosis-based selection criterion, like the one in [START_REF] Vlassis | A kurtosis-based dynamic approach to Gaussian mixture modeling[END_REF], can be used here. We also note that our algorithm executes faster than SMEM, which has a run time of O(nk 3 ). In [START_REF] Figueiredo | Unsupervised learning of finite mixture models[END_REF] it is proposed to start with a large number of mixture components and to successively annihilate components with small mixing weights. This approach can be characterized as 'pruning' a given mixture, whereas our approach can be characterized as 'growing' a mixture. As compared to both these methods, ours does not need an initialization scheme, since the optimal one-component mixture can be found analytically.

Also Bayesian methods have been used to learn Gaussian mixtures. For example in (Richardson and Green, 1997), a reversible jump Markov Chain Monte Carlo (MCMC) method is proposed. There, the MCMC allows jumps between parameter spaces of different dimensionality (i.e. parameter spaces for mixtures consisting of differing number of components). However, the experimental results reported in [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF] indicate that such sampling methods are rather slow as compared to constructive maximum likelihood algorithms. It is reported that about 160 'sweeps' per second are performed on a SUN Sparc 4 workstation. The experiments involve 200.000 sweeps, resulting in about 20 minutes run time. Although it is remarked that the 200.000 sweeps are not needed for reliable results, this contrasts sharply with the 2.8 seconds and 5.7 seconds run time (allowing respectively about 480 and 960 sweeps) of the standard EM and our greedy EM in a similar experimental setting executed also on a SUN Sparc 4 workstation.

Conclusions

We proposed an efficient greedy method to learn mixtures of Gaussian densities that has run time O(k 2 n + kmn) for n data points, a final number of k mixture components and m candidates per component.

The experiments we conducted have shown that our greedy algorithm generally outperforms the k-means initialized EM. In experiments on artificial data we observed that SMEM performs in general comparable to the new method in terms of quality of the solutions found. The results on natural data were comparable to those obtained with SMEM, if SMEM was initialized with the k-means algorithm. An important benefit of the greedy method, both compared to SMEM and 'normal' EM, is the production of a sequence of mixtures, which obviates the need for initialization and facilitates model selection. As compared to VL we note: (i) The O(n 2 k 2 ) time complexity has been reduced by a factor n. (ii) The somewhat arbitrary choice of spherical candidate components with fixed variance has been replaced by a search for candidate components that depends on the data and the current mixture. (iii) Experiments suggest that if the methods yield different performance, then the new method generally outperforms the old one.

Software implementing the algorithm in Matlab is available from the first author via e-mail.

A

A Lower-bound for fast partial EM

In the insertion step we have a two component mixture problem. The first component is the mixture f k which we keep fixed in the partial EM steps. The second component is the new component φ which has a mixing weight α. The initial k + 1 component mixture is given by

f k+1 = αφ + (1 -α)f k .
Using the maximization-maximization view on EM [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF], we can define a lower-bound F on the actual log-likelihood:

i log f k+1 (x i ) ≥ F (Q, θ) = i log f k+1 (x i ) -D KL (Q i p(• | x i )), ( 14
)
where D KL denotes Kullback-Leibler divergence and p(•|x i ) denotes the posterior probability on the mixture components. The distributions Q i are the 'responsibilities' used in EM. We write q i for the responsibility of the new component and (1 -q i ) for the responsibility of f k for data point x i . The above lower bound can be made to match the log-likelihood by setting the responsibilities equal to the posteriors p(s | x i ). However, in our case, when considering a candidate based on component j, we fix the Q i for all points outside A j such that ∀x i ∈ A j : q i = 0. It is easy to see from ( 14), that the Q i for which x i ∈ A j should match the posterior probability on the new component in order to maximize F . Note that the more the assumption that data points outside A j have zero posterior on the new component is true, the tighter the lower bound becomes. If the data outside A j indeed has zero posterior on the new component the lower-bound equals the log-likelihood after the partial E-step.

For the M-step we use the alternative decomposition of F :

F (Q, θ) = i E Qi log p(x i , s) + H(Q i ), (15)
where we used the variable s to range over the two mixture components: f k and φ. Note that the entropies H(Q i ) remain fixed in the M-step since they do not depend on θ. The first term can be written as:

i q i [log φ(x i ) + log α] + (1 -q i )[log f k (x i ) + log (1 -α)] (16) = i∈Aj q i [log φ(x i ) + log α] + i (1 -q i )[log f k (x i ) + log (1 -α)]. (17)
Note that in (17) the new mixture component φ only occurs in terms for data in A j . Setting the derivative of (17) to zero w.r.t. the parameters of the new component (mixing weight α, mean m and covariance matrix C) gives the update equations in Section 4.3. Note that from (17) F can also be evaluated using only the points in A j if we stored the complete data log-likelihood under f k in advance.

Figure 1 :

 1 Figure 1: Images (a)-(d) illustrate the construction of a 4-component mixture distribution.

Figure 3 :

 3 Figure 3: Several Brodatz textures, the white squares indicate the patch size.

  Average conditional entropy for different values of k for the greedy approach, k-means initialized EM and SMEM and VL. The bottom line gives the conditional entropy if the clusters are totally uninformative. Bold face is used to identify the method with minimum conditional entropy for each k.

	k	2	3	4	5	6
	greedy	0.20 0.34 0.48 0.53 0.61
	k-means 0.19 0.46 0.74 0.80 0.94
	SMEM	0.21 0.22 0.36 0.56 0.68
	VL	0.93 1.52 2.00 2.29 2.58
	uniform	1	1.59	2	2.32 2.58
	Figure 4:

Acknowledgments We would like to thank Aristidis Likas for useful discussions and a referee for pointing out the confidence bounds used in (Smola and Schölkopf, 2000). This research is supported by the Technology Foundation STW (project nr. AIF4997) applied science division of NWO and the technology programme of the Dutch Ministry of Economic Affairs.


Aristidis Likas

Nikos Vlassis

Jakob J Verbeek
email: jverbeek@science.uva.nl

The global k-means clustering algorithm

Keywords: Clustering, K-means algorithm, Global optimization, k-d trees, Data mining

We present the global k-means algorithm which is an incremental approach to clustering that dynamically adds one cluster center at a time through a deterministic global search procedure consisting of N (with N being the size of the data set) executions of the k-means algorithm from suitable initial positions. We also propose modifications of the method to reduce the computational load without significantly affecting solution quality. The proposed clustering methods are tested on well-known data sets and they compare favorably to the k-means algorithm with random restarts.

Introduction

A fundamental problem that frequently arises in a great variety of fields such as pattern recognition, image processing, machine learning and statistics is the clustering problem [?]. In its basic form the clustering problem is defined as the problem of finding homogeneous groups of data points in a given data set. Each of these groups is called a cluster and can be defined as a region in which the density of objects is locally higher than in other regions.

The simplest form of clustering is partitional clustering which aims at partitioning a given data set into disjoint subsets (clusters) so that specific clustering criteria are optimized. The most widely used criterion is the clustering error criterion which for each point computes its squared distance from the corresponding cluster center and then takes the sum of these distances for all points in the data set. A popular clustering method that minimizes the clustering error is the k-means algorithm. However, the k-means algorithm is a local search procedure and it is wellknown that it suffers from the serious drawback that its performance heavily depends on the initial starting conditions [?]. To treat this problem several other techniques have been developed that are based on stochastic global optimization methods (eg. simulated annealing, genetic algorithms). However, it must be noted that these techniques have not gained wide acceptance and in many practical applications the clustering method that is used is the k-means algorithm with multiple restarts [?].

In this paper we propose the global k-means clustering algorithm, which constitutes a deterministic effective global clustering algorithm for the minimization of the clustering error that employs the k-means algorithm as a local search procedure. The algorithm proceeds in an incremental way: to solve a clustering problem with M clusters, all intermediate problems with 1, 2, . . . , M -1 clusters are sequentially solved. The basic idea underlying the proposed method is that an optimal solution for a clustering problem with M clusters can be obtained using a series of local searches (using the k-means algorithm). At each local search the M -1 cluster centers are always initially placed at their optimal positions corresponding to the clustering problem with M -1 clusters. The remaining M -th cluster center is initially placed at several positions within the data space. Since for M = 1 the optimal solution is known, we can iteratively apply the above procedure to find optimal solutions for all k-clustering problems k = 1, . . . , M . In addition to effectiveness, the method is deterministic and does not depend on any initial conditions or empirically adjustable parameters. These are significant advantages over all clustering approaches mentioned above.

In the following section we start with a formal definition of the clustering error and a brief description of the k-means algorithm and then we describe the proposed global k-means algorithm. Section 3 describes modifications of the basic method that require less computation at the expense of being slightly less effective. Section 4 provides experimental results and comparisons with the k-means algorithm with multiple restarts. Finally Section 5 provides conclusions and describes directions for future research.

The global k-means algorithm

Suppose we are given a data set X = {x 1 , . . . , x N }, x n ∈ R d . The M -clustering problem aims at partitioning this data set into M disjoint subsets (clusters) C 1 , . . . , C M , such that a clustering criterion is optimized. The most widely used clustering criterion is the sum of the squared Euclidean distances between each data point x i and the centroid m k (cluster center) of the subset C k which contains x i . This criterion is called clustering error and depends on the cluster centers m 1 , . . . , m M :

E(m 1 , . . . , m M ) = N i=1 M k=1 I(x i ∈ C k )||x i -m k || 2 (1)
where I(X) = 1 if X is true and 0 otherwise. The k-means algorithm finds locally optimal solutions with respect to the clustering error. It is a fast iterative algorithm that has been used in many clustering applications. It is a point-based clustering method that starts with the cluster centers initially placed at arbitrary positions and proceeds by moving at each step the cluster centers in order to minimize the clustering error. The main disadvantage of the method lies in its sensitivity to initial positions of the cluster centers. Therefore, in order to obtain near optimal solutions using the k-means algorithm several runs must be scheduled differing in the initial positions of the cluster centers.

In this paper, the global k-means clustering algorithm is proposed, which constitutes a deterministic global optimization method that does not depend on any initial parameter values and employs the k-means algorithm as a local search procedure. Instead of randomly selecting initial values for all cluster centers as is the case with most global clustering algorithms, the proposed technique proceeds in an incremental way attempting to optimally add one new cluster center at each stage.

More specifically, to solve a clustering problem with M clusters the method proceeds as follows. We start with one cluster (k = 1) and find its optimal position which corresponds to the centroid of the data set X. In order to solve the problem with two clusters (k = 2) we perform N executions of the k-means algorithm from the following initial positions of the cluster centers: the first cluster center is always placed at the optimal position for the problem with k = 1, while the second center at execution n is placed at the position of the data point x n (n = 1, . . . , N ). The best solution obtained after the N executions of the k-means algorithm is considered as the solution for the clustering problem with k = 2. In general, let (m * 1 (k), . . . , m * k (k)) denote the final solution for k-clustering problem. Once we have found the solution for the (k -1)clustering problem, we try to find the solution of the k-clustering problem as follows: we perform N runs of the k-means algorithm with k clusters where each run n starts from the initial state (m * 1 (k -1), . . . , m * (k-1) (k -1), x n ). The best solution obtained from the N runs is considered as the solution (m * 1 (k), . . . , m * k (k)) of the k-clustering problem. By proceeding in the above fashion we finally obtain a solution with M clusters having also found solutions for all k-clustering problems with k < M . The latter characteristic can be advantageous in many applications where the aim is also to discover the 'correct' number of clusters. To achieve this, one has to solve the k-clustering problem for various numbers of clusters and then employ appropriate criteria for selecting the most suitable value of k [?]. In this case the proposed method directly provides clustering solutions for all intermediate values of k, thus requiring no additional computational effort.

In what may be a concern of computational complexity, the method requires N executions of the k-means algorithm for each value of k (k = 1, . . . , M ). Depending on the available resources and the values of N and M , the algorithm may be an attractive approach, since, as experimental results indicate, the performance of the method is excellent. Moreover, as we will show later, there are several modifications that can be applied in order to reduce the computational load.

The rationale behind the proposed method is based on the following assumption: an optimal clustering solution with k clusters can be obtained through local search (using k-means) starting from an initial state with

• the k -1 centers placed at the optimal positions for the (k -1)-clustering problem and

• the remaining k-th center placed at an appropriate position to be discovered.

This assumption seems very natural: we expect that the solution of the k-clustering problem to be reachable (through local search) from the solution of (k -1)-clustering problem, once the additional center is placed at an appropriate position within the data set. It is also reasonable to restrict the set of possible initial positions of the k-th center to the set X of available data points. It must be noted that this is a rather computational heavy assumption and several other options (examining fewer initial positions) may also be considered. The above assumptions are also verified experimentally, since in all experiments (and for all values of k) the solution obtained by the proposed method was at least as good as that obtained using numerous random restarts of the k-means algorithm. In this spirit, we can cautiously state that the proposed method is experimentally optimal (although it is difficult to prove theoretically).

Speeding-up execution

Based on the general idea of the global k-means algorithm, several heuristics can be devised to reduce the computational load without significantly affecting the quality of the solution. In the following subsections two modifications are proposed, each one referring to a different aspect of the method.

The fast global k-means algorithm

The fast global k-means algorithm constitutes a straightforward method to accelerate the global k-means algorithm. The difference lies in the way a solution for the k-clustering problem is obtained, given the solution of the (k -1)-clustering problem. For each of the N initial states

(m * 1 (k -1), . . . , m * (k-1) (k -1),
x n ) we do not execute the k-means algorithm until convergence to obtain the final clustering error E n . Instead we compute an upper bound E n ≤ E -b n on the resulting error E n for all possible allocation positions x n , where E is the error in the (k -1)clustering problem and b n as defined in Eq. ( 2). We then initialize the position of the new cluster center at the point x i that minimizes E n , or equivalently that maximizes b n , and execute the k-means algorithm to obtain the solution with k clusters. Formally we have

b n = N j=1 max(d j k-1 -||x n -x j || 2 , 0), ( 2
) i = arg max n b n (3)
where d j k-1 is the squared distance between x j and the closest center among the k -1 cluster centers obtained so far (ie., center of the cluster where x j belongs). The quantity b n measures the guaranteed reduction in the error measure obtained by inserting a new cluster center at position x n .

Suppose the solution of the (k -1)-clustering problem is (m * 1 (k -1), . . . , m * (k-1) (k -1)) and a new cluster center is added at location x n . Then the new center will allocate all points x j whose squared distance from x n is smaller than the distance d j k-1 from their previously closest center. Therefore, for each such data point x j the clustering error will decrease by

d j k-1 -|x n -x j | 2 .
The summation over all such data points x j provides the quantity b n for a specific insertion location x n . Since the k-means algorithm is guaranteed to decrease the clustering error at each step, E -b n upper bounds the error measure that will be obtained if we run the algorithm until convergence after inserting the new center at x n (this is the error measure used in the global k-means algorithm).

Experimental results (see next section) suggest that using the data point that minimizes this bound leads to results almost as good as those provided by the global k-means algorithm. Moreover, the cluster insertion procedure can be efficiently implemented by storing in a matrix all pairwise squared distances between points when the algorithm starts, and using this matrix for directly computing the upper bounds above. A similar 'trick' has been used in the related problems of greedy mixture density estimation using the EM algorithm [?] and principal curve fitting [?].

Finally, we may still apply this method as well as the global k-means algorithm when we do not consider every data point x n (n = 1, . . . , N ) as possible insertion position for the new center, but use only a smaller set of appropriately selected insertion positions. A fast and sensible choice for selecting such a set of positions based on k-d trees is discussed next.

Initialization with k-d trees

A k-d tree [?, ?] is a multi-dimensional generalization of the standard one-dimensional binary search tree, that facilitates storage and search over k-dimensional data sets. A k-d tree defines a recursive partitioning of the data space into disjoint subspaces. Each node of the tree defines a subspace of the original data space and, consequently, a subset containing the data points residing in this subspace. Each nonterminal node has two successors, each of them associated with one of the two subspaces obtained from the partitioning of the parent space using a cutting hyperplane. The k-d tree structure was originally used for speeding up distance-based search operations like nearest neighbors queries, range queries, etc.

In our case we use a variation of the original k-d tree proposed in [?]. There, the cutting hyperplane is defined as the plane that is perpendicular to the direction of the principal component of the data points corresponding to each node, therefore the algorithm can be regarded as a method for nested (recursive) principal component analysis of the data set. The recursion usually terminates if a terminal node (called bucket) is created containing less than a prespecified number of points b (called bucket size) or if a prespecified number of buckets have been created. It turns out that, even if the algorithm is not used for nearest neighbor queries, merely the construction of the tree provides a very good preliminary clustering of the data set. The idea is to use the bucket centers (which are fewer than the data points) as possible insertion locations for the algorithms presented previously.

In Figure 1 average performance results are shown on 10 data sets each one consisting of 300 data points drawn from the same mixture of 15 Gaussian components. The components of the Gaussian mixture are well separated and exhibit limited eccentricity.

We compare the results of three methods to the clustering problem with k = 15 centers: (i) The dashed line depicts the results when using the fast global k-means algorithm with all data points constituting potential insertion locations. The average clustering error over the 10 data sets is 15.7 with standard deviation 1.2. (ii) The solid line depicts results when the standard k-means algorithm is used: one run for each data set was conducted. At each run the 15 cluster centers were initially positioned to the centroids of the buckets obtained from the application of the k-d tree algorithm until 15 buckets were created. The average clustering error over the 10 data sets is 24.4 with standard deviation 9.8. (iii) The solid line (with error bars indicating the standard deviation from the mean value) shows the results when using the fast global k-means algorithm with the potential insertion locations constrained by the centroids of the buckets of a k-d tree. On the horizontal axis we vary the number of buckets for the k-d tree of the last method.

We also computed the 'theoretical' clustering error for each data set, ie., the error computed by using the true cluster centers. The average error value over the 10 data sets was 14.9 with standard deviation 1.3. These results were too close to the results of the standard fast global k-means to include them in the figure .  We can conclude from this experiment that (a) the fast global k-means approach gives rise to performance significantly better than when starting with all centers at the same time initialized using the k-d tree method, and (b) restricting the insertion locations for the fast global k-means to those given by the k-d tree (instead of using all data points) does not significantly degrade performance if we consider a sufficiently large number of buckets in the k-d tree (in general larger than the number clusters).

Obviously, it is also possible to employ the above presented k-d tree approach with the global k-means algorithm.

Experimental results

We have tested the proposed clustering algorithms on several well-known data sets, namely the iris data set [?], the synthetic data set [?] and the image segmentation data set [?]. In all data sets we conducted experiments for the clustering problems obtained by considering only feature vectors and ignoring class labels. The iris data set contains 150 four-dimensional data points, the synthetic data set 250 two-dimensional data points and the for the image segmentation data set we consider 210 six-dimensional data points obtained through PCA on the original 18-dimensional data points. The quality of the obtained solutions was evaluated in terms of the values of the final clustering error.

For each data set we conducted the following experiments:

• one run of the global k-means algorithm for M = 15.

• one run of the fast global k-means algorithm for M = 15.

• the k-means algorithm for k = 1, . . . , 15. For each value of k, the k-means algorithm was executed N times (where N is the number of data points) starting from random initial positions for the k centers, and we computed the minimum and average clustering error as well as its standard deviation.

For each of the three data sets the experimental results are displayed in Figures 2, 3 and 4 respectively. Each figure plot displays the clustering error value as a function of the number of clusters. It is clear that the global k-means algorithm is very effective providing in all cases solutions of equal or better quality with respect to the k-means algorithm. In terms of the fast version of the algorithm, it is very encouraging that, although executing significantly faster, it provides solutions of excellent quality, comparable to those obtained by the original method. Therefore, it constitutes a very efficient algorithm, both in terms of solution quality and computational complexity and can run much faster if k-d trees are employed as explained in the previous section.

In the following subsections we provide extensive experimental results comparing the fast versions of the algorithm to the conventional k-means algorithm with random initialization.

Texture segmentation

In this clustering experiment the objective is to cluster 16 × 16 pixel image patches extracted from a set of 37 Brodatz texture images [?]. Each complete texture image consists of 256 × 256 pixels and 500 patches per texture were extracted by randomly selecting 16 × 16 windows. It is expected that patches originating from the same texture image to form an individual cluster (or clusters).

For every number of textures k = 2, . . . , 6, we randomly constructed 100 data sets. Each data set was created by first randomly selecting k of the 37 textures and then selecting 200 patches for each texture, resulting in 200k patches per data set. All reported results are averages over the 100 data sets.

We compared performance of the three algorithms: (1) k-means initialized using a uniformly selected random subset of the data, (2) fast global k-means, (3) fast global k-means with the insertion locations limited to the top 2k nodes of the corresponding kd-tree. To evaluate the different methods, we considered the mean squared clustering error (MSE) of the patches to their closest mean. In a first experiment the number clusters was considered equal to the number of textures selected to create a data set. The results are given in Figure 5, where we also provide the corresponding execution times (in seconds). Bold values in the Tables indicate the method with best performance in terms of the clustering error.

We also conducted a second series of experiments where we used twice as many clusters as textures. The results are displayed in Figure 6.

It can be observed that (on average) the k-means algorithm with random initialization gives the worst results for almost every number of clusters and textures. It is also interesting to note that the fast global k-means that uses the top 2k nodes of the kd-tree as insertion candidates is not only much faster than the generic fast global k-means algorithm, but it also provides slightly better results in terms of the clustering error.

Artificial data

In this subsection we provide more extensive comparative experimental results using artificially created data sets. The purpose is to compare the randomly initialized k-means algorithm with the fast global k-means algorithm that uses a the top 2k nodes of the corresponding kd-tree as candidate insertion locations.

The data have been drawn from randomly generated Gaussian mixtures [?] and we varied the number of sources (mixture components) k, the dimensionality d of the data space and the 'separation' c between the sources of the mixture. Following [?], the separation c of a Gaussian mixture satisfies:

∀ i =j : µ i -µ j 2 ≥ c • max {i,j} {trace(C i ), trace(C j )}, (4)
where the µ's and C's denote respectively the means and covariance matrices of the mixture components. In our experiments we considered mixtures having a separation in the range between 0.5 and 1 and these values correspond exclusively to weakly separated clusters. The number of data points in each data set was 50k, where k is the number of sources. Some example data sets are shown in Figure 7.

We have considered 120 problems corresponding to any of the following combination of values: k = {2, 7, 12, 17, 22}, d = {2, 4, 6, 8}, c = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0}. For each problem 10 data sets were created. The results are displayed in Figures 8 and9.

On every data set, the 'greedy' (fast global k-means with k-d tree initialization) algorithm was applied first. Then the randomly initialized k-means algorithm was applied as many times as possible in the run time of the greedy algorithm. To evaluate the quality of the solutions found for a specific data set, first the mean clustering error µ is computed for the runs performed using the k-means algorithm. The minimum clustering error value corresponding to the runs with the k-means algorithm is provided in the 'min.' columns as (1 -min /µ) × 100, and the standard deviation in the σ column as σ/µ × 100. The 'gr.' column provides the clustering error E gr obtained with the fast global k-means algorithm in the form: (1 -E gr /µ) × 100. Finally, the 'trials' column indicates how many runs of the k-means algorithm were allowed in the run time of the greedy algorithm. It must be noted that each row of the table displays the averaged results over 10 data sets constructed for the specific values of k, d and c. We used bold values in the Tables to indicate whether the clustering error of the greedy algorithm or the minimum error achieved for the runs with the k-means algorithm was the smallest.

It is clear from the experiments that the benefit of the greedy method becomes larger if there are more clusters, the seperation becomes larger and the dimensionality gets smaller. It is striking that in almost all cases the greedy algorithm gives better results. In cases where the greedy method is not superior, both methods yield results relatively close to the average k-means result. One thing to note is that the clusters in the created data sets are poorly seperable. It is also interesting to note that the number of trials allowed for the random k-means algorithm grows only slowly as the number of clusters increases. From top to bottom, although the number of clusters has increased by a factor 10 the number of trials however has not even doubled on average. The inspection of the number of trials as a function of the k, d and c indicates dependence on the number of clusters only. As Figure 10 indicates, both the average number of trials as well as the corresponding variance increase as the number of clusters increases.

Matlab implementations of the fast global k-means and the k-d tree building algorithms can be downloaded from http://www.science.uva.nl/research/ias.

Discussion and conclusions

We have presented the global k-means clustering algorithm, which constitutes a deterministic clustering method providing excellent results in terms of the clustering error criterion. The method is independent of any starting conditions and compares favorably to the k-means algorithm with multiple random restarts. The deterministic nature of the method is particularly important in cases where the clustering method is used either to specify initial parameter values for other methods (for example RBF training) or constitutes a module in a more complex system. In such a case we can be almost certain that the employment of the global k-means (or any of the fast variants) will always provide sensible clustering solutions. Therefore, one can evaluate the complex system and adjust critical system parameters without having to worry for dependence of system performance on the clustering method employed.

Another advantage of the proposed technique is that in order to solve the M -clustering problem, all intermediate k-clustering problems are also solved for k = 1, . . . , M . This may prove useful in many applications where we seek for the actual number of clusters and the k-clustering problem is solved for several values of k. We have also developed the fast global k-means algorithm, which significantly reduces the required computational effort, while at the same time providing solutions of almost the same quality.

We have also proposed two modification of the method that reduce the computational load without significantly affecting solution quality. These methods can be employed to find solutions to clustering problems with thousands of high-dimensional points and one of our primary aims is to test the techniques on large scale data mining problems.

Another direction of future work is related with the use of parallel processing for accelerating the proposed methods since, for every k, the N executions of the k-means algorithm are independent and can be performed in parallel. Another research direction concerns the application of the proposed method to other types of clustering (for example fuzzy clustering), as well as to topographic methods like SOM. Moreover, an important issue that deserves further study is related with the possible development of theoretical foundations for the assumptions behind the method. Finally, it is also possible to employ the global k-means algorithm as a method for providing effective initial parameter values for RBF networks and data modeling problems using Gaussian mixture models and compare the effectiveness of the obtained solutions with other training techniques for Gaussian mixture models [?, ?].

Figure 1 :

 1 Figure 1: Performance results for data drawn from a Gaussian mixture with 15 components.

Figure 2 :

 2 Figure 2: Performance results for the Iris data set.

Figure 3 :

 3 Figure 3: Performance results for the Synthetic data set.

Figure 4 :

 4 Figure 4: Performance results for the Image Segmentation data set.

Figure 5 :

 5 Figure 5: Results for the Texture Segmentation problem using as many clusters as textures.

Figure 6 :

 6 Figure 6: Results for the Texture Segmentation problem using twice as many clusters as textures.

Figure 7 :

 7 Figure 7: Some example 2D data sets with five sources and separation (from left to right) 0.5, 0.75 and 1.

Figure 8 :

 8 Figure 8: Experimental results on artificial data sets with d=2 and d=4.

Figure 9 :

 9 Figure9: Experimental results on artificial data sets with d=6 and d=8.

Figure 10 :

 10 Figure 10: The number of allowed trials for the randomly initialized k-means as a function of the number of clusters.

Acknowledgements

N. Vlassis and J. J. Verbeek are supported by the Dutch Technology Foundation STW project AIF 4997.

The global k-means clustering algorithm Aristidis Likas, Nikos Vlassis and Jakob J. Verbeek

Summary

In this paper the global k-means algorithm is presented which is an incremental approach to clustering that dynamically adds one cluster center at a time through a deterministic global search procedure consisting of N (with N being the size of the data set) executions of the k-means algorithm from suitable initial positions. The basic idea underlying the proposed method is that an optimal solution for a clustering problem with M clusters can be obtained using a series of local searches (using the k-means algorithm). At each local search the M -1 cluster centers are always initially placed at their optimal positions corresponding to the clustering problem with M -1 clusters. The remaining M -th cluster center is initially placed at several positions within the data space. Since for M = 1 the optimal solution is known, we can iteratively apply the above procedure to find optimal solutions for all k-clustering problems k = 1, . . . , M . The proposed method is deterministic, does not depend on any initial positions for the cluster center and does not contain any empirically adjustable parameters, thus eliminating all problems characterizing the k-means algorithm and its stochastic extensions.

In addition we present modifications of the method to reduce the computational load without significantly affecting solution quality. The first modification called the fast global k-means algorithm defines a fast computed bound on the clustering error that is used instead of local searches. The second modification is related with the partitioning of the data space using a k-d tree structure, in order to reduce the number of examined insertion positions for a new cluster center.

The proposed clustering methods are tested on well-known data sets and they compare favorably to the k-means algorithm with random restarts. In addition, the modified versions lead to implementations that are very fast and exhibit almost equal performance.