08 - K-Means Clustering e clustering gerarchico

Hai una nuvola di punti senza etichette e vuoi sapere se contiene gruppi. Il clustering è il compito di dividere i punti in gruppi in cui i membri si somigliano più di quanto somiglino agli altri, senza che nessuno ti dica prima quali gruppi esistano. Il K-Means lo fa con due mosse ripetute: assegna ogni punto al centro più vicino, poi sposta ogni centro sulla media dei suoi punti.

Converge sempre, ma a un risultato che dipende da dove partono i centri, vuole che tu decida $k$ in anticipo e fallisce in modo sistematico quando i gruppi non sono compatti. Questo articolo, l'ottavo della serie ML Fundamentals, scrive l'algoritmo da zero, misura che cosa cambia con l'inizializzazione, mette a confronto quattro modi di scegliere $k$ e prova K-Means, il clustering gerarchico, DBSCAN e le misture gaussiane su dati dove conosci la risposta giusta.

Per sapere se un clustering è giusto serve una risposta di riferimento, e nella pratica non c'è. Qui la risposta c'è perché i dati sono sintetici (generati con un seme fisso) e le etichette vere si conoscono: per confrontarle con quelle trovate uso l'Adjusted Rand Index (ARI), che vale 1 se le due partizioni coincidono a meno del nome dei gruppi e circa 0 se l'accordo è quello del caso.

Non serve aver letto altro della serie, ma la SVM dell'articolo precedente ti ha già mostrato perché le unità di misura contano: qui contano ancora.

L'algoritmo di Lloyd, scritto da zero

Il K-Means minimizza l'inertia, cioè la somma dei quadrati delle distanze fra ogni punto e il centro (centroide) del suo gruppo: $J = \sum_i \|x_i - \mu_{c(i)}\|^2$. Trovare il minimo assoluto è un problema difficile; l'algoritmo di Lloyd trova un minimo locale alternando due passi.

Nel passo di assegnazione ogni punto va al centroide più vicino, con i centri fissati; nel passo di aggiornamento ogni centroide diventa la media dei punti assegnati, con le assegnazioni fissate. Nessuno dei due passi può far crescere $J$, e siccome le partizioni possibili sono finite, l'algoritmo si ferma: lo afferma anche l'articolo di Arthur e Vassilvitskii, che nota che l'errore totale diminuisce monotonicamente.

Prova su sei numeri sulla retta, 1, 2, 4, 7, 8 e 11, con due centri di partenza scelti male, 1 e 2. Al primo giro il punto 1 va al centro 1 e tutti gli altri al centro 2, che si sposta a 6,4; al secondo giro anche il punto 2 passa al primo gruppo, e così via.

Il listato fa gli stessi conti, stampa $J$ a ogni iterazione e poi confronta la tua funzione lloyd con KMeans di scikit-learn su 600 punti, partendo dagli stessi centri (init con un array e n_init=1), perché solo così i due algoritmi fanno esattamente lo stesso percorso.

# L'algoritmo di Lloyd scritto da zero, provato su 6 numeri e poi confrontato con KMeans
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

def lloyd(X, centri, max_iter=100, verbose=False):
    """Alterna assegnazione ai centri piu' vicini e ricalcolo delle medie."""
    C = np.array(centri, dtype=float)
    lab_prec = None
    for it in range(1, max_iter + 1):
        D = ((X[:, None, :] - C[None, :, :]) ** 2).sum(axis=2)   # distanze^2 punto-centro
        lab = D.argmin(axis=1)
        J = D[np.arange(len(X)), lab].sum()                      # inertia con i centri CORRENTI
        if verbose:
            print(f"iter {it}: etichette={lab.tolist()} centri={np.round(C.ravel(), 3).tolist()} J={J:.3f}")
        if lab_prec is not None and (lab == lab_prec).all():
            return lab, C, J, it
        C = np.array([X[lab == k].mean(axis=0) if (lab == k).any() else C[k] for k in range(len(C))])
        lab_prec = lab
    return lab, C, J, max_iter

# 1) Sei numeri su una retta, due centri iniziali scelti male: 1 e 2
x = np.array([[1.0], [2.0], [4.0], [7.0], [8.0], [11.0]])
lab, C, J, it = lloyd(x, [[1.0], [2.0]], verbose=True)
print("centri finali:", C.ravel().tolist(), " J finale:", J, " iterazioni:", it)

# 2) Stessi centri iniziali su 600 punti 2D: il nostro Lloyd contro sklearn
X, _ = make_blobs(n_samples=600, centers=4, cluster_std=1.1, random_state=3)
C0 = X[:4]
mio_lab, mio_C, mio_J, mio_it = lloyd(X, C0)
km = KMeans(n_clusters=4, init=C0, n_init=1, algorithm="lloyd", tol=0.0, max_iter=100).fit(X)
print("stesse etichette:", bool((mio_lab == km.labels_).all()))
print(f"inertia mia = {mio_J:.6f}  sklearn = {km.inertia_:.6f}")
print("max |centri mio - sklearn| =", f"{np.abs(mio_C - km.cluster_centers_).max():.2e}")

# Output:
# iter 1: etichette=[0, 1, 1, 1, 1, 1] centri=[1.0, 2.0] J=146.000
# iter 2: etichette=[0, 0, 1, 1, 1, 1] centri=[1.0, 6.4] J=30.840
# iter 3: etichette=[0, 0, 0, 1, 1, 1] centri=[1.5, 7.5] J=19.500
# iter 4: etichette=[0, 0, 0, 1, 1, 1] centri=[2.333, 8.667] J=13.333
# centri finali: [2.3333333333333335, 8.666666666666666]  J finale: 13.333333333333332  iterazioni: 4
# stesse etichette: True
# inertia mia = 1454.239467  sklearn = 1454.239467
# max |centri mio - sklearn| = 3.55e-15

L'inertia scende a ogni giro, 146, poi 30,84, poi 19,5, poi 13,333, e l'algoritmo si ferma quando le etichette non cambiano più; i centri finali 2,333 e 8,667 sono le medie di {1, 2, 4} e {7, 8, 11}.

Sui 600 punti le etichette coincidono con quelle di scikit-learn, l'inertia è uguale alla sesta cifra decimale (1454,239467) e i centri differiscono per 3,55e-15: KMeans con algorithm="lloyd" fa esattamente questo, con ottimizzazioni di velocità. L'altra faccia di una convergenza garantita è che il punto di arrivo dipende dal punto di partenza, come mostra la sezione seguente.

L'inizializzazione conta: k-means++ misurato

Il metodo classico sceglie a caso $k$ punti come centri iniziali. L'articolo di Arthur e Vassilvitskii (SODA 2007) propone di scegliere il primo a caso e ciascuno dei successivi con probabilità proporzionale al quadrato della distanza dal centro più vicino già scelto: i centri iniziali tendono così a essere lontani fra loro.

Gli autori dimostrano che l'inertia attesa è al più $O(\log k)$ volte quella ottima e riportano esperimenti in cui k-means++ batte il K-Means classico, spesso con un margine sostanziale. La documentazione di KMeans aggiunge un dettaglio operativo: con init='k-means++' il default n_init='auto' significa una sola partenza, con init='random' dieci.

Per vedere quanto margine c'è ho fatto partire 300 volte, con semi diversi, un K-Means a una sola partenza (n_init=1) con ciascuno dei due metodi su 1.500 punti in due dimensioni che nascono da 10 gruppi veri. Poi ho confrontato la distribuzione delle inertie finali con la migliore trovata nelle 600 partenze, usata come riferimento di «quasi ottimo».

# k-means++ contro inizializzazione casuale: 300 partenze per ciascun metodo, una sola per volta (n_init=1)
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=1500, centers=10, cluster_std=0.9, center_box=(-12, 12), random_state=7)

def prova(init, n=300):
    J, it = [], []
    for seed in range(n):
        km = KMeans(n_clusters=10, init=init, n_init=1, random_state=seed).fit(X)
        J.append(km.inertia_); it.append(km.n_iter_)
    return np.array(J), np.array(it)

J_r, it_r = prova("random")
J_p, it_p = prova("k-means++")
migliore = min(J_r.min(), J_p.min())
print(f"miglior inertia trovata in 600 partenze: {migliore:.2f}")
print(f"{'init':10s} {'media':>9s} {'mediana':>9s} {'massimo':>9s} {'entro 1% dal migliore':>22s} {'iter medie':>11s}")
for nome, J, it in [("random", J_r, it_r), ("k-means++", J_p, it_p)]:
    frac = (J <= migliore * 1.01).mean()
    print(f"{nome:10s} {J.mean():9.1f} {np.median(J):9.1f} {J.max():9.1f} {frac:21.1%} {it.mean():11.1f}")
print(f"inertia media k-means++ rispetto a random: {J_p.mean() / J_r.mean() - 1:+.1%}")

# Output:
# miglior inertia trovata in 600 partenze: 2224.20
# init           media   mediana   massimo  entro 1% dal migliore  iter medie
# random        5530.3    4833.1   19243.9                  6.0%         8.4
# k-means++     2390.9    2435.6    4778.1                 45.7%         5.1
# inertia media k-means++ rispetto a random: -56.8%

L'inertia media con inizializzazione casuale è 5530 e la mediana 4833, contro 2391 e 2436 di k-means++: in media 57% in meno, su questo dataset. Il caso peggiore della partenza casuale arriva a 19.244, più di otto volte il migliore (2224). Conta anche quante partenze finiscono vicino al migliore: 6% per l'inizializzazione casuale, 45,7% per k-means++, che quindi non garantisce l'ottimo, e nemmeno lo trova in una partenza su due.

La cifra del 57% non va generalizzata: dipende dalla separazione dei gruppi e dal loro numero. Il consiglio che regge dopo questa misura è quello della documentazione: k-means++ come default, e n_init=10 invece di 1, tenendo la partenza con inertia minima.

Come scegliere k

Il numero di gruppi non esce dall'algoritmo: lo decidi tu. L'inertia da sola non basta, perché scende sempre quando $k$ cresce (con $k = n$ vale zero). Quattro strumenti usano criteri diversi. Il gomito guarda dove la discesa dell'inertia rallenta. La silhouette di Rousseeuw dà a ogni punto $s = (b - a)/\max(a, b)$, con $a$ distanza media dai punti del proprio gruppo e $b$ distanza media dai punti del gruppo più vicino; va da -1 a 1 e si media sul dataset.

La gap statistic di Tibshirani, Walther e Hastie (2001) confronta il logaritmo dell'inertia osservata con quello atteso su dati di riferimento senza struttura, generati in modo uniforme dentro il riquadro dei dati.

La regola del paper è precisa: si calcola $\mathrm{Gap}(k) = \frac{1}{B}\sum_b \log W^*_{kb} - \log W_k$ su $B$ insiemi di riferimento, con $s_k = \mathrm{sd}_k\sqrt{1 + 1/B}$, e si sceglie il più piccolo $k$ tale che $\mathrm{Gap}(k) \geq \mathrm{Gap}(k+1) - s_{k+1}$. Il listato la implementa con $B = 20$ e la affianca a inertia, silhouette e indice di Davies-Bouldin (più basso è meglio) su quattro gruppi veri, due dei quali sono vicini. Prima verifica a mano la formula della silhouette su otto punti.

# Come scegliere k: inertia, silhouette, Davies-Bouldin e gap statistic (Tibshirani, Walther, Hastie 2001)
import numpy as np
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from sklearn.metrics import davies_bouldin_score, silhouette_samples, silhouette_score

X, _ = make_blobs(n_samples=500, centers=4, cluster_std=1.3, random_state=11)   # 4 gruppi veri, un po' sovrapposti

# --- Controllo a mano della silhouette su 8 punti ---
P = np.array([[0, 0], [0, 1], [1, 0], [5, 5], [5, 6], [6, 5], [3, 3], [9, 9]], float)
lab = np.array([0, 0, 0, 1, 1, 1, 0, 1])
D = np.sqrt(((P[:, None] - P[None]) ** 2).sum(-1))
s_mio = []
for i in range(len(P)):
    stesso = (lab == lab[i]) & (np.arange(len(P)) != i)
    a = D[i, stesso].mean()                      # distanza media nel proprio cluster
    b = D[i, lab != lab[i]].mean()               # distanza media dall'altro cluster
    s_mio.append((b - a) / max(a, b))
print("silhouette a mano = sklearn:", np.allclose(s_mio, silhouette_samples(P, lab)))

# --- Tabella su k = 1..8 ---
rng = np.random.default_rng(0)
lo, hi, B = X.min(axis=0), X.max(axis=0), 20

def W(M, k):
    return KMeans(n_clusters=k, n_init=5, random_state=0).fit(M).inertia_

riferimenti = [rng.uniform(lo, hi, size=X.shape) for _ in range(B)]   # metodo (a) del paper
ks = list(range(1, 9))
gap, s = [], []
for k in ks:
    log_ref = np.log([W(R, k) for R in riferimenti])
    gap.append(log_ref.mean() - np.log(W(X, k)))
    s.append(log_ref.std() * np.sqrt(1 + 1 / B))
gap, s = np.array(gap), np.array(s)

print(f"{'k':>2} {'inertia':>9} {'silhouette':>11} {'Davies-B':>9} {'gap':>7} {'gap(k)-gap(k+1)+s(k+1)':>24}")
for i, k in enumerate(ks):
    inert = W(X, k)
    if k > 1:
        km = KMeans(n_clusters=k, n_init=5, random_state=0).fit(X)
        sil, db = silhouette_score(X, km.labels_), davies_bouldin_score(X, km.labels_)
        ss, dd = f"{sil:11.3f}", f"{db:9.3f}"
    else:
        ss, dd = f"{'-':>11}", f"{'-':>9}"
    crit = f"{gap[i] - gap[i + 1] + s[i + 1]:24.3f}" if k < 8 else f"{'':>24}"
    print(f"{k:2d} {inert:9.1f} {ss} {dd} {gap[i]:7.3f} {crit}")
k_gap = next(k for i, k in enumerate(ks[:-1]) if gap[i] >= gap[i + 1] - s[i + 1])
print("k scelto dalla gap statistic (il piu' piccolo con gap(k) >= gap(k+1) - s(k+1)):", k_gap)

# Output:
# silhouette a mano = sklearn: True
#  k   inertia  silhouette  Davies-B     gap   gap(k)-gap(k+1)+s(k+1)
#  1   21602.0           -         -   0.150                   -0.154
#  2    9250.7       0.540     0.796   0.328                   -0.749
#  3    3010.8       0.675     0.439   1.100                   -0.213
#  4    1649.0       0.636     0.524   1.338                    0.164
#  5    1463.7       0.579     0.664   1.204                    0.128
#  6    1310.4       0.490     0.894   1.102                    0.057
#  7    1174.1       0.450     0.959   1.070                    0.036
#  8    1034.9       0.332     1.123   1.059
# k scelto dalla gap statistic (il piu' piccolo con gap(k) >= gap(k+1) - s(k+1)): 4

La silhouette a mano coincide con quella di scikit-learn. Nella tabella l'inertia ha un gomito fra $k=3$ e $k=4$ (da 3011 a 1649, poi solo 1464), la silhouette è massima a $k=3$ (0,675) e il Davies-Bouldin è minimo a 3 (0,439), mentre la gap statistic sceglie 4, che è il numero vero.

I due centri più vicini distano 4,87 con deviazione standard 1,3: due dei quattro gruppi si toccano, e la silhouette, che premia gruppi ben separati, preferisce fonderli. Non c'è un indicatore che vinca sempre: se i criteri discordano, il dato che hai in più è la domanda di partenza (quanti gruppi puoi usare?) e un grafico dei punti.

Dove K-Means sbaglia

La documentazione dichiara che l'inertia assume gruppi convessi e isotropi (stessa estensione in ogni direzione) e che reagisce male a gruppi allungati o a forme irregolari. Per misurare quanto, ho generato cinque dataset con etichette vere: tre blob sferici; tre strisce orizzontali con i centri sfalsati; tre gruppi di dimensioni molto diverse (500, 50 e 50 punti); due lune; due cerchi concentrici.

Ogni algoritmo riceve il numero giusto di gruppi dove serve ($k$ per K-Means, Ward e single linkage), e per DBSCAN ho preso eps senza guardare le etichette, come 95° percentile della distanza dal quinto punto più vicino (il punto stesso compreso, come in min_samples). Tutti i dati sono standardizzati prima. Il listato ha anche un secondo esperimento: una feature moltiplicata per 1000, come un importo in euro accanto a un conteggio.

# Dove K-Means sbaglia: cinque dataset sintetici con etichette vere note, misura con l'Adjusted Rand Index
import numpy as np
from sklearn.cluster import DBSCAN, HDBSCAN, AgglomerativeClustering, KMeans
from sklearn.datasets import make_blobs, make_circles, make_moons
from sklearn.metrics import adjusted_rand_score as ari
from sklearn.neighbors import NearestNeighbors
from sklearn.preprocessing import StandardScaler

rs = 0
dati = {}
X, y = make_blobs(n_samples=600, centers=3, cluster_std=0.8, random_state=rs); dati["blob sferici"] = (X, y, 3)
def strisce(seed=0, dx=2.0):
    """Tre gruppi allungati in orizzontale (sigma 3 contro 0.39), con i centri sfalsati."""
    rng = np.random.default_rng(seed)
    cov = np.diag([9.0, 0.15])
    medie = [(0, 0), (dx, 1.6), (2 * dx, 3.2)]
    X = np.vstack([rng.multivariate_normal(m, cov, 200) for m in medie])
    return X, np.repeat([0, 1, 2], 200)
X, y = strisce(); dati["tre strisce"] = (X, y, 3)
X, y = make_blobs(n_samples=[500, 50, 50], centers=[[0, 0], [4, 0], [2, 3.5]], cluster_std=[1.3, 0.4, 0.4], random_state=rs)
dati["dimensioni diverse"] = (X, y, 3)
X, y = make_moons(n_samples=500, noise=0.07, random_state=rs); dati["due lune"] = (X, y, 2)
X, y = make_circles(n_samples=500, noise=0.05, factor=0.4, random_state=rs); dati["due cerchi"] = (X, y, 2)

print(f"{'dataset':20s} {'KMeans':>7s} {'Ward':>7s} {'single':>7s} {'DBSCAN':>7s} {'HDBSCAN':>8s}")
for nome, (X, y, k) in dati.items():
    Z = StandardScaler().fit_transform(X)
    # eps senza guardare le etichette: il 95-esimo percentile della distanza dal 5o vicino
    eps = np.quantile(NearestNeighbors(n_neighbors=5).fit(Z).kneighbors(Z)[0][:, -1], 0.95)
    r = [ari(y, KMeans(k, n_init=10, random_state=0).fit_predict(Z)),
         ari(y, AgglomerativeClustering(k, linkage="ward").fit_predict(Z)),
         ari(y, AgglomerativeClustering(k, linkage="single").fit_predict(Z)),
         ari(y, DBSCAN(eps=eps, min_samples=5).fit_predict(Z)),
         ari(y, HDBSCAN(min_cluster_size=20, copy=True).fit_predict(Z))]
    print(f"{nome:20s} " + " ".join(f"{v:7.2f}" for v in r))

# Le unita' di misura: una feature in "euro" (x1000) schiaccia l'altra nella distanza euclidea
X, y = make_blobs(n_samples=600, centers=3, cluster_std=0.8, random_state=rs)
Xu = X.copy(); Xu[:, 1] *= 1000
print("\nblob con la seconda feature moltiplicata per 1000:")
print(f"  KMeans senza scaling: ARI = {ari(y, KMeans(3, n_init=10, random_state=0).fit_predict(Xu)):.2f}")
print(f"  KMeans con StandardScaler: ARI = {ari(y, KMeans(3, n_init=10, random_state=0).fit_predict(StandardScaler().fit_transform(Xu))):.2f}")

# Che cosa fanno single linkage e DBSCAN sui tre blob che si toccano
X, y = make_blobs(n_samples=600, centers=3, cluster_std=0.8, random_state=rs)
Z = StandardScaler().fit_transform(X)
print("\nsingle linkage, dimensioni dei 3 gruppi:", np.bincount(AgglomerativeClustering(3, linkage="single").fit_predict(Z)).tolist())
eps = np.quantile(NearestNeighbors(n_neighbors=5).fit(Z).kneighbors(Z)[0][:, -1], 0.95)
lab = DBSCAN(eps=eps, min_samples=5).fit_predict(Z)
print(f"DBSCAN con eps = {eps:.3f}: cluster trovati = {len(set(lab)) - (1 if -1 in lab else 0)}, punti noise = {(lab == -1).sum()}")

# Output:
# dataset               KMeans    Ward  single  DBSCAN  HDBSCAN
# blob sferici            0.93    0.89    0.00    0.00    0.46
# tre strisce             0.54    0.61   -0.00   -0.00    0.42
# dimensioni diverse      0.29    0.26   -0.01   -0.05    0.62
# due lune                0.48    0.59    1.00    0.98    0.98
# due cerchi             -0.00    0.09    1.00    0.99    1.00
#
# blob con la seconda feature moltiplicata per 1000:
#   KMeans senza scaling: ARI = 0.53
#   KMeans con StandardScaler: ARI = 0.93
#
# single linkage, dimensioni dei 3 gruppi: [597, 1, 2]
# DBSCAN con eps = 0.323: cluster trovati = 1, punti noise = 13

Nessun algoritmo vince su tutte le righe. K-Means e Ward vanno bene sui blob sferici (ARI 0,93 e 0,89) e male altrove: 0,54 sulle strisce, 0,29 sui gruppi di dimensioni diverse, 0,48 sulle lune, -0,00 sui cerchi, dove il K-Means taglia il piano in due semipiani e non può isolare l'anello interno. Il single linkage e il DBSCAN, che seguono la connessione fra punti vicini, risolvono lune e cerchi (ARI 1,00 e 0,98-0,99) ma sui blob sferici valgono zero: i tre blob si toccano, il single linkage stacca solo tre punti isolati (le ultime righe del listato: gruppi di 597, 1 e 2 punti) e il DBSCAN unisce tutto in un solo cluster e scarta 13 punti come rumore. Con densità diverse, sui gruppi di dimensioni diverse, solo HDBSCAN supera 0,6 (0,62).

L'ultimo esperimento riguarda l'unità di misura. Moltiplicare per 1000 una delle due coordinate fa scendere l'ARI del K-Means da 0,93 a 0,53: la distanza euclidea è dominata dalla coordinata grande e l'altra non conta più. Con StandardScaler l'ARI torna a 0,93. Prima di un clustering basato su distanze, la domanda sulle unità di misura va posta sempre.

Il clustering gerarchico

Il clustering agglomerativo parte da un gruppo per punto e fonde a ogni passo i due gruppi più vicini, fino a uno solo. La sequenza delle fusioni, con la distanza a cui avvengono, è un albero detto dendrogramma; tagliarlo a una certa altezza dà una partizione, quindi non devi fissare $k$ prima di vedere l'albero. Cambia il significato di «distanza fra due gruppi», il linkage:

• single: la distanza minima fra un punto dell'uno e uno dell'altro, che concatena punti vicini e segue le forme allungate;
• complete: la distanza massima, che produce gruppi compatti e di diametro limitato;
• average: la media di tutte le distanze fra coppie;
• ward: la fusione che fa crescere di meno l'inertia totale, quindi la più vicina in spirito al K-Means.

Le quattro voci non sono varianti cosmetiche: sullo stesso insieme di punti producono alberi diversi, e la scelta cambia le forme che riesci a trovare (single) o a evitare (complete). La documentazione di scikit-learn presenta Ward come simile alla funzione obiettivo del K-Means, e infatti nella tabella precedente Ward e K-Means ottengono risultati simili sui blob.

Il listato calcola le altezze di fusione dei quattro linkage su cinque numeri, 1, 2, 6, 7 e 15, che si possono ricavare a mano, e controlla che scikit-learn e SciPy diano la stessa partizione.

# Clustering gerarchico: fusioni calcolate a mano su cinque numeri e confrontate con SciPy
import numpy as np
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import pdist
from sklearn.cluster import AgglomerativeClustering
from sklearn.datasets import make_blobs
from sklearn.metrics import adjusted_rand_score

x = np.array([[1.0], [2.0], [6.0], [7.0], [15.0]])
for metodo in ["single", "complete", "average", "ward"]:
    Z = linkage(x, method=metodo)           # righe: [cluster a, cluster b, altezza, dimensione]
    print(f"{metodo:8s} altezze delle 4 fusioni:", np.round(Z[:, 2], 3).tolist())

# Controllo a mano delle altezze:
#   single:   (1,2)=1  (6,7)=1  {1,2}+{6,7}: min = 6-2 = 4   +{15}: min = 15-7 = 8
#   complete: (1,2)=1  (6,7)=1  {1,2}+{6,7}: max = 7-1 = 6   +{15}: max = 15-1 = 14
#   average:  {1,2}+{6,7}: (5+6+4+5)/4 = 5   +{15}: (14+13+9+8)/4 = 11
#   ward:     altezza = sqrt(2 nA nB / (nA + nB)) * |centroide A - centroide B|
for nA, nB, cA, cB in [(2, 2, 1.5, 6.5), (4, 1, 4.0, 15.0)]:
    print("ward a mano:", round(float(np.sqrt(2 * nA * nB / (nA + nB)) * abs(cA - cB)), 3))

# Taglio dell'albero in 2 e in 3 gruppi
Z = linkage(x, method="single")
for k in (2, 3):
    print(f"taglio in {k} gruppi:", fcluster(Z, t=k, criterion="maxclust").tolist())

# sklearn e SciPy fanno la stessa cosa: stesse partizioni su 300 punti
X, _ = make_blobs(n_samples=300, centers=4, cluster_std=1.0, random_state=2)
a = AgglomerativeClustering(n_clusters=4, linkage="ward").fit_predict(X)
b = fcluster(linkage(X, method="ward"), t=4, criterion="maxclust")
print("ARI fra Agglomerative(ward) e linkage+fcluster:", adjusted_rand_score(a, b))

# Il costo in memoria: la matrice delle distanze ha n(n-1)/2 voci da 8 byte
for n in (1_000, 10_000, 100_000):
    voci = n * (n - 1) // 2
    print(f"n = {n:>6}: {voci:>11} distanze, {voci * 8 / 1e6:>7.0f} MB")
print("verifica pdist su 5 punti:", len(pdist(x)), "voci (5*4/2 = 10)")

# Output:
# single   altezze delle 4 fusioni: [1.0, 1.0, 4.0, 8.0]
# complete altezze delle 4 fusioni: [1.0, 1.0, 6.0, 14.0]
# average  altezze delle 4 fusioni: [1.0, 1.0, 5.0, 11.0]
# ward     altezze delle 4 fusioni: [1.0, 1.0, 7.071, 13.914]
# ward a mano: 7.071
# ward a mano: 13.914
# taglio in 2 gruppi: [1, 1, 1, 1, 2]
# taglio in 3 gruppi: [1, 1, 2, 2, 3]
# ARI fra Agglomerative(ward) e linkage+fcluster: 1.0
# n =   1000:      499500 distanze,       4 MB
# n =  10000:    49995000 distanze,     400 MB
# n = 100000:  4999950000 distanze,   40000 MB
# verifica pdist su 5 punti: 10 voci (5*4/2 = 10)

Le altezze coincidono con il conto a mano scritto nei commenti. Con single linkage le fusioni sono (1,2) e (6,7) a distanza 1, poi i due gruppi a distanza 4 (da 2 a 6) e infine il 15 a distanza 8 (da 7 a 15); con complete diventano 6 e 14, perché prendono la distanza massima.

Per Ward vale la formula $\sqrt{2 n_A n_B / (n_A + n_B)}\,\|c_A - c_B\|$, e le altezze 7,071 e 13,914 calcolate a mano nel listato sono uguali a quelle di SciPy. Il costo sta nell'ultima parte dell'uscita: la matrice delle distanze a coppie ha $n(n-1)/2$ voci da 8 byte, 400 MB a diecimila punti e 40.000 MB a centomila. Per questo il clustering gerarchico classico resta uno strumento da qualche migliaio di punti.

DBSCAN e HDBSCAN: la densità

DBSCAN (Ester, Kriegel, Sander e Xu, KDD 1996) non cerca centri, cerca regioni dense. Un punto è un punto core se nel raggio eps ci sono almeno min_samples punti; un cluster è l'insieme dei punti core collegati fra loro, più i punti di bordo vicini a un core; il resto è rumore (etichetta -1). Non serve $k$ e i punti isolati vengono dichiarati rumore invece di essere forzati in un gruppo, ma tutto dipende da eps.

HDBSCAN, disponibile in scikit-learn come sklearn.cluster.HDBSCAN, estende l'idea a più scale di densità e non richiede eps. Il listato fa variare eps su due lune da 500 punti, con min_samples=5, e stampa quanti cluster e quanto rumore risultano; calcola anche i quantili della distanza dal quinto punto più vicino, che sono l'informazione del «k-distance plot» senza il grafico.

# DBSCAN e il raggio eps: sulle due lune il risultato cambia drasticamente con un solo parametro
import numpy as np
from sklearn.cluster import DBSCAN, HDBSCAN
from sklearn.datasets import make_moons
from sklearn.metrics import adjusted_rand_score as ari

X, y = make_moons(n_samples=500, noise=0.07, random_state=0)
print(f"{'eps':>5s} {'cluster':>8s} {'punti noise':>12s} {'ARI':>6s}")
for eps in [0.02, 0.05, 0.1, 0.2, 0.3, 0.5, 1.0]:
    lab = DBSCAN(eps=eps, min_samples=5).fit_predict(X)
    n_cl = len(set(lab)) - (1 if -1 in lab else 0)
    print(f"{eps:5.2f} {n_cl:8d} {(lab == -1).sum():12d} {ari(y, lab):6.2f}")

# Il k-distance plot, senza il grafico: la distanza dal 5o vicino, ordinata
from sklearn.neighbors import NearestNeighbors
d = NearestNeighbors(n_neighbors=5).fit(X).kneighbors(X)[0][:, -1]   # d[:, 0] e' il punto stesso
q = np.quantile(d, [0.5, 0.9, 0.95, 0.99, 1.0])
print("distanza dal 5o vicino, quantili 50/90/95/99/100%:", np.round(q, 3).tolist())

h = HDBSCAN(min_cluster_size=20, copy=True).fit(X)
print(f"HDBSCAN senza eps: {h.labels_.max() + 1} cluster, {(h.labels_ == -1).sum()} noise, ARI = {ari(y, h.labels_):.2f}")

# Output:
#   eps  cluster  punti noise    ARI
#  0.02        0          500   0.00
#  0.05       30          253   0.02
#  0.10        3           18   0.71
#  0.20        2            0   1.00
#  0.30        2            0   1.00
#  0.50        1            0   0.00
#  1.00        1            0   0.00
# distanza dal 5o vicino, quantili 50/90/95/99/100%: [0.062, 0.1, 0.117, 0.148, 0.194]
# HDBSCAN senza eps: 2 cluster, 5 noise, ARI = 0.98

Con eps=0.02 nessun punto ha abbastanza vicini e tutti i 500 sono rumore; con 0.05 il DBSCAN spezza le lune in 30 frammenti; con 0.10 trova 3 cluster e 18 punti di rumore (ARI 0,71); con 0.20 e 0.30 trova le due lune (ARI 1,00); da 0.50 in su le due lune si uniscono in un solo cluster.

L'intervallo buono è largo, ma ha bordi netti. I quantili aiutano a posizionarlo: il 95% dei punti ha il quinto punto più vicino (se stesso compreso) a meno di 0,117, quindi un eps poco sopra quel valore è un buon candidato di partenza, non un valore garantito. HDBSCAN, senza eps, trova 2 cluster, 5 punti di rumore e ARI 0,98. Quando le densità sono simili conviene il DBSCAN; quando sono molto diverse, o non vuoi cercare eps, HDBSCAN.

Le misture gaussiane: gruppi ellittici e appartenenza sfumata

Una mistura gaussiana (Gaussian Mixture Model, GMM) descrive i dati come somma di $k$ distribuzioni normali, ciascuna con la propria media, covarianza e peso. Si addestra con l'algoritmo EM, che ha la stessa struttura a due passi di Lloyd ma assegna a ogni punto una probabilità di appartenenza a ciascun gruppo invece di un'etichetta sola.

Con covarianza 'full' i gruppi possono essere ellissi inclinate, e il numero di componenti si sceglie con il BIC (Bayesian Information Criterion, più basso è meglio), che penalizza i parametri aggiunti. Il listato riprende le tre strisce del confronto precedente, standardizzate, e confronta K-Means e GMM; poi sceglie $k$ con il BIC e mostra le probabilità del punto più ambiguo.

# Gaussian Mixture: cluster ellittici e appartenenze probabilistiche, k scelto con il BIC
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score as ari
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)                          # tre strisce orizzontali con i centri sfalsati
cov = np.diag([9.0, 0.15])
X = np.vstack([rng.multivariate_normal(m, cov, 200) for m in [(0, 0), (2, 1.6), (4, 3.2)]])
y = np.repeat([0, 1, 2], 200)
X = StandardScaler().fit_transform(X)                    # come nel confronto precedente

km = KMeans(3, n_init=10, random_state=0).fit(X)
gm = GaussianMixture(3, covariance_type="full", n_init=5, random_state=0).fit(X)
print(f"ARI K-Means = {ari(y, km.labels_):.2f}   ARI GMM full = {ari(y, gm.predict(X)):.2f}")

bic = [GaussianMixture(k, covariance_type="full", n_init=5, random_state=0).fit(X).bic(X) for k in range(1, 7)]
print("BIC per k = 1..6:", np.round(bic, 1).tolist(), "-> minimo a k =", int(np.argmin(bic)) + 1)

# Il punto piu' ambiguo: probabilita' di appartenenza (soft) contro etichetta unica (hard)
P = gm.predict_proba(X)
i = int(np.argsort(P.max(axis=1))[0])
print("punto piu' incerto:", np.round(X[i], 2).tolist(), " probabilita':", np.round(P[i], 3).tolist(), " etichetta hard:", int(gm.predict(X[[i]])[0]))
print("pesi delle componenti:", np.round(gm.weights_, 3).tolist())

# Output:
# ARI K-Means = 0.54   ARI GMM full = 0.93
# BIC per k = 1..6: [3281.8, 3224.3, 3051.9, 3089.8, 3115.4, 3154.3] -> minimo a k = 3
# punto piu' incerto: [-0.25, -0.58]  probabilita': [0.0, 0.479, 0.521]  etichetta hard: 2
# pesi delle componenti: [0.329, 0.337, 0.334]

Sulle tre strisce il K-Means resta a 0,54 e la mistura sale a 0,93, perché la covarianza completa descrive la forma allungata. Il BIC ha il minimo a $k=3$ (3051,9, contro 3224,3 di $k=2$ e 3089,8 di $k=4$), cioè il numero vero.

Il punto più incerto ha probabilità 0,479 e 0,521 per due componenti: un'etichetta dura lo assegna alla seconda, senza dire che il margine è di quattro punti percentuali. La mistura costa di più, e con pochi punti per componente la covarianza completa può diventare instabile: in quel caso covariance_type='diag' o 'tied' riduce i parametri.

Riepilogo di quando usare ciascun metodo, ricavato dagli esperimenti precedenti e non da una regola generale:

• K-Means o Ward, se i gruppi sono compatti, simili per dimensione, e le feature sono scalate;
• mistura gaussiana, se i gruppi sono ellittici o servono probabilità di appartenenza;
• DBSCAN o HDBSCAN, se i gruppi hanno forma qualsiasi e vuoi scartare i punti isolati come rumore;
• clustering gerarchico, se vuoi vedere l'albero delle fusioni e hai al massimo qualche migliaio di punti.

Nessuna voce è un verdetto: la tabella dell'ARI mostra che con dati diversi ciascun metodo scende a zero. Il criterio che conta è se le ipotesi del metodo (compattezza, forma, densità uniforme) sono ragionevoli per i tuoi dati, e per controllarlo serve guardarli, ad esempio proiettati su due dimensioni con la PCA dell'articolo successivo.

Esercizio svolto

Esegui a mano il K-Means con $k=2$ sui numeri 1, 3, 4, 10, 12, con centri iniziali 1 e 3. Per ogni iterazione scrivi le etichette, i nuovi centri e il valore di $J$. Dopo quante iterazioni ti fermi e quanto vale l'inertia finale? Poi controlla con il listato.

# Esercizio: K-Means (k=2) su 1, 3, 4, 10, 12 con centri iniziali 1 e 3
import numpy as np

x = np.array([1.0, 3.0, 4.0, 10.0, 12.0])
c = np.array([1.0, 3.0])
prec = None
for it in range(1, 10):
    lab = np.abs(x[:, None] - c[None, :]).argmin(axis=1)
    J = ((x - c[lab]) ** 2).sum()
    print(f"iter {it}: centri {np.round(c, 3).tolist()}  etichette {lab.tolist()}  J = {J:.3f}")
    if prec is not None and (lab == prec).all():
        break
    c = np.array([x[lab == k].mean() for k in range(2)])
    prec = lab

# Output:
# iter 1: centri [1.0, 3.0]  etichette [0, 1, 1, 1, 1]  J = 131.000
# iter 2: centri [1.0, 7.25]  etichette [0, 0, 0, 1, 1]  J = 43.125
# iter 3: centri [2.667, 11.0]  etichette [0, 0, 0, 1, 1]  J = 6.667

Che cosa sai fare adesso

Dopo questo articolo sai fare cinque cose, tutte con una verifica numerica a portata di mano:

• scrivere l'algoritmo di Lloyd e dimostrare che l'inertia non cresce mai;
• scegliere l'inizializzazione e n_init sapendo quanto pesano (57% di inertia in meno su un dataset misurato);
• scegliere $k$ incrociando gomito, silhouette, Davies-Bouldin e gap statistic, e spiegare perché possono discordare;
• riconoscere i dati su cui K-Means fallisce e provare gerarchico, DBSCAN o mistura gaussiana;
• controllare le unità di misura prima di qualunque metodo a distanze.

Le cinque capacità ruotano attorno a un'idea: i metodi di clustering fanno ipotesi sulla forma dei gruppi e non te le dicono, quindi sta a te metterle alla prova su dati dove la risposta è nota. Nell'articolo successivo ti serve la stessa disciplina per ridurre le dimensioni: la PCA. Per valutare un modello supervisionato con le metriche giuste c'è invece l'articolo sulle metriche; un altro taglio sul clustering è quello di K-Means e DBSCAN senza etichette, nella serie Machine Learning Engineer. Anche le reti neurali si addestrano migliorando passo dopo passo una funzione di perdita, ma con il gradiente al posto dell'alternanza di Lloyd; e un clustering in uso va monitorato come ogni modello, come mostra l'articolo sul deployment.