Hartigan's K-Means Versus Lloyd's K-Means - Is It Time for a Change?

1 Introduction

The goal of cluster analysis is to partition a given set of data items into clusters such that similar items are assigned to the same cluster whereas dissimilar ones are not. Perhaps the most popular clustering formulation is K-means in which the goal is to maximize the expected similarity between data items and their associated cluster centroids. The classical optimization heuristic for K-means is Lloyd's algorithm [Lloyd, 1982; MacQueen, 1967; Forgy, 1965]. In this algorithm, one typically starts with a random partition of the data and then proceeds by alternating between two steps - an assignment step where each item is assigned with the cluster represented by the nearest centroid; and an update step in which the clusters' centroids are updated given the partition obtained by the previous step. This algorithm is proven to converge in the sense that after a finite number of such trials the assignments of the data items no longer change.

Awell known problemwith this procedure is that it is sensitive to the initial partition being used. Namely, the algorithm typically converges to a fixed-point that represents only a local minimum in the space of all possible partitions, and various heuristics have been proposed in the literature to address this issue (cf. [Bradley and Fayyad, 1998]). However,

Koby Crammer

Dept. Elec. Engineering The Technion Haifa, Israel koby@ee.technion.ac.il

Telgarsky and Vattani [Telgarsky and Vattani, 2010] have recently highlighted that the problem of converging to a local minima is less severe if an alternative optimization heuristic - termed Hartigan's method in [Telgarsky and Vattani, 2010] - is being used. Specifically, in this alternative heuristic, a single data item is being examined and then optimally re-assigned by the algorithm [Hartigan, 1975]. Interestingly, since this algorithm takes into account the motion of the centroids resulting from the re-assignment step, a data item is not always assigned with the cluster with the nearest centroid. As a result, the set of local minima of Hartigan's algorithm is a strict subset of those of Lloyd's method, implying that the algorithm is less sensitive to the choice of the initial random partition of the data [Telgarsky and Vattani, 2010].

Here, we start by demonstrating how Hartigan's method can be easily applied with any Bregman divergence [Bregman, 1967]. Our formulation includes as special cases Hartigan's method with the Euclidean norm that was analyzed in detail in [Telgarsky and Vattani, 2010], as well as the sequential Information Bottleneck (sIB) algorithm, originally proposed in [Slonim, 2002]. In addition, we provide a systematic quantitative estimation of the difference between Hartigan's algorithm and Lloyd's algorithm. In particular, we consider the classical example of data generated by a mixture of Gaussians as well as several real world data. In all cases we characterize the number and quality of the local minima obtained by both algorithms as a function of various parameters of the problem. Our results reveal a substantial difference between both algorithms. Specifically, we characterize a wide range of problems for which any initial partition represents a local minimum for the Lloyd's algorithm. Nonetheless, for the same problems, the number of local minima associated with Hartigan's algorithm is small, and furthermore, these local minima correspond well with the true partition of the data.

2 Lloyd's Algorithm

Let X denote a finite set of data items, X = {x1, ..., xn}, each represented by some vector, vx ∈ Rm. Let C = {c1, ..., cK} denote some partition of X into K disjoint nonempty subsets, or clusters, and let n(c) denote the number of data items assigned with a cluster c. The standard quality measure of the partition C, known as the K-means formulation [MacQueen, 1967], is given by the mean distortion,

where d(·, ·) is a pre-specified non-negative distortion measure, and every c ∈ C is represented by a centroid vector in vc ∈ Rm. Given K < n, the goal is to find a partition C into K clusters that minimizes D(C).

The most popular optimization heuristic to this end is Lloyd's algorithm [Lloyd, 1982; MacQueen, 1967; Forgy, 1965]. In most applications, d(vx, vc) = 1 2 ‖vx - vc‖2 , where ‖ · ‖ denotes the Euclidean norm, and one starts from some random partition of the data into K clusters. The algorithm is then an iterative two-step process. In the assignment step, each data item is assigned with the cluster that corresponds to its nearest centroid. In the update step, the centroids are updated through (1), based on the new partition. It is easy to verify that each of these steps can only reduce D(C). Since D(C) ≥ 0 the algorithm is guaranteed to converge to a stable fixed-point that corresponds to a local minimum of the expected distortion.

An important feature of Lloyd's algorithm is that the parameters' updates are performed in batch mode; first, all the assignments are re-estimated and only then all centroids are updated. Hence, we will refer to this classical algorithm as batch K-means and to its fixed-point partitions as batch local minima.

A known variant of this algorithm is the online K-means (cf. [Har-peled and Sadri, 2005]). In this variant, in each trial a single item is selected at random and re-assigned to the cluster with the nearest centroid; next, the relevant centroids are updated accordingly. Clearly, this algorithm as well converges to a fixed-point partition that may depend on the initial partition being used and on the particular item selections per trial. We will henceforth term the fixed-point partitions of this algorithm as online local minima.

Observation 2.1 C is a batch local minimum or an online local minimum, if and only if ∀ x ∈ X, ∀c ∈ C, d(vx, vc(x)) ≤ d(vx, vc) .

Corollary 2.2 For any given data, the set of batch local minima is identical to the set of online local minima.

Thus, we will henceforth refer to these local minima simply as Lloyd's local minima. Notice, that the above corollary do not imply that for a given initial partition the batch algorithm and the online algorithm converge to the same local minimum. In fact, the online algorithm is more stochastic in nature since for a given initial partition it may converge to various local minima, depending on the random selections it performs per trial.

3 Hartigan's Method

Hartigan's algorithm is an alternative heuristic to Lloyd's algorithm, also aiming to optimize the K-means cost function - (1). The basic idea of this algorithm is rather simple. Given a partition C, the algorithm first draws at random a single item, x, out of its cluster; e.g., if x ∈ c, after drawing x from c it is represented as a singleton cluster with centroid vx, while vc is updated accordingly as the average vx′ taken across all x′ ∈ c, x′ ̸= x. Next, the algorithm finds the optimal cluster c∗ to which x should be re-assigned, in terms of minimizing D, assigns x with c∗ , and updates vc∗ accordingly. By construction, in each trial the algorithm can only decrease D, hence it is guaranteed to converge to a fixed-point partition after a finite number of trials. In the following we refer to these fixed-point partitions as Hartigan's local minima.

To gain further intuition we examine the impact over D due to merging the singleton cluster consisted of x to one of the K clusters. Let c+ denote the cluster c after the addition of x. Then, obviously, x /∈ c, x ∈ c+, n(c+) = n(c) + 1 , and it is easy to verify that the increment in D due to this merger is given by

Thus, ∆D(x, c) is consisted of a trade-off between two terms. First, we would like to assign x to c such that d(vx, vc+) will be minimized; however, at the same time, we would like the resulting new centroid, vc+ , to remain a relatively good representative for all x′ ∈ c, namely to minimize the second term in (2).

To summarize, although reminiscent of online K-means, Hartigan's clustering is different in two important aspects. First, when considering to what cluster to re-assign x, one first draws x out of its original cluster, c(x), and updates vc(x) accordingly. Thus, while in online (and batch) Lloyd's algorithm there is some bias in favor of re-assigning x back to c(x) - since the examined vc(x) includes the impact of vx - in Hartigan's algorithm this bias vanishes. Second, Hartigan's algorithm considers a wider horizon of the impact of assigning x to c ∈ C, as it not only considers how similar is x to c, but further takes into account the general impact over the cost function, D, that will result from updating vc to vc+. Thus, Hartigan's algorithm will assign x to c such that D will be minimized, and as implied by (2), this is not necessarily the cluster with the nearest centroid. A recent insightful theoretical analysis of the difference between Hartigan's algorithm and Lloyd's algorithm is given in [Telgarsky and Vattani, 2010].

4 Hartigan's Method with Bregman Divergence

A possible concern regarding Hartigan's algorithm is that its potential advantages come at the cost of higher complexity. In each trial the algorithm need to calculate ∆D(x, c) for all c ∈ C; assuming d(·, ·) is calculated in O(m) and using the straightforward derivation in (2), we find that the complexity of calculating each ∆D(x, c) is O(n(c)·m), hence the overall complexity of a single trial is O(n · m). However, as we show next, for an important family of distortion measures, the complexity of a single trial can be reduced to O(K · m), which is identical to the complexity of a single trial in the online variant of Lloyd's algorithm.

Recall that each trial of Hartigan's algorithm starts by drawing a particular x out of its cluster, c(x) and updating vc(x) accordingly. Denoting by c(x)- the cluster c(x) after x's removal, it is easy to verify that

hence, this part can be computed in O(m). The trial is completed by merging x to a cluster c and updating vc accordingly. Denoting by c+ the cluster obtained after adding x to c, it is again easy to verify that the updated centroid is given by

hence, this part as well can be computed in O(m). Thus, the remaining question is how to efficiently identify the cluster c∗ to which x should be assigned such that D will be minimized. In our context, this maps to calculating ∆D(x, c) for all c ∈ C. Next, we provide a closed-form expression for ∆D(x, c) for any distortion measure that corresponds to a Bregman divergence [Bregman, 1967].

Definition 4.1 Let Y ⊂ Rm be some closed, convex set. Let F : Y → R be a strictly convex function. Then, for any v, w ∈ Y the Bregman divergence [Bregman, 1967] is defined as 1

Thus, dF measures the difference between F and its firstorder Taylor expansion about w, evaluated at v. Although dF does not necessarily satisfy the triangle inequality nor symmetry, it is always non-negative and equals to zero if and only if its two arguments are identical. The Bregman divergence generalizes some commonly studied distortion measures. For example, for F(v) = 0.5‖v‖2 we have dF (v|w) = 0.5‖v - w‖2 , namely the squared Euclidean distance. Alternatively, if Y = Rm+ and F(v) = ∑mr=1 v(r) log(v(r) ), then dF (v|w) is the unnormalized KL divergence, dF (v|w) = ∑mr=1 (v(r) log v(r) w(r) +w(r) -v(r) ). In the special case where Y is the m-th dimensional simplex, v and w are normalized distributions, F(v) is Shannon entropy, and dF (v|w) reduces to the conventional KL divergence [Cover and Thomas, 2006].

Proposition 4.2 Given the above notations, if d(·, ·) is a Bregman divergence, then

Thus, the effect of merging x to c is related not only to the similarity of vx and vc+, but also to the induced change on vc, as quantified by the second term in (6) (in this context, see also [Telgarsky and Dasgupta, 2012]). In particular, using (6), the complexity of computing ∆D(x, c) is O(m), as required. A Pseudo-code for a general Hartigan's algorithm is given in Fig. 1.

A few special cases are worth considering in detail. First, if Y is the simplex of Rm and dF (v|w) is the KL divergence,

1 In principle, a Bregman function F needs to fulfill some additional technical constraints which we omit here to ease the presentation (see [Censor and Zenios, 1997] for details).

Input

data items: X = {vx1 , ..., vxn } ⊂ Rm, K.

Output

A Hartigan's fixed point partition of X into K clusters. Initialization

C ← random partition of X into K clusters.

Main Loop

While not Done

Done ← TRUE .

Scan X by some random order and ∀x ∈ X Remove x from c(x) and update vc(x). c∗ = argminc∈C∆D(x, c). If c∗ ̸= c(x), Done ← FALSE. Merge x into c∗ and update vc∗ .

Figure 1: Pseudo-code of the Hartigan's algorithm. If d(·, ·) is a Bregman divergence, ∆D(x, c) can be computed efficiently via (6); otherwise, ∆D(x, c) can be computed directly from (2).

the derivation above yields precisely the sequential Information Bottleneck (sIB) algorithm, that was originally proposed in the context of document clustering [Slonim et al., 2002; Slonim, 2002]. However, since dF can be applied to unnormalized nonnegative vectors in Rm, our derivation further yields an extension of the sIB algorithm to cluster unnormalized vectors, which we intend to explore in future work. Another important special case emerges when dF (v|w) is the squared Euclidean distance, d(vx, vc) = 0.5‖vx - vc‖2 . In this case Proposition 4.2 reduces to

which is precisely the algorithm proposed in [Hartigan, 1975], that was more recently analyzed in detail in [Telgarsky and Vattani, 2010]. Due to the wide use of this distortion measure in cluster analysis, in the remaining of this work we focus on this special case.

5 Theoretical Analysis

For completeness, we first repeat Theorem 2.2 from [Telgarsky and Vattani, 2010].

Theorem 5.1 (Telgarsky-Vattani) For any given data X = {x1, x2, . . . , xn} the set of Hartigan's local minima is a - possibly strict - subset of the set of Lloyd's local minima.

Corollary 5.2 Lloyd's algorithm can never improve a Hartigan's local minimum while Hartigan's algorithm might improve a Lloyd's local minimum.

Proof: The first part is proven by Theorem 5.1. The second part may be demonstrated using a simple example. Consider the following one-dimensional data items {-5, 0, 0, 0, 0, 0, 1}. For these data, the partition C (1) = {{-5, 0, 0, 0, 0, 0}, {1}} is a Lloyd's local minimum, with associated cost D ≈ 1.49. However, C (1) is not a Hartigan's minimum. Specifically, drawing {-5} out of its cluster, the Hartigan's algorithmwill assign it to the other cluster to obtain the partition C (2) = {{0, 0, 0, 0, 0}, {-5, 1}} with an improved cost of D ≈ 1.29; next, the Hartigan's algorithm can re-assign {1} to end up with a partition C (3) = {{1, 0, 0, 0, 0, 0}, {-5}} with a cost of D ≈ 0.06, which is the globally optimal partition in this example. 2

The example outlined in the proof of Corollary 5.2 demonstrates another important difference between Hartigan's algorithm and Lloyd's algorithm. Specifically, at each iteration Lloyd's algorithm effectively partitions the samples using a Voronoi diagram according to current centroids. Thus it is restricted to solely explore partitions that correspond to Voronoi diagrams, and is prohibited from crossing through partitions such as C (2) which allows the Hartigan's algorithm to escape the local minimum. A detailed theoretical analysis in this context is given in [Telgarsky and Vattani, 2010]. In addition, since the number of Hartigan's minima is always bounded by the number of Lloyd's minima, it is reasonable to expect that the Hartigan's algorithm will typically be less sensitive to different initializations, namely, more stable with respect to the initial conditions. Moreover, given the output of Lloyd's algorithm, one can only gain by using it as the input to the Hartigan's algorithm, while the other direction is vain.

Given the above qualitative observations, we now outline two results that provide additional quantitative insights. In particular, we show that there exists a wide range of settings where any possible partition C is a local minimum of Lloyd's algorithm, rendering this algorithm useless for those settings. The crux of the analysis is based on Observation 2.1. In order to improve a partition C Lloyd's algorithm must find an element x such that the distance ||vx - vc(x)||2 is larger than ||vx - vc||2 for some cluster c. However, vc(x) is biased towards vx since the calculation of the centroid takes vx into account, thus vx tends to be closer to its own cluster centroid vc(x) than to any other centroid.

Lemma 5.3 Let A = {a1, a2, . . . , an1 } and B = {b1, b2, . . . , bn2 } be two sets of scalar values sampled i.i.d. from some distribution F with variance σ2 . Let va and vb be their respective mean values. If S is the random variable defined by S = (a1 - vb)2 - (a1 - va)2 , then E(S) =

2 Additional examples can be generated and illustrated via a Web interface provided by Telgarsky and Vattani in this context at http://cseweb.ucsd.edu/∼mtelgars/htv/.

Proof: Considering the random variable a1 - vb, since vb is the mean of i.i.d. samples out of F and a1 is independently sampled from F, we have E(a1 - vb) = 0 and V ar(a1 - vb) = (1 + 1 n2 )σ2 . Thus E((a1 - vb)2 ) = V ar(a1 - vb) = (1 + 1 n2 )σ2 . Similar arguments hold for a1 - va, but first let us rewrite it as n1-1 n1 (a1 -v′a) where v′a is the mean of A excluding a1. Therefore v′a is independent of a1, and E(a1 - v′a) = 0 and V ar(a1 - v′a) = (1+ 1 n1-1 )σ2 . Simple algebra yields E((a1 - va)2 ) = V ar(a1 - va) = (1- 1 n1 )σ2 . Thus we may conclude E(S) = ( 1 n1 + 1 n2 )σ2 .

Theorem 5.4 Let v1, v2, . . . , vn be n data items such that vi ∈ Rm1+m2 . Let us assume that for each vi the first m1 components were drawn from some multivariate distribution F1 while each of the remaining m2 components were drawn i.i.d. from some univariate distribution F2 with variance σ2 . Let C be a partition of the n vectors chosen with uniform distribution over all possible partitions. If LC denotes the event that C is a Lloyd's local minimum, then P(LC) → 1 with m2 → ∞.

Proof: Let us randomly pick a pair of clusters c1 and c2 with centroids vc1 and vc2 . Let vx be a randomly chosen vector in c1. We now explore the probability of vx being closer to vc2 than to vc1 by investigating the random variable S = ||vx - vc2 ||2 - ||vx - vc1 ||2 . Let Sj be the contribution of the j'th component of each vector to S, i.e., Sj = (vjx - vj c2 )2 - (vjx - vj c1 )2 and S = ∑m1+m2 j=1 Sj. The contribution of the first m1 components, ∑m1 j=1 Sj is a random variable with some distribution derived from F1, with some expectation µ1 and variance σ2 1. All the remaining Sj's for j > m1 are i.i.d. variables with some expectation µ2 and variance σ2 2. From Lemma 5.3, for j > m1 we have E(Sj|C) = ( 1 n(c1) + 1 n(c2) )σ2 , and since 1 n(c1) + 1 n(c2) > 2n we obtain µ2 > 2n σ2 . From this it follows that E(S) > µ1 +m2 2n σ2 and V ar(S) = σ2 1 + m2σ2 2 . The expectation grows faster with m2 than the standard deviation √σ2 1 + m2σ2 2 . Thus using Chebyshev's inequality [DeGroot and Schervish, 2002] the probability P(S < 0) can be made arbitrarily close to 0 for sufficiently large m2. If C is not a Lloyd's local minimum then necessarily there exists some vx whose distance to another cluster's centroid is smaller than its own. Thus the probability of such an event is bounded by n · K · P(S < 0), which concludes the proof.

We now proceed to show that this potential problem of local minima overflow cannot occur for Hartigan's algorithm. Proposition 5.5 proves that the frequency of Hartigan's minima among all possible partitions is bounded by 1K , irrespective of the number of dimensions, m.

Proposition 5.5 Let C be the set of all possible valid partitions of X = {x1, . . . , xn} into 1 ≤ k ≤ K clusters. Then, if no two partitions have the exact same cost, the number of Hartigan's local minima is bounded by |C| K .

Proof: Let R ⊂ C×C be a binary relation such that (Ci, Cj) ∈ R if and only if Ci and Cj are identical except in the cluster assignment of x1, i.e., the Hartigan's algorithm may cross from Ci to Cj by re-assigning x1. Clearly this is an equivalence relation therefore it partitions C into equivalence classes. An equivalence class may contribute at most one local minimum since we assume no two partitions have the same cost. Most equivalence classes size is exactly K. However an equivalence class may be smaller. This may only occur if all its partitions have less than K clusters. However, a partition with less than K clusters cannot be a Hartigan's minimum since it must have a cluster with at least 2 elements of which one can be extracted to form a new singleton cluster with reduced cost. Thus, since only equivalence classes of size K contribute Hartigan's minima and each contributes at most one, the total number of Hartigan's minima is bounded by |C| K .

6 Experimental Results

Systematic Exploration of Toy Examples: We start with a synthetic toy example for which the space of all possible partitions can be explored systematically. Specifically, in this example n = 16, K = 2, and m = 1, 2, . . . , 50. The 16 data points are sampled as follows. For the first dimension, 8 points are sampled from a normal distribution with µ1 = -1, σ = 1, while the remaining 8 points are sampled from a normal distribution with µ = 1, σ = 1. For each of the remaining m- 1 dimensions, all 16 points are sampled from a normal distribution with µ = 0, σ = 1. Hence, increasing m amounts to adding noise features, and according to Theorem 5.4 is expected to increase the number of Lloyd's minima. For each m value we drew 100 samples of 16 points each, and for each of these samples we explored the entire space of all 216 = 65, 536 partitions, where for simplicity permutations of the same partition and partitions with empty clusters were not removed.

In Fig. 2(a) we depict for each value of m the observed number of Lloyd's and Hartigan's local minima, normalized by the total number of partitions, 65, 536. As expected, the number of Lloyd's local minima is increasing as m increases, getting close to covering the entire partition space. Note that due to Corollary 2.2 these results apply to Lloyd's online algorithm as well. In contrast, the number of Hartigan's local minima seems to converge to a small constant.

In Fig. 2(b) we evaluate the clustering results using the cost function. We executed each of the algorithms until it converged and measured the cost of the resulting partition, normalized by the cost of the global minimum. We repeated this process using every possible partition as a starting point to obtain an average cost. For Lloyd's online and Hartigan's algorithms, which allow multiple possible steps in each iteration, we further averaged over all possible choices. The error bars in the plot represent the standard deviation of this result when repeated 100 times for different data samples for each value of m. The graph shows that Hartigan's algorithm has a stable and significant advantage in terms of cost over the online/batch Lloyd's algorithms, and that it consistently converges to a cost close to the global minimum.

In Fig. 2(c) we explored the probability of reaching the global minimum. We define this measure as the probability of reaching a partition which is a global minimum when the initial partition is selected with uniform distribution over all possible partitions and whenever the (online Lloyd's or Hartigan's) algorithm allows multiple choices, a particular one is selected randomly with uniform distribution. The error bars again represent the standard deviation for 100 different data sets for each m. Evidently, even in the most difficult case considered the Hartigan's algorithm has probability larger than 0.5 to reach the global minimum in a single run, while the increasing number of Lloyd's local minima shown in Fig. 2(a), swiftly blocks the batch/online Lloyd's algorithms from reaching it.

Probabilistic Exploration of Larger Synthetic Examples:

We adopted a different simulation strategy to explore additional synthetic examples with larger n and m values. The data samples were drawn as in the previous sub-section, except for the first dimension µ = -5 for the first n/2 data points and µ = 5 for the remaining n/2 data points. We have tested the domain of n = 10, 20, . . . , 400, m = 100, 200, . . . , 4000, K = 2, drawing 40 data sets for each (n,m) combination. In order to cope with these large n values we estimated the number of local minima instead of directly counting them. This was done by randomly choosing 1000 partitions with uniform distribution and counting how many of them correspond to Lloyd's minima and/or to Hartigan's minima. Fig. 3(a) depicts the average number of Lloyd's local minima thus found for every combination of n and m. As can be seen, the result was 1000 out of 1000 for a large part of the domain. Since the 1000 partitions were chosen with repetitions, this can be viewed as 1000 samples of a Bernoulli random variable with parameter p, denoting the frequency of Lloyd's minima. The binomial distribution can be used to calculate a confidence interval for p. For the result of 1000 the 99% confidence interval for p is 0.9947 < p < 1. In conjunction with Theorem 5.4 we conclude that for these values of n and m almost all partitions are Lloyd's local minima. In contrast, the number of Hartigan's local minima was dramatically lower. For n = 10 there were an average of 4.81 Hartigan's local minima found out of 1000 partitions (and specifically 6.5 for m = 4000). For n = 20 the observed average dropped to 0.012 out of 1000 and for n > 20 no further Hartigan's local minima were found.

We then measured the quality of the results of the Hartigan's algorithm in those cases. This was done by executing it 10 times until convergence, starting with random initial partitions. We measured the dependency between the obtained Hartigan's minimum and the true partition of the two Gaussians from which the data was sampled by calculating the normalized mutual information [Cover and Thomas, 2006] between both partitions, where a result of 1.0 means the two are identical. The average information results are depicted in Fig. 3(b). For the most part of this graph the resulting average is a solid 1.0, including a large overlap with the range where the number of Lloyd's local minima is a solid 1000 out of 1000. Thus, while batch/online Lloyd's algorithms utterly fail, the Hartigan's algorithm consistently finds the true solution at once.

Results for real-world data: Finally, we examine the performance of the algorithms on several real world data. These include 9 document subsets of the 20NG corpus and their associated topic labels [Lang, 1995; Slonim et al., 2002]; the SP500 data that include the day-to-day fractional changes in price of the stocks in the Standard and Poor's 500 list during 2003, and their associated Industry Classification (GICS) labels [Slonim et al., 2005]; the EachMovie data that include movie ratings provided by many viewers and the associated genre labels [Slonim et al., 2005], where we focused on movies labeled to a single genre, and the 5000 viewers that rated the maximal number of movies; and the OlivettiFaces data that include 400 images of 40 different people in different poses, available at http://cs.nyu.edu/ roweis/data.html. The details of all datasets are described in Table 1(a). In all cases we used standard pre-processing that included using the tf-idf counts representation for the text data, and further normalizing each input vector to 1 under the Euclidean L2 norm.

For each of these datasets we executed the Lloyd's batch and online algorithms, as well as the Hartigan's algorithm. In all runs we used 500 different random partitions as the initial partitions and executed the three algorithms until they converge to a local minimum. In Table 1(b) we report the average final D(C) values obtained per algorithm. In addition, for the Lloyd's batch and online algorithms we report in parenthesis the - log10(p) score, where p is the p-value of a t-test measuring the statistical significance of the difference between the mean D(C) values associated with Lloyd's batch and online minima to the mean D(C) value associated with the Hartigan's minima, respectively. Evidently, the Hartigan's minima are typically associated with significantly lower D(C) values, especially compared to Lloyd's batch algorithm. Furthermore, at least for the text datasets, the gap monotonically increases with K. To gain further perspective regarding the quality of the obtained minima, in Table 1(c) we report the averaged normalized mutual information between the obtained minima and the true partition of the data, according to the associated objective labels, that were not used during the clustering process. The Hartigan's minima are typically associated with partitions that better capture the data structure as reflected by the objective labels. This is further depicted in Fig. 3(c) and Fig. 3(d), where similar results were obtained for the other datasets.

7 Discussion

In [Banerjee et al., 2005] it was shown that the classical Lloyd's K-means algorithm can be defined and analyzed for any Bregman divergence. Here, we derived Hartigan's Kmeans for any Bregman divergence. Our derivation captures as special cases Hartigan's algorithm with the Euclidean norm [Hartigan, 1975; Telgarsky and Vattani, 2010] as well as the sIB algorithm [Slonim et al., 2002].

We characterize a range of problems for which any random partition will represent a local minimum for Lloyd's algorithm while Hartigan's algorithm easily converges to the true solution. Our experimental results, e.g., Fig. 3(a), provide additional empirical evidence regarding the clear relation between the problem parameters and the associated number of Lloyd's minima and Hartigan's minima, in accord with Theorem 3.1 in [Telgarsky and Vattani, 2010]. In addition, in Proposition 5.5 we provide an upper bound over the number of Hartigan's minima for any given data. Overall, these results suggest that in the presence of many irrelevant features, relying on Hartigan's algorithm might be an attractive alternative to trying to learn the weights of all different dimensions as done in many previous works (e.g., [Huang et al., 2005]).

The theoretical and empirical results of the current work and [Telgarsky and Vattani, 2010] highlight important potential advantages of Hartigan's algorithm over Lloyd's algorithm for K-means clustering. Since the complexity of both algorithms is similar, and since both are equally trivial to implement, one might wonder why is it that Lloyd's algorithm is so prevalent while Hartigan's algorithm is scarcely used in practice. Perhaps it is time for a change.