Hierarchical Clustering

The Story That Explains Hierarchical Clustering

Hierarchical clustering is one of the oldest and most interpretable clustering algorithms in unsupervised machine learning. Unlike K-Means — where you must declare the number of clusters before you begin — hierarchical clustering builds a full nested structure of groupings, letting you choose the depth of detail after you have seen the data.

Two Flavours — Agglomerative vs Divisive

There are two fundamental directions you can build a hierarchy: from the bottom up, or from the top down.

Step-by-Step: How Agglomerative Clustering Works

Linkage Criteria — How to Measure Cluster Distance

Once we merge two points into a cluster, how do we compute the distance from that cluster to all other clusters? This is the linkage criterion and it is arguably the most important design choice in hierarchical clustering.

Distance Metrics — What Is "Similar"?

Before we can measure cluster distances, we need a way to measure point-to-point distances. The choice of metric dramatically changes the result.

The Dendrogram — Reading It Like a Pro

Time & Space Complexity

Hierarchical clustering is elegant but has computational limits. Understanding them helps you decide when to use it vs K-Means or DBSCAN.

Evaluating Cluster Quality

Since hierarchical clustering is unsupervised (no ground truth labels), we need internal validation metrics to judge how good the clusters are.

Python Implementation — Full Walkthrough

Below is a complete, annotated implementation using both scipy (for dendrogram visualisation) and sklearn (for production use).

Step 1 — Import Libraries and Generate Data

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# scipy — dendrogram plotting
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from scipy.spatial.distance   import pdist

# sklearn — clustering + evaluation
from sklearn.cluster           import AgglomerativeClustering
from sklearn.preprocessing     import StandardScaler
from sklearn.metrics           import silhouette_score, davies_bouldin_score
from sklearn.datasets          import make_blobs

# Reproducibility
np.random.seed(42)

# Generate 3-cluster synthetic data
X, y_true = make_blobs.make_blobs(
    n_samples=150,
    centers=3,
    cluster_std=0.8,
    random_state=42
)

print(f"Dataset shape: {X.shape}")

Step 2 — Scale Features

# Always scale before distance-based algorithms
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

print("Mean per feature:  ", X_scaled.mean(axis=0).round(4))
print("Std  per feature:  ", X_scaled.std(axis=0).round(4))

Step 3 — Build Linkage Matrix and Plot Dendrogram

# Build the full linkage matrix using Ward's method
Z = linkage(X_scaled, method='ward', metric='euclidean')

# Z is an (n-1) x 4 matrix; each row:
# [cluster_i, cluster_j, distance, count_in_new_cluster]
print("Linkage matrix shape:", Z.shape)
print("First 3 merges:\n", Z[:3].round(3))

# ── Plot the dendrogram ─────────────────────────────────────
plt.figure(figsize=(14, 6))
plt.title('Hierarchical Clustering Dendrogram (Ward Linkage)', fontsize=14)
plt.xlabel('Sample Index')
plt.ylabel('Distance (Ward)')

dendrogram(
    Z,
    color_threshold=5.0,     # cut at height 5 → colour-codes clusters
    leaf_font_size=8,
    above_threshold_color='grey'
)

# Draw a horizontal cut line at height = 5
plt.axhline(y=5.0, color='red', linestyle='--', label='Cut threshold')
plt.legend()
plt.tight_layout()
plt.show()

Step 4 — Extract Clusters and Evaluate

# ── Method A: scipy fcluster (cut dendrogram at K=3) ───────
labels_scipy = fcluster(Z, t=3, criterion='maxclust')

# ── Method B: sklearn AgglomerativeClustering ──────────────
# Production-grade, integrates with sklearn pipelines
hac = AgglomerativeClustering(
    n_clusters=3,
    linkage='ward',
    metric='euclidean'
)
labels_sklearn = hac.fit_predict(X_scaled)

# ── Evaluation metrics ──────────────────────────────────────
sil  = silhouette_score(X_scaled, labels_sklearn)
dbi  = davies_bouldin_score(X_scaled, labels_sklearn)

print(f"Silhouette Score  : {sil:.4f}  (higher = better, max=1.0)")
print(f"Davies-Bouldin    : {dbi:.4f}  (lower  = better, min=0.0)")
print(f"Cluster sizes     : {np.bincount(labels_sklearn)}")

Step 5 — Silhouette Plot to Find Optimal K

# Cut dendrogram at K = 2 to 8 and score each
k_range = range(2, 9)
sil_scores = []

for k in k_range:
    labels = fcluster(Z, t=k, criterion='maxclust')
    sil_scores.append(silhouette_score(X_scaled, labels))

best_k  = k_range[int(np.argmax(sil_scores))]
best_sc = max(sil_scores)
print(f"Best K = {best_k}  →  Silhouette = {best_sc:.4f}")

for k, sc in zip(k_range, sil_scores):
    bar = '█' * int(sc * 50)
    marker = ' ← BEST' if k == best_k else ''
    print(f"  K={k}  {bar}  {sc:.4f}{marker}")

Hierarchical vs K-Means — When to Use Which

Common Pitfalls and How to Avoid Them

🏪 Customer Segmentation for a Retail Bank

Step 1 — Load & Explore

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics import silhouette_score
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import matplotlib.pyplot as plt

np.random.seed(42)

# Simulated bank customer data (800 customers)
n = 800
df = pd.DataFrame({
    'age':         np.random.randint(22, 70, n),
    'annual_income': np.random.normal(55000, 25000, n).clip(15000, 200000),
    'credit_score':  np.random.randint(450, 850, n),
    'num_products':  np.random.randint(1, 6, n),
    'monthly_txns':  np.random.randint(2, 80, n)
})

print(df.describe().round(1))

Step 2 — Preprocess & Choose Optimal K

# Scale all features to zero mean, unit variance
scaler  = StandardScaler()
X_scaled = scaler.fit_transform(df)

# Build Ward linkage matrix
Z = linkage(X_scaled, method='ward')

# Test K from 2 to 8
results = []
for k in range(2, 9):
    labels = fcluster(Z, t=k, criterion='maxclust')
    sil    = silhouette_score(X_scaled, labels)
    results.append({'k': k, 'silhouette': round(sil, 4)})

best = max(results, key=lambda x: x['silhouette'])
print(f"\n{'K':>4}  {'Silhouette':>12}")
for r in results:
    star = ' ← BEST' if r['k'] == best['k'] else ''
    print(f"  {r['k']:>2}  {r['silhouette']:>12}{star}")

Step 3 — Apply Final Clustering (K = 4) & Profile Segments

# Apply final model with K=4
hac = AgglomerativeClustering(n_clusters=4, linkage='ward')
df['segment'] = hac.fit_predict(X_scaled)

# Profile each segment with mean values
profile = df.groupby('segment').agg(
    count=('age', 'size'),
    avg_age=('age', 'mean'),
    avg_income=('annual_income', 'mean'),
    avg_credit=('credit_score', 'mean'),
    avg_products=('num_products', 'mean'),
    avg_txns=('monthly_txns', 'mean')
).round(1)

print(profile.to_string())

Step 4 — Interpret & Name the Segments

💬 Gene Expression Clustering in Bioinformatics

Step 1 — Simulate Gene Expression Data

import numpy as np
import pandas as pd
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from scipy.spatial.distance   import pdist, squareform
from sklearn.metrics           import silhouette_score
import matplotlib.pyplot       as plt
import seaborn                  as sns

np.random.seed(99)

conditions = ['D0', 'D1', 'D2', 'D4', 'D7', 'D14', 'D21', 'D28']

# Create 3 co-expression groups of ~17 genes each
# Group A: genes that peak early then drop off
base_A = np.array([3.0, 5.2, 7.8, 9.1, 6.4, 3.5, 1.8, 1.2])
# Group B: genes that remain steadily elevated
base_B = np.array([1.0, 1.2, 1.5, 4.8, 6.2, 7.1, 7.4, 7.9])
# Group C: genes that are constitutively suppressed
base_C = np.array([8.5, 7.2, 5.1, 3.3, 2.1, 1.5, 1.0, 0.8])

rows = []
for i in range(17):
    rows.append(base_A + np.random.normal(0, 0.4, 8))
for i in range(17):
    rows.append(base_B + np.random.normal(0, 0.4, 8))
for i in range(16):
    rows.append(base_C + np.random.normal(0, 0.4, 8))

gene_names = [f"GENE_{i+1:02d}" for i in range(50)]
expr = pd.DataFrame(rows, index=gene_names, columns=conditions)

print("Expression matrix shape:", expr.shape)
print(expr.head(3).round(2))

Step 2 — Compute Correlation Distance Matrix

# 1 − Pearson correlation as distance between genes
corr_dist = pdist(expr.values, metric='correlation')
corr_matrix = squareform(corr_dist)

print("Distance matrix shape:", corr_matrix.shape)
print(f"Min dist: {corr_dist.min():.4f}  |  Max dist: {corr_dist.max():.4f}")
print(f"Median dist: {np.median(corr_dist):.4f}")

Step 3 — Cluster Genes with Average Linkage

# Average linkage on precomputed correlation distances
Z_genes = linkage(corr_dist, method='average')

# Find optimal K via silhouette
best_k, best_sil = 2, -1
for k in range(2, 8):
    labels = fcluster(Z_genes, t=k, criterion='maxclust')
    sil    = silhouette_score(corr_matrix, labels, metric='precomputed')
    if sil > best_sil:
        best_k, best_sil = k, sil
    print(f"  K={k}  silhouette={sil:.4f}")

print(f"\nOptimal K = {best_k}  (silhouette = {best_sil:.4f})")

# Apply final clustering
gene_labels = fcluster(Z_genes, t=best_k, criterion='maxclust')
expr['cluster'] = gene_labels

Step 4 — Biological Interpretation

for c in sorted(expr['cluster'].unique()):
    genes_in = expr[expr['cluster'] == c].drop('cluster', axis=1)
    mean_expr = genes_in.mean()
    print(f"\n── Cluster {c}  ({len(genes_in)} genes) ──────────────")
    print("Mean expression per timepoint:")
    for cond, val in mean_expr.items():
        bar = '■' * int(val)
        print(f"  {cond:>4}: {bar:<12} {val:.2f}")
    print(f"  Genes: {', '.join(genes_in.index[:5].tolist())} ...")

sklearn API Quick Reference

from sklearn.cluster import AgglomerativeClustering

# Full parameter reference
model = AgglomerativeClustering(
    n_clusters=4,          # int or None (if distance_threshold used)
    metric='euclidean',    # 'euclidean','manhattan','cosine','precomputed'
    linkage='ward',        # 'ward','complete','average','single'
    distance_threshold=None, # set this instead of n_clusters to cut by distance
    connectivity=None,     # connectivity matrix (e.g. from kneighbors_graph)
    compute_full_tree='auto',# 'auto', True, False
    compute_distances=False  # set True if you need .distances_ attribute
)

labels = model.fit_predict(X_scaled)

# Useful post-fit attributes
print("n_clusters_ : ", model.n_clusters_)    # actual number of clusters
print("n_leaves_   : ", model.n_leaves_)      # number of leaves (= n samples)
print("n_connected_: ", model.n_connected_components_) # connectivity components

Golden Rules — Hierarchical Clustering