Hierarchical Cluster Analysis

Hierarchical cluster analysis constructs a nested hierarchy of clusters. This page focuses on agglomerative HCA. Divisive hierarchical clustering is a related family. It is not derived here.

Agglomerative HCA begins with each observation as its own cluster and repeatedly merges clusters according to a specified agglomeration criterion. Ward described hierarchical grouping as a successive reduction of $n$ sets to $n-1$ mutually exclusive sets by selecting, at each stage, a union of a pair according to an objective function, continuing until one group remains.

HCA does not use class labels when constructing the hierarchy. Known sample categories may be compared with the resulting tree afterward. They must not be used to construct an unsupervised HCA model. HCA is not a classifier. Unlike LDA, PLS-DA, SVM, or Random Forest, it does not learn from predefined class labels.

Let the data matrix be $X\in {\mathbb{R}}^{n\times p}$, with rows as samples and columns as variables. Initialization places each of the $n$ observations in its own cluster. At every agglomeration step two current clusters are selected by the linkage criterion and merged. For a full binary agglomerative hierarchy that reduction requires $n-1$ merges and ends in a single cluster containing all observations.

The primary output is therefore not one flat partition. It is a nested sequence of merges, equivalently a binary tree. A flat clustering is obtained only after an explicit cut or equivalent stopping rule. Murtagh and Contreras describe that hierarchy as a set of $n-1$ partitions, a dendrogram, and an ultrametric structure on the $n$ objects.

Pairwise dissimilarity $d(x_i,x_j)$ defines what similar means for two observations. The distance metric and the linkage rule are different modeling choices. Not every HCA analysis uses Euclidean distance. Euclidean geometry is required for the Ward formulation taught below and for the Ward implementations in the software examples.

The Euclidean distance used here is

Squared Euclidean distance $d_{\mathrm{E}}(x_i,x_j{)}^2$ omits the square root. Those two quantities are not interchangeable labels. Murtagh and Legendre show that Ward algorithms historically called by the same name can differ according to whether squared distances enter the update.

Distances depend on the numerical representation of the predictors. If variables are expressed on different numerical scales, Euclidean distances can be dominated by variables with larger numerical dispersion. That is a geometric statement. It is not the claim that HCA always requires autoscaling. Scaling and spectroscopic preprocessing are modeling choices determined by the measurement and the scientific question. See Normalization & Scaling. Correlation-based dissimilarities also exist in software, with more than one convention. This page does not adopt a universal correlation distance.

After observation-level distances are defined, agglomerative HCA still needs a rule for the dissimilarity between clusters. That rule is the linkage criterion. Current SciPy and MATLAB documentation use the following three pairwise-linkage definitions.

Single linkage takes the nearest cross-cluster pair:

Complete linkage takes the farthest cross-cluster pair:

Average linkage on this page is the unweighted pair-group arithmetic mean (UPGMA):

That is SciPy method="average" and MATLAB 'average'. It is not WPGMA, centroid, or median linkage. Different linkage rules can produce different hierarchies from the same distances. scikit-learn's clustering guide notes that agglomerative merging can produce uneven cluster sizes, and in that comparison describes single linkage as the most uneven of the four strategies it discusses. That is a software-guide comparison. It is not a law that single linkage always forms chains.

Ward linkage is not a min, max, or mean of pairwise distances. Under the Euclidean error-sum-of-squares formulation, each cluster $C$ has centroid ${\mu }_C$ and within-cluster sum of squares

Merging clusters $A$ and $B$ increases that total by

The selected pair is

Murtagh and Legendre write the same increment as

That quantity is the merge cost in the within-cluster sum of squares. It is not automatically the number stored as dendrogram height in every software package. Murtagh and Legendre distinguish Ward1 and Ward2 implementations that can differ when the same unsquared distance matrix is supplied, and they report that MATLAB implements Ward2. Current SciPy Ward updates use a square-root recurrence. Current MATLAB documents a factor chosen so that the Ward distance between two singleton clusters equals the Euclidean distance. Those dendrogram heights are therefore not claimed to equal $\Delta (A,B)$.

Current scikit-learn AgglomerativeClustering with linkage="ward" accepts only metric="euclidean". SciPy states that Ward is correctly defined only for a Euclidean pairwise metric. Do not combine Ward with cosine or Manhattan in the sklearn examples on this page.

HCA does not have one universal equation. It is a distance definition, a linkage or agglomeration criterion, and iterative merging.

The procedure below is the generic agglomerative algorithm. Ward-specific merge costs are not inserted into the generic loop. When the selected linkage is Ward, the pair that is merged is the pair that minimizes the increase in the within-cluster sum of squares under the verified Euclidean implementation. That is not the minimum Euclidean distance between cluster members.

A dendrogram is a graphical representation of the hierarchical merge structure. Leaves correspond to individual observations in the usual sample-clustering display. Internal joins correspond to merges. The axis that displays linkage height depends on plot orientation. Dendrogram height is the stored linkage value at that merge. It is not always raw Euclidean distance. It depends on the linkage criterion and on the software convention, especially for Ward.

The left-to-right order of leaves is not a unique scientific ordering. Multiple leaf permutations can represent the same hierarchy by swiveling subtrees. Adjacent leaves in a plot are not, by themselves, proof of closest similarity. SciPy exposes an optional optimal_ordering argument that reorders the linkage matrix so that the distance between successive leaves is minimized. That option is not used in the code below and is not required to define the hierarchy.

A hierarchical tree can be cut to obtain a flat partition. SciPy fcluster can request a maximum number of clusters or a cophenetic-distance threshold. scikit-learn AgglomerativeClustering can specify n_clusters or, alternatively, distance_threshold with n_clusters=None. MATLAB cluster can request a maximum number of clusters. There is no universal scientifically correct number of clusters determined by HCA alone. The dendrogram does not reveal a uniquely correct $k$.

Selecting or validating a partition is a separate problem from constructing the hierarchy. Internal cluster-validity indices are not derived on this page. If known chemical classes are available, they may be inspected after clustering as external information. They must not be used to tune an analysis that is still described as unsupervised.

HCA is used exploratorily with spectral profiles and other multivariate analytical measurements. Yu applied hierarchical cluster analysis together with PCA to synchrotron FTIR microspectra of feed tissues and reported that both methods could discriminate inherent structures and molecular-chemistry differences in that study. That is an application result. It does not make HCA inherently spectroscopic, and it does not prove that dendrogram branches are chemically meaningful classes in general.

Spectroscopic preprocessing such as SNV, MSC, derivatives, baseline correction, normalization, or scaling changes the geometry supplied to the distance and linkage rules. PCA constructs continuous latent directions. HCA constructs a hierarchical grouping. Neither substitutes for the other. Murtagh and Legendre note that Ward clustering and PCA share a Euclidean variance geometry and are complementary rather than equivalent. PCA scores may be used as the representation passed to HCA. That is optional. The number of retained components and any scaling change the distances. PCA is not a required preprocessor.

This page clusters samples, that is, rows of $X$. HCA can also cluster variables if the objects represented by rows are variables. Do not silently transpose $X$. High-dimensional, highly correlated spectral channels can change distance structure. HCA does not automatically solve that geometry. A cluster tree can reflect chemical structure, scale, preprocessing, batch or instrument effects, noise, metric, linkage, or combinations of these.

Relative to k-means, HCA builds a hierarchy by sequential merging. k-means seeks a flat partition of specified size according to its own objective. k-means is not derived here. Ordinary classification accuracy is not an HCA performance metric unless labels are used only as post-hoc external information under an explicitly stated question. HCA does not, by itself, assign labels to new observations in the sense of a fitted supervised classifier. An additional assignment rule would be required.

The first Python listing implements Euclidean distances, the three pairwise linkages, the Ward merge cost, and a small agglomerative loop for single, complete, and average linkage. It does not implement optimized Ward. Production hierarchies use SciPy linkage. Flat partitions use fcluster or sklearn AgglomerativeClustering. pdist returns a condensed distance vector. Pass that vector, or the observation matrix, to linkage. Do not pass a square distance matrix. The sklearn parameter is metric, not deprecated affinity. No scaler is applied automatically. Cluster integer labels from SciPy and sklearn may be permutations of the same partition.

MATLAB examples use official pdist, linkage, and cluster. Rows are observations. dendrogram is the MathWorks plotting function for the tree. It is not called in the public code. SciPy, sklearn, and MATLAB Ward hierarchies are not claimed to be numerically identical in stored merge heights. Compare partitions after a matched cut, allowing for label permutation, and only under matched Euclidean Ward conventions.

import numpy as np  def _validate_matrix(X):    X = np.asarray(X, dtype=float)    if X.ndim != 2:        raise ValueError("X must be a 2D array of samples by variables.")    if min(X.shape) == 0:        raise ValueError("X must have at least one row and one column.")    if not np.all(np.isfinite(X)):        raise ValueError("X must contain only finite values.")    return X  def euclidean_distance(x_i, x_j):    x_i = np.asarray(x_i, dtype=float).reshape(-1)    x_j = np.asarray(x_j, dtype=float).reshape(-1)    if x_i.shape[0] != x_j.shape[0]:        raise ValueError("x_i and x_j must have the same number of variables.")    return float(np.sqrt(np.sum((x_i - x_j) ** 2)))  def pairwise_euclidean(X):    X = _validate_matrix(X)    n = X.shape[0]    D = np.zeros((n, n), dtype=float)    for i in range(n):        for j in range(i + 1, n):            d = euclidean_distance(X[i], X[j])            D[i, j] = d            D[j, i] = d    return D  def _pair_distances(D, members_a, members_b):    D = np.asarray(D, dtype=float)    a = np.asarray(list(members_a), dtype=int)    b = np.asarray(list(members_b), dtype=int)    if a.size == 0 or b.size == 0:        raise ValueError("Cluster member sets must be non-empty.")    return D[np.ix_(a, b)]  def single_linkage(D, members_a, members_b):    return float(np.min(_pair_distances(D, members_a, members_b)))  def complete_linkage(D, members_a, members_b):    return float(np.max(_pair_distances(D, members_a, members_b)))  def average_linkage(D, members_a, members_b):    pair = _pair_distances(D, members_a, members_b)    return float(np.sum(pair) / (pair.shape[0] * pair.shape[1]))  def cluster_centroid(X, members):    X = _validate_matrix(X)    idx = np.asarray(list(members), dtype=int)    if idx.size == 0:        raise ValueError("Cluster member sets must be non-empty.")    return X[idx].mean(axis=0)  def within_cluster_sse(X, members):    X = _validate_matrix(X)    idx = np.asarray(list(members), dtype=int)    centroid = cluster_centroid(X, idx)    return float(np.sum((X[idx] - centroid) ** 2))  def ward_merge_cost(X, members_a, members_b):    X = _validate_matrix(X)    a = np.asarray(list(members_a), dtype=int)    b = np.asarray(list(members_b), dtype=int)    w_ab = within_cluster_sse(X, np.concatenate([a, b]))    return w_ab - within_cluster_sse(X, a) - within_cluster_sse(X, b)  def ward_merge_cost_closed_form(X, members_a, members_b):    X = _validate_matrix(X)    mu_a = cluster_centroid(X, members_a)    mu_b = cluster_centroid(X, members_b)    n_a = float(len(list(members_a)))    n_b = float(len(list(members_b)))    return (n_a * n_b / (n_a + n_b)) * float(np.sum((mu_a - mu_b) ** 2))  def agglomerative_merges(X, method):    method = str(method)    if method not in {"single", "complete", "average"}:        raise ValueError("method must be single, complete, or average.")    X = _validate_matrix(X)    n = X.shape[0]    if n < 2:        raise ValueError("Agglomerative clustering requires at least two observations.")    D = pairwise_euclidean(X)    members = {i: {i} for i in range(n)}    active = list(range(n))    merges = []    next_id = n    while len(active) > 1:        best = None        for ia, a in enumerate(active):            for b in active[ia + 1 :]:                if method == "single":                    d = single_linkage(D, members[a], members[b])                elif method == "complete":                    d = complete_linkage(D, members[a], members[b])                else:                    d = average_linkage(D, members[a], members[b])                key = (d, min(a, b), max(a, b))                if best is None or key < best[0]:                    best = (key, a, b, d)        _, a, b, d = best        new_id = next_id        members[new_id] = set(members[a]) | set(members[b])        merges.append((a, b, d, len(members[new_id])))        active = [c for c in active if c not in (a, b)] + [new_id]        next_id += 1    return merges 

1. 1.

   Ward, J. H., Jr. (1963). Hierarchical Grouping to Optimize an Objective Function. Journal of the American Statistical Association, 58(301), 236-244.

   doi:10.1080/01621459.1963.10500845
2. 2.

   Murtagh, F., & Legendre, P. (2014). Ward's Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward's Criterion? Journal of Classification, 31, 274-295.

   doi:10.1007/s00357-014-9161-z
3. 3.

   Murtagh, F., & Contreras, P. (2011). Methods of Hierarchical Clustering. arXiv:1105.0121.

4. 4.

   Yu, P. (2005). Applications of Hierarchical Cluster Analysis (CLA) and Principal Component Analysis (PCA) in Feed Structure and Feed Molecular Chemistry Research, Using Synchrotron-Based Fourier Transform Infrared (FTIR) Microspectroscopy. Journal of Agricultural and Food Chemistry, 53(18), 7115-7127.

   doi:10.1021/jf050959b
5. 5.

   SciPy Developers (n.d.). scipy.cluster.hierarchy.linkage, fcluster, and scipy.spatial.distance.pdist. SciPy 1.17 documentation.

6. 6.

   scikit-learn Developers (n.d.). AgglomerativeClustering and the hierarchical clustering user guide. scikit-learn 1.6 documentation.

7. 7.

   The MathWorks, Inc. (n.d.). pdist, linkage, and cluster. MATLAB documentation.