Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?

Fionn MurtaghPierre Legendre

article2011Journal of Classification3,532 citations

Resolves widespread discrepancies in how Ward's hierarchical clustering method is formulated and implemented across standard software packages, establishing clear mathematical guidelines for developers and data analysts.

Listen

Two common but subtly different algorithms in statistical software both claim to implement Ward’s hierarchical clustering, which groups observations to minimize the increase in within-cluster sum of squares at each step. The resulting confusion matters because Ward’s method is widely used alongside principal component analysis and K-means partitioning, and incorrect implementations can produce non-optimal or non-comparable cluster trees.

The paper set out to determine which published algorithms and software packages actually optimize Ward’s criterion, to explain why apparently similar code yields different results, and to show users how to obtain correct output.

The authors reviewed the Lance-Williams recurrence formulas that drive most agglomerative clustering routines, derived the precise input and scaling requirements for each variant, and tested the implementations on reproducible simulated data. They compared the R functions hclust and agnes, a custom extension of hclust that offers both variants, and six commercial packages, checking both dendrogram topology and the numerical values of the clustering criterion.

The analysis shows that only one variant—commonly called Ward2—minimizes the Ward criterion when Euclidean distances are supplied directly. The other variant, Ward1, requires squared Euclidean distances as input and produces squared criterion values; taking square roots of its dendrogram heights then recovers identical results. When both algorithms are applied to the same unsquared distance matrix, they generate different trees and different criterion values. Experiments confirmed that agnes (method “ward”) and the “ward.D2” option in an updated hclust implement Ward2 correctly, while the original hclust (method “ward”) implements Ward1. Among commercial packages examined in 2012, Statistica, Systat, and SPSS used Ward1, whereas Matlab, SAS, and JMP used Ward2.

These differences matter for any analysis that relies on the numerical values of the clustering criterion or on consistency with K-means or ordination methods. Using the wrong input or failing to scale output heights can produce trees that do not minimize within-cluster variance and that cannot be compared directly with results from other software.

Users should verify which algorithm their package applies and adjust inputs or post-process outputs accordingly: supply squared distances to Ward1 implementations or use a confirmed Ward2 routine. Software maintainers are advised to document the required steps clearly and, where feasible, to offer both options with appropriate warnings. Further verification on larger or domain-specific data sets would increase practical confidence, but the mathematical identity of the correctly configured algorithms is established.

arXiv: 1111.6285

No sufficiently relevant recommendations were found.

No sufficiently relevant recommendations were found.

Cover for Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?

Abstract

The Ward error sum of squares hierarchical clustering method has been very widely used since its first description by Ward in a 1963 publication. It has also been generalized in various ways. However there are different interpretations in the literature and there are different implementations of the Ward agglomerative algorithm in commonly used software systems, including differing expressions of the agglomerative criterion. Our survey work and case studies will be useful for all those involved in developing software for data analysis using Ward's hierarchical clustering method.

Table of Contents

  • 1. Introduction
  • 2. Applications
  • 3. Ward's Agglomerative Hierarchical Clustering Method
  • 3.1 Some Definitions
  • 3.2 Alternative Expressions for the Variance
  • 3.3 Lance-Williams Dissimilarity Update Formula
  • 3.4 Generalizing Lance-Williams
  • 4. Implementations of Ward’s Method
  • 4.2 Implementation of Ward: Ward1
  • 4.3 Implementation of Ward: Ward2
  • 5. Case Studies: Ward Implementations and Their Relationships
  • 5.1 Experiment 1: Ward2 Implementation in hclust.PL and agnes
  • 5.2 Experiment 2: hclust and Ward1 Implementation in hclust.PL
  • 5.3 Experiment 3: Non-Ward Result Produced by hclust and
  • 5.4 Experiment 4: Modifying Inputs and Options so that Ward1 Output is Identical to Ward2 Output
  • 6. Discussion
  • 6.1 Where Ward1 and Ward2 Implementations Lead to an Identical Result
  • 6.2 How and Why the Ward1 and Ward2 Implementations Can Differ
  • 6.3 Other Software: A Look at Six Other Packages
  • 6.4 Other Implementations Based on the Lance-Williams Update Formula
  • 7. Conclusions
  • References

Knowls

  1. Knowl 1 — Equivalence and Distinction Between Ward1 and Ward2 Hierarchical Clustering Algorithms

    theoretical result

    Two distinct algorithmic formulations exist in hierarchical agglomerative clustering that both claim to implement Ward's minimum variance criterion:

    1. Ward1: Uses the linear Lance-Williams recurrence applied directly to input dissimilarities:

    D(i∪i′,i′′)=∣i∣+∣i′′∣∣i∣+∣i′∣+∣i′′∣D(i,i′′)+∣i′∣+∣i′′∣∣i∣+∣i′∣+∣i′′∣D(i′,i′′)−∣i′′∣∣i∣+∣i′∣+∣i′′∣D(i,i′)D(i \cup i', i'') = \frac{|i| + |i''|}{|i| + |i'| + |i''|} D(i, i'') + \frac{|i'| + |i''|}{|i| + |i'| + |i''|} D(i', i'') - \frac{|i''|}{|i| + |i'| + |i''|} D(i, i')

    where the input matrix must be squared Euclidean distances, D(i,i′)=∥xi−xi′∥2D(i, i') = \|x_i - x_{i'}\|^2.

    1. Ward2: Uses unsquared Euclidean distances as input δ(i,i′)=∥xi−xi′∥\delta(i, i') = \|x_i - x_{i'}\| and updates distances according to:

    δ(i∪i′,i′′)=∣i∣+∣i′′∣∣i∣+∣i′∣+∣i′′∣δ2(i,i′′)+∣i′∣+∣i′′∣∣i∣+∣i′∣+∣i′′∣δ2(i′,i′′)−∣i′′∣∣i∣+∣i′∣+∣i′′∣δ2(i,i′)\delta(i \cup i', i'') = \sqrt{\frac{|i| + |i''|}{|i| + |i'| + |i''|} \delta^2(i, i'') + \frac{|i'| + |i''|}{|i| + |i'| + |i''|} \delta^2(i', i'') - \frac{|i''|}{|i| + |i'| + |i''|} \delta^2(i, i')}

    When Ward1 receives squared Euclidean distances as input and Ward2 receives unsquared Euclidean distances as input, both algorithms optimize the exact same criterion (minimizing the increase in within-cluster error sum of squares) and produce identical dendrogram clustering topologies. They differ solely in the scale of the dendrogram node heights: Ward1 node heights represent squared Euclidean distances, whereas Ward2 node heights represent the square root of Ward1 heights (hWard2=hWard1h_{\text{Ward2}} = \sqrt{h_{\text{Ward1}}}), expressing agglomeration levels on the original distance scale.

  2. Knowl 2 — Ward Minimum Variance Agglomerative Criterion

    equation

    In Ward's hierarchical agglomerative clustering, the pairwise agglomeration chosen at each step minimizes the increase in total within-cluster error sum of squares (the variance index). The objective value resulting from the fusion of two disjoint clusters c1c_1 and c2c_2 is:

    D(c1,c2)=∣c1∣∣c2∣∣c1∣+∣c2∣∥c1−c2∥2D(c_1, c_2) = \frac{|c_1| |c_2|}{|c_1| + |c_2|} \|c_1 - c_2\|^2

    where:

    • ∣c1∣|c_1| and ∣c2∣|c_2| denote the cardinalities (or total masses) of clusters c1c_1 and c2c_2.
    • c1c_1 and c2c_2 inside the norm denote the mean vectors (centroids) of clusters c1c_1 and c2c_2 in multivariate Euclidean space Rp\mathbb{R}^p.
    • ∥c1−c2∥2=∑j=1p(c1j−c2j)2\|c_1 - c_2\|^2 = \sum_{j=1}^p (c_{1j} - c_{2j})^2 is the squared Euclidean distance between the two cluster centroids.

    Minimizing this fusion cost at each step minimizes the increase in within-cluster variance and maximizes the between-cluster variance of the resulting partition.

  3. Knowl 3 — Pitfall of Ward1 Algorithm Applied to Unsquared Euclidean Distances

    empirical result

    When the Ward1 Lance-Williams update formula is executed with an unsquared Euclidean distance matrix as input (as occurred historically in standard calls such as R's hclust(dist(x), method="ward")), the algorithm fails to optimize Ward's minimum variance criterion.

    Consequences of supplying unsquared distances to Ward1 include:

    • The agglomeration process does not minimize the within-cluster error sum of squares.
    • The resulting clustering topology (tree structure and cluster memberships) differs substantially from the correct Ward clustering solution.
    • The dendrogram fusion heights are not monotonically related to true Ward fusion heights.

    To obtain correct Ward clustering when using a Ward1 implementation, input distances must be explicitly squared (e.g., dist(x)^2). To recover fusion levels on the original Euclidean distance scale, the square roots of the resulting node heights must then be extracted.

  4. Knowl 4 — Pairwise Distance Identity for Within-Cluster Sum of Squares

    theoretical result

    For any cluster qq containing ∣q∣|q| observations with centroid q∗=1∣q∣∑i∈qiq^* = \frac{1}{|q|} \sum_{i \in q} i, the within-cluster error sum of squares around the centroid can be calculated directly from all pairwise squared Euclidean distances among cluster members without explicitly computing the centroid:

    ∑i∈qd2(i,q∗)=1∣q∣∑i,i′∈q, i<i′d2(i,i′)\sum_{i \in q} d^2(i, q^*) = \frac{1}{|q|} \sum_{i, i' \in q, \, i < i'} d^2(i, i')

    where:

    • d2(i,q∗)=∥i−q∗∥2=∑j(ij−qj∗)2d^2(i, q^*) = \|i - q^*\|^2 = \sum_j (i_j - q^*_j)^2 is the squared Euclidean distance from observation ii to the centroid q∗q^*.
    • d2(i,i′)=∥i−i′∥2=∑j(ij−ij′)2d^2(i, i') = \|i - i'\|^2 = \sum_j (i_j - i'_j)^2 is the squared Euclidean distance between observations ii and i′i'.

    This identity provides the theoretical bridge allowing agglomerative clustering to minimize variance criteria using only pairwise proximity matrices via Lance-Williams recurrences.

  5. Knowl 5 — Huygens' Variance Decomposition in Agglomerative Clustering

    equation

    By Huygens' theorem, the total variance (inertia) V(I)V(I) of a dataset of observations II with total mass mI=∑i∈Imim_I = \sum_{i \in I} m_i and overall centroid g=1mI∑i∈Imiig = \frac{1}{m_I} \sum_{i \in I} m_i i decomposes additively into between-cluster variance V(Q)V(Q) and summed within-cluster variances for any partition QQ:

    V(I)=V(Q)+∑q∈QmqmIV(q)V(I) = V(Q) + \sum_{q \in Q} \frac{m_q}{m_I} V(q)

    where:

    • V(Q)=∑q∈QmqmI∥q∗−g∥2V(Q) = \sum_{q \in Q} \frac{m_q}{m_I} \|q^* - g\|^2 is the between-cluster variance (dispersion of cluster centroids q∗q^* around the overall mean gg).
    • V(q)=1mq∑i∈qmi∥i−q∗∥2V(q) = \frac{1}{m_q} \sum_{i \in q} m_i \|i - q^*\|^2 is the internal variance of cluster qq.
    • mq=∑i∈qmim_q = \sum_{i \in q} m_i is the mass of cluster qq.

    Because the total variance V(I)V(I) is fixed for a given dataset, maximizing the between-cluster variance V(Q)V(Q) at each agglomerative step is mathematically identical to minimizing the within-cluster variance ∑q∈QmqmIV(q)\sum_{q \in Q} \frac{m_q}{m_I} V(q).

  6. Knowl 6 — Taxonomy of Ward Implementations Across Statistical Software

    empirical result

    Major statistical software packages implement Ward's clustering method using two distinct conventions, leading to potential discrepancies in output topology depending on whether input dissimilarities are squared:

    • Ward1 packages (require squared Euclidean distances to yield true Ward clustering; produce squared-distance fusion heights):

      • Statistica
      • Systat
      • SPSS (implements Ward1 and warns users to supply squared Euclidean distances)
      • R stats::hclust with method="ward" (later named method="ward.D")
    • Ward2 packages (take unsquared Euclidean distances as input, square them internally during update, and produce fusion heights on the original distance scale):

      • Matlab
      • SAS
      • JMP
      • R cluster::agnes with method="ward"
      • R stats::hclust with method="ward.D2"

    Supplying unsquared distances to Ward1 packages produces a non-Ward clustering topology, whereas Ward2 packages correctly minimize the error sum of squares criterion directly from unsquared distances.

  7. Knowl 7 — Reducibility and Quadratic Complexity of Ward's Agglomeration

    theoretical result

    Ward's minimum variance agglomerative clustering criterion satisfies Bruynooghe's reducibility property: the fusion of two clusters ii and jj never decreases the distance to any external cluster kk below min⁡(d(i,k),d(j,k))\min(d(i, k), d(j, k)). This ensures that the agglomeration sequence is monotonic and strictly free of reversals or inversions in dendrogram heights.

    Because reducibility holds, Ward's clustering can be computed using the nearest-neighbor chain algorithm or reciprocal nearest-neighbor algorithm. These algorithms reduce the worst-case computational complexity of the hierarchical agglomeration of nn observations to O(n2)O(n^2), applicable to both stored-dissimilarity implementations and stored-data (on-the-fly distance computation) implementations.

Coverage note — None was omitted; the full methodological comparison, mathematical formulation, software survey, and algorithmic complexity results are captured.

References

  1. 1.ANDERBERG, M.R. (1973), Cluster Analysis for Applications, New York: Academic.
  2. 2.BATAGELJ, V. (1988), “Generalized Ward and Related Clustering Problems”, in Classification and Related Methods of Data Analysis, ed. H.H. Bock, Amsterdam: North-Holland, pp. 67–74.
  3. 3.BENZÉCRI, J.P. (1976), L’Analyse des Données, Tome 1, La Taxinomie (2nd ed.), Paris: Dunod.
  4. 4.BRUYNOOGHE, M. (1977), “Méthodes Nouvelles en Classification Automatique des Données Taxinomiques Nombreuses”, Statistique et Analyse des Données, 3, 24–42.
  5. 5.CAILLIEZ, F., and PAGÈS, J.-P. (1976), Introduction à l’Analyse des Données, SMASH (Société de Mathématiques Appliquées et Sciences Humaines).
  6. 6.FISHER, R.A. (1936), “The Use of Multiple Measurements in Taxonomic Problems”, Annals of Eugenics, 7, 179–188.
  7. 7.GOWER, J.C. (1966), “Some Distance Properties of Latent Root and Vector Methods Used in Multivariate Analysis”, Biometrika, 53, 325–338.
  8. 8.JAIN, A.K., and DUBES, R.C. (1988), Algorithms for Clustering Data, Englewood Cliffs NJ: Prentice-Hall.
  9. 9.JAMBU, M. (1978), Classification Automatique pour l’Analyse des Données. I. Méthodes et Algorithmes, Paris: Dunod.
  10. 10.JAMBU, M. (1989), Exploration Informatique et Statistique des Données, Paris: Dunod.
  11. 11.KAUFMAN, L., and ROUSSEEUW, P.J. (1990), Finding Groups in Data: An Introduction to Cluster Analysis, New York: Wiley.
  12. 12.LANCE, G.N., and WILLIAMS, W.T. (1967), “A General Theory of Classificatory Sorting Strategies. 1. Hierarchical Systems”, Computer Journal, 9(4), 373–380.
  13. 13.LEGENDRE, P., and FORTIN, M.-J. (2010), “Comparison of the Mantel Test and Alternative Approaches for Detecting Complex Relationships in the Spatial Analysis of Data”, Molecular Ecology Resources, 10, 831–844.
  14. 14.LEGENDRE, P. (2011), “const.clust Space- and Time-Constrained Clustering Package”, http://adn.biol.umontreal.ca/~numericalecology/Rcode/
  15. 15.LEGENDRE, P., and LEGENDRE, L. (2012), Numerical Ecology (3rd. English ed.), Amsterdam: Elsevier.
  16. 16.LE ROUX, B., and ROUANET, H. (2004), Geometric Data Analysis: From Correspondence Analysis to Structured Data Analysis, Dordrecht: Kluwer.
  17. 17.MURTAGH, F. (1983), “A Survey of Recent Advances in Hierarchical Clustering Algorithms”, The Computer Journal, 26, 354–359.
  18. 18.MURTAGH, F. (1985), Multidimensional Clustering Algorithms, Vienna: Physica-Verlag.
  19. 19.MURTAGH, F. (1992), “Comments on: Parallel Algorithms for Hierarchical Clustering and Cluster Validity, IEEE Transactions on Pattern Analysis and Machine Intelligence, 14, 1056–1057.
  20. 20.MURTAGH, F. (2000), “Multivariate Data Analysis Software and Resources”, http://www.classification-society.org/csna/mda-sw.
  21. 21.MURTAGH, F. (2005), Correspondence Analysis and Data Coding with R and Java, Boca Raton FL: Chapman & Hall/CRC.
  22. 22.ORLÓCI, L. (1967), “An Agglomerative Method for Classification of Plant Communities, Journal of Ecology, 55, 193–206.
  23. 23.SZÉKELY, G.J., and RIZZO, M.L. (2005), “Hierarchical Clustering Via Joint Between-Within Distances: Extending Ward’s Minimum Variance Method”, Journal of Classification, 22(2), 151–183.
  24. 24.WARD, J.H. (1963), “Hierarchical Grouping to Optimize an Objective Function”, Journal of the American Statistical Association, 58, 236–244.
  25. 25.WISHART, D. (1969), “An Algorithm for Hierachical Classifications”, Biometrics 25, 165–170.
  26. 26.XPLORE (2007), Version 4.8, Collaborative Research Center 649, Humboldt-Universität, Berlin, Germany, http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore.php.

Citation

MLA
Murtagh, F., and P. Legendre. “Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?”. Journal of Classification, vol. 31, no. 3, 2014, pp. 274–95, https://doi.org/10.1007/s00357-014-9161-z.
APA
Murtagh, F., & Legendre, P. (2014). Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?. Journal of Classification, 31(3), 274–295. https://doi.org/10.1007/s00357-014-9161-z
Chicago
Murtagh, F., and P. Legendre. 2014. “Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?”. Journal of Classification 31 (3): 274–95. https://doi.org/10.1007/s00357-014-9161-z.
Harvard
Murtagh, F. and Legendre, P. (2014) “Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?”, Journal of Classification, 31(3), pp. 274–295. Available at: https://doi.org/10.1007/s00357-014-9161-z.
Vancouver
1. Murtagh F, Legendre P (2014) Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?. Journal of Classification 31:274–295

BibTeX

@article{Murtagh_2014, title={Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion?}, volume={31}, ISSN={1432-1343}, url={http://dx.doi.org/10.1007/s00357-014-9161-z}, DOI={10.1007/s00357-014-9161-z}, number={3}, journal={Journal of Classification}, publisher={Springer Science and Business Media LLC}, author={Murtagh, Fionn and Legendre, Pierre}, year={2014}, month=Oct, pages={274–295} }
Metadata:Crossref

Access the Paper

This paper is available from its original source. Click below to access the PDF.

Open PDF