Clustering e k-means

Il clustering è il secondo grande compito non supervisionato, dopo il rilevamento di anomalie (Anomaly detectionUn'anomalia (outlier) è un'osservazione che si discosta tanto dalle altre da far pensare che sia generata da un meccanismo diverso. Rilevarle serve come pulizia dei dati (solo se sono errori o rumore, non per migliorare artificialmente le metriche), come obiettivo finale (frodi, guasti, cybersicurezza) e per il monitoraggio di un modello in produzione. Metodi semplici: box plot (oltre $1{,}5,\mathrm{IQR}$), carte di controllo univariate ($\mu\pm3\sigma$) e la statistica multivariata di Hotelling $T^2=(x-\bar x)^TS^{-1}(x-\bar x)$ con soglia $\chi^2_{p,1-\alpha}$, valida per dati gaussiani e unimodali. Metodi non supervisionati multivariati danno un anomaly score: l'isolation forest isola ogni punto con split casuali (le anomalie hanno cammini corti) e calcola $s(x,n)=2^{-E(h(x))/c(n)}$, con soglia scelta dalla contaminazione. Senza etichette si valuta con esperti, eventi noti o anomalie sintetiche. Approfondimento: non nel programma di Telecomunicazioni.Anomaly detection →): raggruppare i dati in modo che i membri di uno stesso gruppo siano più simili tra loro che a quelli degli altri, senza usare etichette né categorie predefinite (Lezione 20 · Clustering e k-means, Lezione 21 · Laboratorio anomaly detection e clustering). Gli esercizi sono Esercizio - k-means a mano e gap statistic su cinque gruppi, Esercizio - Clustering gerarchico agglomerativo e divisivo, Esercizio - Regressione lineare a tratti con k-means e Esercizio - k-means online e con buffer a memoria limitata.

A che cosa serve

• Come preprocessing. Quando i dati sono eterogenei o complessi, si può usare un modello per ogni cluster (modelli locali): può migliorare molto le prestazioni (per esempio un reattore chimico che lavora in tre punti di lavoro diversi, ciascuno con una relazione lineare propria: Esercizio - Regressione lineare a tratti con k-means). Svantaggi: più complessità, e meno dati per ogni modello, con possibile peggioramento se i dati sono pochi.
• Come obiettivo finale. Marketing e pubblicità digitale (segmentazione dei clienti), organizzazione di documenti (topic modelling), analisi degli utenti di un sito, produzione industriale (dopo il clustering si può installare uno strumento di classificazione).

I due approcci visti sono il k-means e il clustering gerarchico.

K-means

Definizione (k-means). Se si conoscessero i centroidi (i centri) di $K$ gruppi, basterebbe assegnare ogni punto al centroide più vicino. L'algoritmo li trova insieme alle assegnazioni:

1. Scegliere $K$ (iperparametro): il numero di cluster.
2. Inizializzare i centroidi: si scelgono a caso $K$ dei punti dati come centri iniziali.
3. Assegnare ogni punto al centroide più vicino (distanza euclidea o altra).
4. Aggiornare ogni centroide alla media dei punti assegnati al suo cluster.
5. Ripetere 3-4 fino a convergenza (i centroidi non cambiano più, o cambiano meno di una tolleranza, o si è raggiunto il numero massimo di iterazioni).

centroidi = K punti casuali del dataset
ripeti:
    etichetta[i] = indice del centroide più vicino a x_i
    nuovo_centroide[k] = media dei punti con etichetta k
    se i centroidi non cambiano (entro tol): stop

Che cosa minimizza

Per definire la «convergenza» serve una metrica: l'errore quadratico medio entro i cluster, $$\mathrm{MSE}_{\text{within}}=\frac1N\sum_{k=1}^K\sum_{x_i\in C_k}\lVert x_i-\mu_k\rVert^2,$$ dove $N$ è il numero totale di punti, $K$ il numero di cluster, $C_k$ l'insieme dei punti del cluster $k$ e $\mu_k$ il suo centroide (la media). Si usa anche la dispersione entro i cluster non divisa per $N$: $W_K=\sum_k\sum_{x\in C_k}\lVert x-\mu_k\rVert^2$ (inerzia).

Perché l'algoritmo funziona (non sale mai). Ciascuno dei due passi fa diminuire (o lascia uguale) $W$:

• Assegnazione: spostare un punto dal suo cluster a quello con centroide più vicino non può aumentare la somma delle distanze al quadrato.
• Aggiornamento: fissate le assegnazioni, per un cluster il punto $\mu$ che minimizza $\sum_{x\in C}\lVert x-\mu\rVert^2$ è la media: si annulla il gradiente rispetto a $\mu$, $\sum_{x\in C}2(\mu-x)=0\Rightarrow|C|\,\mu=\sum_{x\in C}x\Rightarrow\mu=\frac1{|C|}\sum_{x\in C}x$ (Derivate parziali e derivate direzionaliLa derivata parziale ∂f/∂xi(x0) è la derivata della funzione di una variabile che si ottiene fissando tutte le altre variabili: in pratica si deriva in xi trattando le altre come costanti. Il gradiente ∇f raccoglie le derivate parziali. La derivata direzionale lungo un versore v è il limite di [f(x0 + hv) − f(x0)]/h. In più variabili avere tutte le derivate (anche direzionali) non garantisce la continuità.Derivate parziali e derivate direzionali →, Massimi e minimi liberi in più variabiliSe f è differenziabile e ha un estremo relativo in un punto interno x0, allora ∇f(x0) = 0 (Fermat): gli estremi interni vanno cercati tra i punti critici. Se f è C², in un punto critico: hessiana definita positiva → minimo, definita negativa → massimo, indefinita → sella, semidefinita → il test non decide e si studia il segno di f(x) − f(x0). In due variabili: det H > 0 e fxx > 0 minimo, det H > 0 e fxx < 0 massimo, det H < 0 sella.Massimi e minimi liberi in più variabili →; analogo statistico: la media campionaria, Valore attesoIl valore atteso E[X] = Σ x p_X(x) è la media dei valori di X pesata con le loro probabilità (esiste se la serie converge assolutamente); per una funzione g vale E[g(X)] = Σ g(x) p_X(x) senza trovare la legge di g(X), ed E è lineare: E[aX + bY + c] = aE[X] + bE[Y] + c.Valore atteso →). Si tratta di un minimo perché la funzione è una paraboloide. Poiché $W$ non cresce e il numero di partizioni dei punti è finito, l'algoritmo converge. Ma a un minimo locale, non necessariamente globale.

Esempio svolto. Sette punti: A$(1,1)$, B$(1{,}5;2)$, C$(3,4)$, D$(5,7)$, E$(3{,}5;5)$, F$(4{,}5;5)$, G$(3{,}5;4{,}5)$ con $K=2$ e centroidi iniziali A e D.

• Assegnazione 1. Distanze da $c_1=(1,1)$ e $c_2=(5,7)$: A $(0;\,7{,}21)\to1$; B $(1{,}12;\,6{,}10)\to1$; C $(3{,}61;\,3{,}61)$ pareggio, assegnato al primo $\to1$; D $(7{,}21;\,0)\to2$; E $(4{,}72;\,2{,}50)\to2$; F $(5{,}32;\,2{,}06)\to2$; G $(4{,}30;\,2{,}92)\to2$.
• Aggiornamento 1. $c_1=\text{media}(A,B,C)=(\frac{1+1{,}5+3}3,\frac{1+2+4}3)=(1{,}833;\,2{,}333)$; $c_2=\text{media}(D,E,F,G)=(\frac{5+3{,}5+4{,}5+3{,}5}4,\frac{7+5+5+4{,}5}4)=(4{,}125;\,5{,}375)$. $W=12{,}21$.
• Assegnazione 2. C è ora più vicino a $c_2$ (distanza da $c_1$: $\sqrt{1{,}167^2+1{,}667^2}=2{,}03$; da $c_2$: $\sqrt{1{,}125^2+1{,}375^2}=1{,}78$). Nuovi cluster: $\{A,B\}$ e $\{C,D,E,F,G\}$. Centroidi $(1{,}25;\,1{,}5)$ e $(3{,}9;\,5{,}1)$, $W=8{,}525$.
• Assegnazione 3: nessun cambiamento, convergenza.

Sensibilità all'inizializzazione. I centroidi iniziali casuali possono portare a minimi locali diversi. Esempio: i quattro punti $(0,0),(0,1),(4,0),(4,1)$ con $K=2$. Se i centroidi iniziali sono $(0,0)$ e $(0,1)$, il punto $(4,0)$ è più vicino al primo e $(4,1)$ al secondo: i cluster sono $\{(0,0),(4,0)\}$ e $\{(0,1),(4,1)\}$ con centroidi $(2;0)$ e $(2;1)$, e l'algoritmo è già fermo con $W=4\cdot4=16$. La soluzione migliore (sinistra e destra) ha centroidi $(0;0{,}5)$ e $(4;0{,}5)$ e $W=4\cdot0{,}25=1$. Per questo si fa partire l'algoritmo più volte con inizializzazioni diverse e si tiene la migliore (o si fissa un seme); esiste anche l'inizializzazione k-means++, che sceglie centri iniziali lontani tra loro.

Ipotesi implicite. Il k-means separa bene gruppi compatti e di forma sferica, con dimensioni simili; funziona male per gruppi allungati o a forma di anelli (vedi il gerarchico con linkage single). Come ogni metodo con distanze, richiede feature standardizzate (Classificazione e k-nearest neighborsNella classificazione l'uscita $y$ è una categoria (con $C$ classi; $C=2$ è il caso binario). Il classificatore più semplice è il k-nearest neighbors: una nuova osservazione prende la classe più frequente (voto di maggioranza) tra i suoi $k$ vicini più prossimi nel training, con distanza euclidea $\sqrt{\sum(A_i-B_i)^2}$ o di Manhattan $\sum|A_i-B_i|$ (per la regressione si fa la media dei vicini). $k$ è un iperparametro: $k$ piccolo dà bordi frastagliati e overfitting, $k$ grande underfitting. È un metodo basato su istanze e «pigro» (nessun addestramento, costo alla predizione), sensibile a scala e feature irrilevanti e alla maledizione della dimensionalità; gli ingressi categorici si codificano con one-hot. Approfondimento: non nel programma di Telecomunicazioni.Classificazione e k-nearest neighbors →, Statistica per il machine learningI dati di un problema ML si organizzano nella matrice di progetto $X$ ($n$ osservazioni, $p$ variabili). La statistica serve a capirli, ripulirli e prepararli: i momenti (media $\mu$, varianza $\sigma^2$, asimmetria, curtosi), i quartili con lo scarto interquartile $\mathrm{IQR}=Q_3-Q_1$ (all'esame senza interpolazione), la moda per i dati categorici. Con queste quantità si imputano i dati mancanti (media o mediana), si eliminano le variabili costanti e si standardizza con lo z-score $z=(x-\mu)/\sigma$, usando sempre media e deviazione standard del solo training set.Statistica per il machine learning →). Il k-means è anche usabile per classificare un nuovo punto: lo si assegna al centroide più vicino.

Versioni a memoria limitata

Il k-means standard richiede di conservare tutti i punti per ricalcolare le medie. Con la media incrementale si può aggiornare un centroide senza i punti: se il cluster ha $n$ punti e centroide $\mu$ e arriva un nuovo punto $x$ assegnato a quel cluster, $$\mu_{\text{new}}=\frac{n\,\mu+x}{n+1}=\mu+\frac{x-\mu}{n+1},\qquad n\leftarrow n+1,$$ (la seconda forma si ottiene scrivendo $n\mu+x=(n+1)\mu+(x-\mu)$). Basta quindi tenere in memoria solo i $K$ centroidi e i $K$ contatori $n_k$ (k-means online). Una variante intermedia, il k-means a buffer, tiene al massimo $M=m\times K$ punti campionati e ricalcola i centroidi solo da quelli. Esempio: centroide $(1;1)$ con $n=3$ punti e nuovo punto $(5;5)$: $\mu_{\text{new}}=\frac{3(1;1)+(5;5)}4=(2;2)$.

Quanti cluster?

Il numero $K$ non è noto a priori.

Metodo del gomito

Al crescere di $K$ la dispersione entro i cluster diminuisce sempre (più cluster adattano meglio i dati: per $K=N$ è zero), ma dopo un certo punto il guadagno è marginale e i nuovi cluster si limitano a sovradattare. Il «gomito» del grafico $W_K$ (o dell'errore medio) contro $K$ è il punto in cui il miglioramento comincia ad appiattirsi.

Esempio (laboratorio: 1000 punti in cinque gruppi gaussiani, deviazione standard $1{,}5$, centri $(-10,-15)$, $(0,0)$, $(10,-15)$, $(25,-5)$, $(10,15)$). $\mathrm{MSE}_{\text{within}}$ per $K=1,\dots,8$: $264{,}5;\ 136{,}7;\ 67{,}8;\ 36{,}0;\ 4{,}28;\ 3{,}95;\ 3{,}66;\ 3{,}45$. Dopo $K=5$ il calo è minimo: gomito a $K=5$.

Gap statistic

Il gomito è un criterio visivo e soggettivo. La gap statistic confronta il clustering ottenuto con quello che ci si aspetterebbe da dati senza struttura: «il clustering osservato è migliore di quello su dati casuali?»

1. Per ogni $K$ (per esempio da 1 a 10) si esegue il k-means sui dati reali e si calcola la dispersione $W_K=\sum_{k=1}^K\sum_{x\in C_k}\lVert x-\mu_k\rVert^2$.
2. Si generano $n_{\text{ref}}$ dataset di riferimento con lo stesso numero di campioni e dimensioni e gli stessi limiti dei dati originali, ma con valori estratti da una distribuzione uniforme (nessuna struttura); servono più dataset per non dipendere da una singola realizzazione.
3. Si esegue il k-means su ciascuno e si calcola $W_K^{\text{ref}}$.
4. $\displaystyle\mathrm{Gap}(K)=\mathbb E\big[\log W_K^{\text{ref}}\big]-\log W_K$: più alto è il gap, migliore è il clustering rispetto al caso. Si usa il logaritmo perché interessa il rapporto tra le due dispersioni.
5. Si sceglie il più piccolo $K$ tale che $\mathrm{Gap}(K)\ge\mathrm{Gap}(K+1)-s_{K+1}$, dove $s_{K+1}$ è la deviazione standard dei $\log W^{\text{ref}}$ (tiene conto dell'incertezza dei dati casuali).

Esempio (stessi cinque gruppi). Gap per $K=1,\dots,6$: $0{,}04;\ 0{,}16;\ 0{,}41;\ 0{,}64;\ 2{,}60;\ 2{,}47$ (laboratorio). Per $K=4$: $0{,}64\ge2{,}60-s$? No. Per $K=5$: $2{,}60\ge2{,}47-s_6$ sì: si sceglie $K=5$. (Una replica indipendente con 10 dataset di riferimento ha dato $0{,}02;\ 0{,}13;\ 0{,}41;\ 0{,}62;\ 2{,}56;\ 2{,}45$.)

Clustering gerarchico

Definizione (clustering gerarchico). Costruisce un albero di cluster, il dendrogramma, in due modi:

• agglomerativo (bottom-up): ogni punto parte come cluster a sé e si fondono cluster passo dopo passo;
• divisivo (top-down): si parte da un unico cluster che contiene tutto e lo si divide passo dopo passo.

Un dendrogramma è un diagramma ad albero che mostra il processo di fusione (l'altezza di ciascuna fusione è la distanza tra i cluster fusi); lo si «taglia» a un certo livello per decidere quanti cluster ottenere.

Agglomerativo, algoritmo. (1) Ogni punto è un cluster. (2) Si calcolano le distanze tra tutti i cluster. (3) Si fondono i due più vicini. (4) Si aggiornano le distanze. (5) Si ripete finché resta un solo cluster (o $m$ cluster). Con $n$ campioni e $m$ cluster desiderati servono $n-m$ iterazioni, perché a ogni fusione il numero di cluster diminuisce di uno.

Distanza tra cluster (linkage). Per due cluster $C_i,C_j$:

• single linkage (distanza minima): $d(C_i,C_j)=\min_{a\in C_i,\,b\in C_j}d(a,b)$;
• complete linkage (distanza massima): $\max_{a,b}d(a,b)$;
• average linkage (distanza media): $\frac1{|C_i||C_j|}\sum_{a\in C_i}\sum_{b\in C_j}d(a,b)$;
• Ward: fonde la coppia che aumenta meno la varianza totale entro i cluster.

Esempio svolto (average e single linkage). Sei punti $P_0=(1,2)$, $P_1=(5,8)$, $P_2=(1,1)$, $P_3=(2,3)$, $P_4=(1{,}5;1{,}8)$, $P_5=(5;6{,}5)$ e distanza euclidea. Ordine delle fusioni con average linkage (distanza di fusione tra parentesi): $\{P_0,P_4\}$ $(0{,}539)$; poi $\{P_0,P_4\}$ con $P_2$ $(0{,}972)$; $\{P_1,P_5\}$ $(1{,}50)$; $\{P_0,P_4,P_2\}$ con $P_3$ $(1{,}65)$; infine i due gruppi $\{P_0,P_2,P_3,P_4\}$ e $\{P_1,P_5\}$ $(6{,}44)$. Con single linkage le distanze sono $0{,}539;\ 0{,}943;\ 1{,}30;\ 1{,}50;\ 4{,}61$: stessi gruppi, altezze diverse. Il salto grande nell'ultima fusione indica che due cluster sono la scelta naturale.

Il laboratorio memorizza il processo con una matrice di assegnazione: una riga per iterazione e una colonna per campione; l'elemento $(i,j)$ è l'etichetta del cluster del campione $j$ all'iterazione $i$ (la riga $i$ contiene $n-i$ cluster); la riga iniziale è $[0,1,\dots,n-1]$. I centroidi si ottengono come baricentri dei cluster.

Forme dei gruppi. Con linkage single le fusioni seguono catene di punti vicini e si riescono a separare forme allungate (due «lune», due anelli concentrici); con average o complete (e come il k-means) i cluster tendono a essere compatti e le forme non convesse non sono separate bene. Lo svantaggio del single è l'effetto «catena» (un ponte di punti unisce due gruppi).

Divisivo. In ogni passo si sceglie il cluster da dividere (per esempio quello con la maggiore distanza media tra le coppie di punti, cioè la dispersione più alta) e lo si divide in due con un criterio (per esempio con il k-means con $K=2$). Per ottenere $m$ cluster servono $m-1$ divisioni (Esercizio - Clustering gerarchico agglomerativo e divisivo).

Costo. La versione ingenua confronta tutte le coppie di cluster a ogni passo (almeno quadratico per ogni passo, circa cubico in totale): è molto più costosa del k-means per $n$ grande.

Codice

import numpy as np
def kmeans(X, k, max_iters=100, tol=1e-6, rng=np.random):
    centroids = X[rng.choice(len(X), k, replace=False)]         # K punti a caso
    for _ in range(max_iters):
        d = np.linalg.norm(X[:, None] - centroids[None], axis=2)   # (n, k) distanze
        labels = d.argmin(axis=1)                                  # assegnazione
        new = np.array([X[labels == j].mean(axis=0) for j in range(k)])   # medie
        if np.all(np.abs(new - centroids) < tol): break
        centroids = new
    return labels, centroids

def within_dispersion(X, labels, centroids):                     # W_K per gomito e gap
    return sum(((X[labels == j] - centroids[j]) ** 2).sum() for j in range(len(centroids)))

Attenzione: un cluster che resta senza punti dà una media non definita (nan); in tal caso si reinizializza quel centroide. Con scikit-learn: KMeans(n_clusters=5, n_init=10).fit(X) e AgglomerativeClustering(n_clusters=3, linkage="average").

Errori tipici

• Scegliere $K$ minimizzando l'errore entro i cluster: cala sempre, e per $K=N$ vale zero.
• Lanciare il k-means una sola volta senza guardare l'inizializzazione.
• Usare il k-means senza standardizzare le feature, o su gruppi non sferici.
• Scrivere nella gap statistic $\mathrm{Gap}$ con il segno invertito: si sottrae $\log W_K$ dei dati reali dalla media dei $\log W^{\text{ref}}$.
• Contare male le iterazioni agglomerative: servono $n-m$ fusioni per $m$ cluster.
• Scordare che il dendrogramma si taglia a un'altezza: il numero di cluster lo decide chi usa l'algoritmo.