AN EFFECTIVE VARIANT OF THE HARTIGAN K-MEANS ALGORITHM
arXiv:2604.21798v1 [cs.LG] 23 Apr 2026
FRANÇOIS CLÉMENT AND STEFAN STEINERBERGER
Abstract. The k-means problem is perhaps the classical clustering problem and often synonymous with Lloyd’s algorithm (1957). It has become clear that Hartigan’s algorithm (1975) gives better results in almost all cases, TelgarskyVattani note a typical improvement of 5% – 10%. We point out that a very minor variation of Hartigan’s method leads to another 2% – 5% improvement; the improvement tends to become larger when either dimension or k increase.
1. Introduction 1.1. k-means. The k−means problem is as follows: given x1 , . . . , xn ∈ Rd and k ∈ N, the goal is to partition the n points into k clusters S1 , . . . , Sk such that k X X i=1 x∈Si
2
1 X x− xj |Si |
→ min .
j∈Si
P The center of mass |Si |−1 j∈Si xj is also sometimes called the centroid µi ∈ Rd . Geometrically, we are asked to partition the set into k sets that are all as close as possible to their center of mass. This problem has a long history: it was proposed, independently, by Steinhaus [19] in 1956, Lloyd [12] in 1957, Ball and Hall [3] in 1965 and MacQueen [13] in 1967 (see Jain [9]). The problem is NP-hard even for k = 2 clusters [1], the best one can hope for are approximations to the true solution.
=⇒
Figure 1. Thousands points in R2 are clustered into k = 3 sets.
1.2. Lloyd’s algorithm. The most commonly used way to find an approximate solution is the algorithm proposed by Lloyd (and, independently, Forgy [6]), with a variant by MacQueen [13]. For any given cluster partition S1 , . . . , Sk with associated centroids µ1 , . . . , µk , the algorithm goes through every point and assigns xj to the cluster Sℓ whose centroid µℓ is closest to xj (breaking ties in an arbitrary manner). After that, the centroids are updated to account for new cluster assignments 1
2
and the procedure is repeated. The k−means functional can only decrease under one step of this algorithm, the algorithm stops when no point is moved anymore. Since the functional is monotonically decreasing and bounded from below, the algorithm stops eventually. Lloyd’s algorithm is almost trivial to implement, easy to explain and still very widely used today (for example in Python’s sklearn.cluster.KMeans). 1.3. Hartigan’s algorithm. Hartigan’s algorithm [7] was proposed in 1975. It seems to have not been considered much until the early 2010’s when it was again popularized by Telgarsky-Vattani [21], Slonim-Aharoni-Crammer [18] and others. The algorithm has an easy motivation: suppose x is a point currently assigned to cluster 1 and ∥x − µ1 ∥ = ∥x − µ2 ∥. Then Lloyd’s algorithm would be indifferent about moving x, it is already connected to a centroid of closest distance. However, note that if we were to move x over the cluster 2, then µ2 would move in the direction of x since x would then be factored into how µ2 is computed and the k-means functional would decrease.
µ1
µ2
x
Figure 2. A point x assigned to cluster 1 but having the same distance to µ2 that it has from µ1 .
A more formal explanation is as follows (see also §2 for pseudocode): if Si is a cluster, then we abbreviate its contribution to the k−means functional via |Si |
ϕ(Si ) =
X x∈Si
1 X x− y |Si | y∈Si
2
=
X
2
∥x − µ(Si )∥ .
x∈Si
Suppose now that x ∈ Si . Then the decrease in the k−means functional when moving x from Si to Sj is given by ∆(x, Si , Sj ) = ϕ(Si ) + ϕ(Sj ) − ϕ(Si \ {x}) − ϕ(Sj ∪ {x}) and a bit of algebra shows that this can be rewritten as ∆(x, Si , Sj ) =
|Sj | |Si | ∥µ(Sj ) − x∥2 − ∥µ(Si ) − x∥2 . |Sj | + 1 |Si | − 1
Hartigan’s method now picks a point x ∈ Si , checks whether there exists j ̸= i such that ∆(x, Si , Sj ) > 0 and, if so, then it removes x from Si and assigns it to the cluster indexed by arg maxℓ ∆(x, Si , Sℓ ). As pointed out by TelgarskyVattani [21], every local minimum of Hartigan’s algorithm is also a local minimum of Lloyd’s algorithm but not necessarily the other way around: Hartigan’s algorithm may further improve a local Lloyd-minimum. Telgarsky-Vattani [21] summarize an empirical comparison by saying that on average, Hartigan’s method provides an improvement of roughly 5-10% over Lloyd’s method. It is becoming more popular: kmeans() in R calls the Hartigan-Wong implementation [8].
3
2. A variation of the Hartigan algorithm 2.1. The Algorithm. We now present a simple variation of Hartigan’s method.1 We describe a particular formulation of Hartigan’s method (where points are evaluated in the order given by a random permutation which slightly outperforms iid random sampling); our new algorithm is identical except in a single line where the difference is made explicit. Algorithm 1 Hartigan/Smartigan Algorithm. Require: A set {x1 , . . . , xn } ⊂ Rd , number of clusters k, max iterations Niter ∈ N Initialize cluster C1 , . . . , Ck in some way Change=False niter = 0 while Change=False and niter < Niter do Take a random permutation π of {1, . . . , n} Change = True for i = 1 to n do Consider the point xπ(i) , currently associated with cluster Cr . if |Cr | = 1 then Pass else Find cluster Cj (different from Cr ) minimizing ∆j =
|Cj | ||xπ(i) − µ(Cj )||2 |Cj | + 1
if |Cr | ||xπ(i) − µ(Cr )||2 · ∆j ≤ |Cr | − 1
( 1 3 1 niter 2 − 2 Niter
(Hartigan) (Smartigan)
then Assign xπ(i) to Cj Update µ(Cj ) and µ(Cr ) Change = False end if end if end for niter = niter + 1 end while Return sum of squared distances
2.2. Remarks. Several remarks are in order. (1) Hartigan’s method is the obvious thing once one already is close to a good clustering. However, in the beginning, one may only have a very vague notion of the underlying cluster structure: Smartigan encourages exploration. 1The authors kept referring to it, tongue-in-cheek, as Smartigan (because it is both a good idea and also extremely close to Hartigan’s method), slowly got used to the nickname and now cannot bear to part with it.
4
(2) Smartigan, asymptotically, turns into Hartigan; one might specify the algorithm to run Hartigan at the very end which would ensure that one inherits all the guarantees that one has for a Hartigan minimizer. (3) The choice of constants in 3/2 − niter /(2Niter ) is motivated by experiments; the main idea suggests that one should choose a monotonically decreasing function in niter that approaches 1 as niter approaches Niter . Many such functions are conceivable and many seem to lead to good improvements; we picked the linear function for the sake of concreteness, simplicity and performance in practice. (4) When it comes to actual performance, there are relatively few theoretical results in the literature; the difference between Lloyd’s algorithm and Hartigan’s method is seen through numerical experimentation. Moreover, see Telgarsky-Vattani [21], the supremacy of Hartigan’s method is not subtle but very clear and easily observable. Likewise, we will argue that Smartigan outperforms Hartigan in a manner that is equally clear (but smaller in scale than the Lloyd → Hartigan improvement). 2.3. Theoretical guarantees. Very few things are rigorously known for any of these algorithms. It is easy to see that Lloyd-stable assignment, a cluster assignment that remains unchanged under Lloyd’s algorithm, is the weakest form of guarantee: every Hartigan-stable configuration is also Lloyd-stable (since Hartigan is more prone to changing cluster assignments). Similarly, we may deduce that any Smartigan-stable configuration is Hartigan-stable and thus Lloyd-stable. Smartigan-stable
⊆
Hartigan-stable
⊆
Lloyd-stable
We emphasize that Smartigan-stability, the guarantee that cluster assignments remain unchanged independently of how many times the Smartigan algorithm is applied, can be seen as a very powerful form of Hartigan-stability (with an additional safety margin of 50%). One may artificially define Smartigan∗ as Smartigan followed by Hartigan in which case one trivially recovers the guarantees of Hartigan’s algorithm – this, however, would be missing the main point which is the greater exploration of configuration space that occurs early on. 2.4. Un commentaire sociologique. There is a curious discrepancy in the literature that deserves a short sociological comment. The importance of k−means in the literature is beyond doubt; however, there is a gap between how well-known the k−means problem is and by how much Lloyd’s algorithm is reliably and substantially outperformed by Hartigan’s method. This is sometimes, but rarely, hinted at in the literature. Slonim-Aharoni-Crammer mention that 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 [18]. We have no explanation. It is conceivable that the simplicity of Lloyd’s algorithm, it being taught at a basic level, and its easily available implementations give it a distinguished position in the literature that is never questioned. We want to emphasize that the improved performance of both Hartigan and then the further improvement by Smartigan, both easily validated, suggest that it is conceivable that the problem of actually minimizing the k−means functional may have never received the attention it deserves.
5
3. Numerical Results 3.1. Real Data: Low Dimensions. We start with a classic example: the Fisher Iris Data Set comprised of 150 points in R4 describing three different subtypes of the flower Iris, each of them represented 50 times. It is known to be an imperfect example which k−means will not solve with perfect accuracy when k = 3. Algorithm
k=3
k=4
k=5
k = 10
k = 20
Lloyd Hartigan Smartigan
96.88 78.85 79.23
80.82 77.02 57.97 51.37 59.16 49.28
55.61 29.86 28.18
21.68 17.38 16.76
Table 1. Average performance on the Fisher Iris Dataset (random initialization, averaged over 500 runs each) Another reasonably generic data set is Fisher’s cat data set; for each of the 144 cats, the gender, total weight (in kg) and weight of the heart (in g) is recorded. We dropped the gender and worked with the remaining data. The picture is again quite consistent: for k = 2, there seems to be no difference between Hartigan and Smartigan, there is a mild improvement for larger values of k and substantial improvements for k = 10 and k = 20. Algorithm
k=2
k=3
k=5
k = 10
k = 20
Lloyd Hartigan Smartigan
351.05 309.128 309.128
228.472 184.231 182.738
163.07 83.89 83.86
64.37 34.10 31.00
34.49 19.97 16.94
Table 2. Average performance on Fisher’s Cat Dataset (random initialization, averaged over 500 runs each) 3.2. Real Data: High Dimensions. The next example is the Breast Cancer Wisconsin dataset [20] containing 569 points in R30 (30 features of 569 tumors). Algorithm
k=2
k=3
k=5
k = 10
k = 20
Lloyd Hartigan Smartigan
12.93 7.79 7.79
8.08 5.05 4.95
6.56 2.12 2.06
5.42 1.12 1.06
3.46 0.88 0.86
Table 3. Average performance on BCW (random initialization, averaged over 100 runs each; all numbers ×107 ). The final examples come from Lederman et al. [10]. We took the exact implementation and test cases [11], and modified precisely three lines of code: the cluster change condition on which Smartigan differs from Hartigan, and the function definition to include niter . This allows us to replicate perfectly their results from 4 datasets presented in Table 1 in [10]: the Olivetti faces dataset [16], and three datasets from the 20 newsgroups dataset [14]. Table 4 shows that Smartigan gives
6
comparable results for the k-means loss, with usually better NMI values (correlation between the output clustering and the true labels, the closer to 1 the better). Lederman et al. [10] also compared the Hartigan method to the SDP algorithm of Peng-Wei [15] as well as the spectral clustering method of Shi-Malik [17] with Hartigan leading to superior results both in terms of the functional and NMI, we omit these results for the sake of brevity. Dataset Parameters n d K Olivetti 20NG-A 20NG-B 20NG-C
400 200 500 1000
4096 5000 5000 5000
40 2 5 10
k-Means loss Hartigan Smartigan 8.11 193.46 481.72 951.96
8.00 193.46 481.90 953.41
NMI Hartigan Smartigan 0.77 0.54 0.44 0.31
0.78 0.62 0.49 0.25
Table 4. Results obtained by taking the implementation from [10] and changing 3 lines to obtain Smartigan. 3.3. Synthetic Data Sets. It is easy to generate synthetic data. We consider two different types of examples. (1) In the small distance examples, we sample ks centers from [0, 3]d and define them to be Gaussians with covariance matrix 0.3 · Idd×d . (2) In the large distance examples, we sample kl centers from [0, 5]d and consider Gaussians with covariance matrix 0.1 · Idd×d . The small distance problem is more complicated than the large distance problem. Clusters are not guaranteed to have the same number of points, we assign a point to a random cluster before generating the point. d=2
ks = 2
ks = 10
n = 250 > −0.1% −2.3% n = 500 > −0.1% −1.5% n = 1 000 > −0.1% −1.2%
ks = 25
kl = 2
kl = 10
kl = 25
−4.4% > −0.1% −4.2% −4.3% −2.9% > −0.1% −2.7% −3.0% −2.1% > −0.1% −2.7% −2.1%
Table 5. Average difference between Hartigan/Smartigan in 2 dimensions (100 different point sets, k-means++ initialization and averaged over 20 runs each). Negative means Smartigan is better, while > 0.1% indicates the difference is negligible. Any random element of the algorithm is done in the same way for both Hartigan and Smartigan in all these tests: they start with the same point sets and cluster assignments, and the order in which the points are considered is exactly the same. The difference between the two algorithms is already quite noteworthy in two dimensions. When there are only 2 clusters, the problem is easy enough to be solved by either method and we do not see a measurable difference between the two methods (they basically both solve the problem perfectly). However, as soon as there are more clusters, Smartigan leads to consistently better results. The results remain consistent in higher dimensions with the actual improvements becoming larger and larger. In all three cases, dimension d = 2, dimension d = 5
7
d=5
ks = 2
ks = 10
ks = 25
kl = 2
kl = 10
kl = 25
n = 250 < 0.1% −1.7% n = 500 < 0.1% −1.3% n = 1 000 < 0.1% −0.9%
−2.2% −1.6% −1.2%
< 0.1% < 0.1% < 0.1%
−11.8% −9.3% −8.6% −8.4% −8.4% −7.1%
Table 6. Average percentage difference between Hartigan and Smartigan in 5 dimensions (100 different point sets, each with kmeans++ initialization and averaged over 20 runs each). and dimension d = 20, the case of two clusters is solved in a way leading to very comparable scores by both methods; the moment the number of clusters increases, Smartigan gains a definitive advantage. d = 20
ks = 2
ks = 10
ks = 25
kl = 2
kl = 10
kl = 25
n = 250 n = 500 n = 1 000
0% 0% 0%
−6.3% −6.1% −5.2%
−5.5% −4.8% −4.3%
0% 0% 0%
−16.0% −24.1% −12.3% −19.6% −9.8% −17.4%
Table 7. Average percentage difference between Hartigan and Smartigan in 20 dimensions (100 different point sets, each with k-means++ initialization and averaged over 20 runs each). Acknowledgment. We are grateful to Roy Lederman for insightful discussions. References [1] D. Aloise, A. Deshpande, P. Hansen and P. Popat, NP-hardness of Euclidean sum-of-squares clustering. Machine learning, 75 (2009), p. 245-248. [2] D. Arthur and S. Vassilvitskii, k-means++: The advantages of careful seeding, Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms. Society for Industrial and Applied Mathematics Philadelphia, PA, USA. pp. 1027–1035. [3] G. Ball and D. Hall, ISODATA, a novel method of data anlysis and pattern classification. Technical report NTIS AD 699616, 1965, Stanford Research Institute, Stanford, CA. [4] R. A. Fisher, The use of multiple measurements in taxonomic problems, Annals of Eugenics 7 (1936), p. 179–188. [5] R. A. Fisher, The analysis of covariance method for the relation between a part and the whole, Biometrics 3 (1947), 65–68 [6] E. Forgy, Cluster analysis of multivariate data: efficiency versus interpretability of classifications, Biometrics. 21 (1965), 768–769. [7] J. A. Hartigan. Clustering algorithms. Wiley series in probability and mathematical statistics: Applied probability and statistics. Wiley, 1975. [8] H. Hartigan and M. Wong, Algorithm AS 136: A k-Means Clustering Algorithm. Journal of the Royal Statistical Society, Series C. 28 (1979), p. 100-108. [9] A. Jain, Data clustering: 50 years beyond K-means, Pattern Recognition Letters 31 (2010), p. 651–666 [10] R. R. Lederman, D. Silva-Sánchez, Z. Chen, G. Mordant, A. Balanov and T. Bendory, The Catastrophic Failure of the k-means Algorithm in High Dimensions and How Hartigan’s Algorithm Avoids it. ArXiV https://arxiv.org/abs/2602.09936 (2026). [11] R. R. Lederman, https://github.com/Lederman-Group/Catastrophic_Failure_KMeans, accessed April 17, 2026 [12] S. Lloyd, Least squares quantization in PCM. IEEE Trans. Inform. Theory 28 (1982), p. 129–137. Originally as an unpublished Bell laboratories Technical Note (1957).
8
[13] J. MacQueen, Some methods for classification and analysis of multivariate observations. In: Fifth Berkeley Symposium on Mathematics. Statistics and Probability, University of California Press, 1967, pp. 281–297. [14] T. Mitchell, Twenty Newsgroups. UCI Machine Learning Repository 1997. DOI: https://doi.org/10.24432/C5C323. [15] J. Peng and Y. Wei, Approximating k-means-type clustering via semidefinite programming. SIAM journal on optimization, 18(1):186–205, 2007. [16] F. S. Samaria and A. C. Harter, Parameterisation of a stochastic model for human face identification. In Proceedings of 1994 IEEE workshop on applications of computer vision, IEEE (1994), pp 138-142. [17] J. Shi and J. Malik, Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(8):888–905, August 2000. ISSN 1939-3539. doi: 10.1109/34.868688. [18] N. Slonim, E. Aharoni and K. Crammer, Hartigan’s K-means vs. Lloyd’s K means–is it time for a change?. In Proceedings of the 23rd International Joint Conference on Artificial Intelligence (IJCAI), 2013. [19] H. Steinhaus, Sur la division des corps matériels en parties. Bull. Acad. Polon. Sci. 4 (1956): p. 801–804. [20] W. Street, W. Wolberg and O. Mangasarian, Nuclear feature extraction for breast tumor diagnosis. In Biomedical image processing and biomedical visualization, SPIE, 1905 (1993), pp. 861-870. [21] M. Telgarsky and A. Vattani, Hartigan’s method: k-means clustering without voronoi. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings (2010), pp. 820-827. Department of Mathematics, University of Washington, Seattle Email address: [email protected] Department of Mathematics and Department of Applied Mathematics, University of Washington, Seattle Email address: [email protected]