39 - Riduzione di dimensionalità: PCA da zero e t-SNE

Una lampada proietta sul muro l'ombra di un oggetto tridimensionale. Se la lampada è puntata dal lato sbagliato, l'ombra di una bottiglia è un cerchio e non capisci cosa sia; se la punti dal lato giusto, l'ombra conserva quasi tutta la forma. La riduzione di dimensionalità cerca il lato giusto da cui «proiettare» dati con molte feature su poche dimensioni, perdendo il meno possibile. Serve a comprimere, a velocizzare altri modelli, a togliere ridondanza e, con 2 o 3 dimensioni, a guardare i dati.

Due strumenti, due filosofie. La PCA è lineare: ruota gli assi per mettere la maggior varianza sui primi. La t-SNE è non lineare e punta a un solo obiettivo, la visualizzazione: conservare i vicinati. Le implementiamo entrambe a mano e vediamo anche gli errori di interpretazione che inducono.

Gli ingredienti, uno per uno

• Varianza: quanto i valori di una feature si disperdono attorno alla loro media. Una feature con varianza zero è una costante e non porta informazione.
• Covarianza: quanto due feature variano insieme. Positiva se quando una cresce cresce anche l'altra, negativa se cresce l'una e cala l'altra, zero se sono scorrelate. La matrice di covarianza C raccoglie tutte le covarianze a coppie; sulla diagonale ci sono le varianze.
• Centrare i dati: sottrarre a ogni colonna la sua media, così la nuvola di punti è centrata nell'origine.
• Autovettore e autovalore: una matrice trasforma i vettori ruotandoli e stirandoli. Un autovettore è una direzione speciale che la matrice non ruota, ma solo allunga o accorcia; l'autovalore è il fattore di allungamento. Per la matrice di covarianza, gli autovettori sono gli assi naturali della nuvola di dati e gli autovalori dicono quanto la nuvola è «lunga» lungo ciascuno.

Le quattro voci si concatenano in un solo procedimento: centri i dati, calcoli la matrice di covarianza, ne cerchi gli autovettori, e gli autovalori ti dicono quanta varianza porta ciascun asse. Se capisci questa sequenza, la PCA è solo la sua applicazione; il resto della lezione la verifica con i numeri.

Il disegno: una nuvola allungata in diagonale. La prima componente principale (PC1, rossa) è la direzione lungo cui la nuvola è più estesa; la seconda (PC2, verde) è perpendicolare. Ridurre a una dimensione significa proiettare ogni punto su PC1: si perde solo la dispersione lungo PC2, che è piccola.

PCA: trovare la direzione di massima varianza

Dati X con n righe e d colonne, già centrati (sottratta la media di ogni colonna). Cerchiamo una direzione, cioè un vettore unitario w (‖w‖ = 1), lungo la quale i dati proiettati abbiano varianza massima. La proiezione di un punto è x·w e la varianza delle proiezioni è

Var(Xw) = (1/(n−1)) · (Xw)ᵀ(Xw) = wᵀ C w          con C = XᵀX / (n−1)   (matrice di covarianza)

Il problema è: massimizzare wᵀCw con il vincolo wᵀw = 1. Con un moltiplicatore di Lagrange λ:

L(w, λ) = wᵀCw − λ (wᵀw − 1)
∂L/∂w = 2Cw − 2λw = 0     →     C w = λ w

Le direzioni stazionarie sono quindi gli autovettori di C. E la varianza lungo ciascuna è wᵀCw = wᵀ(λw) = λ: l'autovalore è la varianza. La prima componente principale è l'autovettore con l'autovalore maggiore; la seconda quello successivo (ortogonale alla prima, perché C è simmetrica), e così via. La quota di varianza spiegata dalla componente i è λ_i / Σλ.

Per ridurre a k dimensioni: Z = X W_k, dove W_k ha come colonne i primi k autovettori. Per ricostruire: X̂ = Z W_kᵀ (+ la media). Si può dimostrare che questa scelta minimizza anche l'errore di ricostruzione quadratico: i due punti di vista (massima varianza, minimo errore) coincidono. Qui lo verifichiamo numericamente invece di dimostrarlo.

Implementazione a mano sulle cifre scritte a mano

Il dataset Digits contiene 1797 immagini 8×8 di cifre, cioè 64 feature per immagine.

import numpy as np
from sklearn.datasets import load_digits, load_wine
from sklearn.decomposition import PCA

X, y = load_digits(return_X_y=True)
print("forma:", X.shape)

mu = X.mean(axis=0)
Xc = X - mu                                         # 1) centrare
C = Xc.T @ Xc / (len(X) - 1)                        # 2) matrice di covarianza (64 x 64)
valori, vettori = np.linalg.eigh(C)                 # 3) autovalori/autovettori (C simmetrica)
ordine = np.argsort(valori)[::-1]                   #    dal più grande
valori, vettori = valori[ordine], vettori[:, ordine]
quota = valori / valori.sum()

sk = PCA().fit(X)
print("quote di varianza uguali a sklearn:", np.allclose(quota, sk.explained_variance_ratio_))
print("differenza massima tra componenti (in valore assoluto):", np.abs(np.abs(vettori[:, :5].T) - np.abs(sk.components_[:5])).max())
print("prime 5 quote di varianza:", np.round(quota[:5], 3))

forma: (1797, 64)
quote di varianza uguali a sklearn: True
differenza massima tra componenti (in valore assoluto): 3.469446951953614e-15
prime 5 quote di varianza: [0.149 0.136 0.118 0.084 0.058]

Nel listato si centrano i dati (Xc), si costruisce la matrice di covarianza C 64×64 con la formula XcᵀXc/(n−1), e si risolve il problema agli autovalori con np.linalg.eigh, la variante per matrici simmetriche (più veloce e stabile di eig, e garantisce autovalori reali e autovettori ortogonali). eigh restituisce gli autovalori in ordine crescente, perciò si riordina argsort(...)[::-1]. La quota di varianza della componente i è valori[i] / valori.sum(). Le tre stampe confrontano con sklearn.decomposition.PCA: le quote coincidono e le componenti differiscono di ordine 10⁻¹⁵ in valore assoluto.

Il confronto usa il valore assoluto perché un autovettore è definito a meno del segno: w e −w sono la stessa direzione, e librerie diverse possono restituire l'uno o l'altro. Non è un errore, ma può sorprenderti se confronti i numeri a occhio.

Quante componenti tenere

cum = np.cumsum(quota)
for soglia in (0.80, 0.90, 0.95, 0.99):
    print("componenti per spiegare il %2d%% della varianza: %d su 64" % (soglia * 100, np.searchsorted(cum, soglia) + 1))

for k in (2, 5, 10, 20, 40):
    W = vettori[:, :k]
    errore = ((Xc - Xc @ W @ W.T) ** 2).sum(axis=1).mean()
    teorico = (valori[k:].sum()) * (len(X) - 1) / len(X)
    print("k=%2d  errore di ricostruzione medio %8.2f   varianza scartata %8.2f" % (k, errore, teorico))

componenti per spiegare il 80% della varianza: 13 su 64
componenti per spiegare il 90% della varianza: 21 su 64
componenti per spiegare il 95% della varianza: 29 su 64
componenti per spiegare il 99% della varianza: 41 su 64
k= 2  errore di ricostruzione medio   858.94   varianza scartata   858.94
k= 5  errore di ricostruzione medio   546.72   varianza scartata   546.72
k=10  errore di ricostruzione medio   314.51   varianza scartata   314.51
k=20  errore di ricostruzione medio   126.99   varianza scartata   126.99
k=40  errore di ricostruzione medio    14.17   varianza scartata    14.17

Con np.searchsorted(cum, soglia) si trova la prima posizione in cui la varianza cumulata supera la soglia (da cui +1 per contare le componenti). Poi per ogni k si costruisce W con i primi k autovettori (matrice 64×k), si proietta (Xc @ W: n×k), si ricostruisce (@ W.T: di nuovo n×64) e si misura l'errore quadratico medio per campione; lo si confronta con la somma degli autovalori scartati, moltiplicata per (n−1)/n per allineare la normalizzazione.

L'errore di ricostruzione coincide con la varianza delle componenti scartate: è la prova numerica che massima varianza e minimo errore sono la stessa cosa. Con 29 componenti su 64 si conserva il 95% della varianza: le cifre hanno molta ridondanza tra pixel adiacenti.

S = [ 5  2 ]
    [ 2  2 ]

PCA come passo di una pipeline

Ridurre le dimensioni non peggiora sempre un classificatore, e può accelerarlo. Misuriamolo con validazione incrociata, con PCA dentro la pipeline (così viene rifatta su ogni fold, senza leakage):

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

for k in (2, 5, 10, 20, 40, 64):
    pipe = make_pipeline(StandardScaler(), PCA(n_components=k), LogisticRegression(max_iter=2000))
    print("componenti %2d  accuratezza CV %.3f" % (k, cross_val_score(pipe, X, y, cv=5).mean()))

componenti  2  accuratezza CV 0.535
componenti  5  accuratezza CV 0.771
componenti 10  accuratezza CV 0.840
componenti 20  accuratezza CV 0.899
componenti 40  accuratezza CV 0.914
componenti 64  accuratezza CV 0.920

La Pipeline concatena standardizzazione, PCA con n_components=k e regressione logistica (un classificatore lineare, descritto in regressione logistica); cross_val_score con cv=5 la addestra e valuta cinque volte su suddivisioni diverse, e .mean() riporta la media. Ogni volta lo scaler e la PCA vedono solo il training di quel fold.

Con 2 componenti l'accuratezza è 0.535, con 10 sale a 0.840, con 20 a 0.899 e con 40 a 0.914, contro 0.920 con tutte le 64: le ultime 24 componenti aggiungono mezzo punto, il resto era in gran parte ridondanza. Il numero di componenti è un iperparametro, da scegliere in validazione.

La trappola della scala

La PCA cerca la varianza, e la varianza dipende dall'unità di misura. Una colonna in migliaia di unità ha una varianza gigantesca e si prende la prima componente, qualunque cosa dica. Sul dataset Wine:

Xw, yw = load_wine(return_X_y=True)
print("senza scaling, quota delle prime 3 componenti:", np.round(PCA().fit(Xw).explained_variance_ratio_[:3], 4))
print("con standardizzazione:                        ", np.round(PCA().fit(StandardScaler().fit_transform(Xw)).explained_variance_ratio_[:3], 4))

senza scaling, quota delle prime 3 componenti: [9.981e-01 1.700e-03 1.000e-04]
con standardizzazione:                         [0.362  0.1921 0.1112]

Il listato confronta le quote di varianza spiegata delle prime tre componenti con e senza standardizzazione. La differenza è il segno che, senza scaling, la PCA sta ordinando le feature per unità di misura, non per contenuto informativo.

Senza scaling la prima componente «spiega» il 99,8% della varianza, ma è solo la colonna proline (valori nell'ordine delle migliaia) travestita da componente. Dopo la standardizzazione la varianza si distribuisce tra le componenti (36%, 19%, 11%). La regola: standardizza quando le feature hanno unità diverse; non farlo solo quando le unità sono omogenee e la scala ha significato (per esempio i pixel di un'immagine).

t-SNE: conservare i vicinati per visualizzare

La PCA è lineare: se i dati stanno su una superficie curva, la proiezione lineare la schiaccia. La t-SNE (t-distributed Stochastic Neighbor Embedding) costruisce una mappa 2D in cui i punti vicini nello spazio originale restano vicini. L'idea ha tre passi.

1. Vicinati nello spazio originale. Per ogni punto i, la probabilità che j sia suo vicino è una gaussiana centrata su i:

p(j|i) = exp(−‖x_i − x_j‖² / 2σ_i²) / Σ_{k≠i} exp(−‖x_i − x_k‖² / 2σ_i²)
p_ij = ( p(j|i) + p(i|j) ) / 2n

La larghezza σ_i si sceglie per ogni punto, in modo che il punto abbia un numero «efficace» di vicini pari alla perplexity: Perp = 2^H(P_i), con H l'entropia di Shannon della distribuzione p(·|i). Nelle zone dense σ_i è piccola, in quelle rade è grande. Si trova con una ricerca binaria.

2. Vicinati nella mappa. Per i punti y_i della mappa si usa una distribuzione t di Student con un grado di libertà, a code pesanti:

q_ij = (1 + ‖y_i − y_j‖²)⁻¹ / Σ_{k≠l} (1 + ‖y_k − y_l‖²)⁻¹

Le code pesanti risolvono il «problema dell'affollamento»: in 2D non c'è spazio per mettere vicini tutti i vicini moderati di un punto; con la t, punti moderatamente lontani possono stare molto lontani nella mappa senza pagare molto.

3. Ottimizzazione. Si sposta y per minimizzare la divergenza di Kullback-Leibler tra P e Q, KL(P‖Q) = Σ p_ij log(p_ij/q_ij), con la discesa del gradiente. Il gradiente ha una forma compatta:

∂KL/∂y_i = 4 Σ_j (p_ij − q_ij) (y_i − y_j) (1 + ‖y_i − y_j‖²)⁻¹

Ogni coppia agisce come una molla: attrae i punti se p_ij > q_ij (sono vicini nei dati ma lontani nella mappa), li respinge nel caso opposto.

Implementazione a mano (500 cifre)

Il listato ha due funzioni: affinita costruisce la matrice P con la ricerca binaria di σ_i per ogni punto, e tsne ottimizza le posizioni 2D con il gradiente scritto sopra. Alla fine confronta la nostra mappa con quella di scikit-learn e con la semplice PCA a due componenti, usando la trustworthiness (quanto i vicini nella mappa sono vicini anche nei dati originali).

from sklearn.manifold import TSNE, trustworthiness

def affinita(D, perplexity):
    n = len(D); P = np.zeros((n, n)); target = np.log(perplexity)
    for i in range(n):
        d = np.delete(D[i], i); d = d - d.min()                 # traslare non cambia la distribuzione
        lo, hi, beta = 0.0, np.inf, 1.0 / d.mean()              # beta = 1 / (2 sigma^2)
        for _ in range(60):                                     # ricerca binaria su beta
            pr = np.exp(-beta * d); s = pr.sum()
            H = np.log(s) + beta * (d * pr).sum() / s           # entropia (in nat)
            if H > target: lo, beta = beta, (beta * 2 if np.isinf(hi) else (beta + hi) / 2)
            else:          hi, beta = beta, (beta + lo) / 2
        P[i, np.arange(n) != i] = np.exp(-beta * d) / np.exp(-beta * d).sum()
    return np.maximum((P + P.T) / (2 * n), 1e-12)

def tsne(X, perplexity=30, n_iter=500, lr=200.0, seed=0):
    n = len(X)
    D = ((X[:, None, :] - X[None, :, :]) ** 2).sum(axis=2)
    P = affinita(D, perplexity)
    rng = np.random.default_rng(seed)
    Y = rng.normal(0, 1e-4, (n, 2)); V = np.zeros_like(Y)
    for it in range(n_iter):
        esagerazione = 12.0 if it < 100 else 1.0                # i primi passi: P ingrandita per formare i cluster
        num = 1.0 / (1.0 + ((Y[:, None, :] - Y[None, :, :]) ** 2).sum(axis=2))
        np.fill_diagonal(num, 0.0)
        Q = np.maximum(num / num.sum(), 1e-12)
        PQ = (esagerazione * P - Q) * num
        grad = 4 * (PQ.sum(axis=1)[:, None] * Y - PQ @ Y)
        V = (0.5 if it < 250 else 0.8) * V - lr * grad          # discesa con momento
        Y = Y + V; Y = Y - Y.mean(axis=0)
    return Y, (P * np.log(P / Q)).sum()

idx = np.random.default_rng(0).choice(len(X), 500, replace=False)
Xs, ys = X[idx], y[idx]

Y_mia, kl = tsne(Xs)
Y_sk = TSNE(perplexity=30, init="pca", random_state=0).fit(Xs)
Y_pca = PCA(n_components=2).fit_transform(Xs)

print("KL finale a mano: %.3f    sklearn: %.3f" % (kl, Y_sk.kl_divergence_))
for nome, Y in (("PCA 2D", Y_pca), ("t-SNE a mano", Y_mia), ("t-SNE sklearn", Y_sk.embedding_)):
    print("%-14s trustworthiness (10 vicini): %.3f" % (nome, trustworthiness(Xs, Y, n_neighbors=10)))

KL finale a mano: 0.550    sklearn: 0.519
PCA 2D         trustworthiness (10 vicini): 0.830
t-SNE a mano   trustworthiness (10 vicini): 0.981
t-SNE sklearn  trustworthiness (10 vicini): 0.983

Nel listato affinita costruisce la matrice P: per ogni punto i esegue la ricerca binaria su beta (= 1/2σ²) finché l'entropia della distribuzione dei vicini eguaglia il logaritmo della perplexity; poi simmetrizza e normalizza. tsne inizializza i punti 2D con rumore minuscolo, e per 500 iterazioni calcola Q con la distribuzione t, il gradiente nella forma compatta, e aggiorna le posizioni con la discesa del gradiente con momento (il momento accumula una frazione del passo precedente per attraversare le zone piatte).

Nelle prime 100 iterazioni P è moltiplicata per 12 (early exaggeration): ingrandire le attrazioni aiuta i cluster a formarsi prima di affinarli. Il confronto finale usa 500 cifre estratte a caso.

La trustworthiness misura quanto i vicini nella mappa lo sono anche nello spazio originale (1 = perfetto). La proiezione lineare della PCA a 2 dimensioni ottiene 0.830, la t-SNE a mano 0.981 e quella di scikit-learn 0.983: conservare i vicinati è esattamente lo scopo della t-SNE. La nostra implementazione arriva vicina a quella di scikit-learn (KL finale 0.550 contro 0.519), che inizializza i punti con la PCA (init="pca"), un'opzione che la guida di scikit-learn raccomanda per preservare meglio la struttura globale. Per i dettagli di questo esperimento: 500 cifre invece di 1797, per tenere i tempi del calcolo a mano accettabili.

Come NON leggere una t-SNE

La mappa è bella e ingannevole. Gli errori di lettura più comuni:

• Le distanze tra cluster non significano niente: la t-SNE conserva i vicinati locali, non le distanze globali. Due cluster lontani nella mappa non sono necessariamente «più diversi» di due vicini.
• Le dimensioni dei cluster non significano niente: la normalizzazione per densità (σ_i) espande i cluster densi e comprime quelli radi.
• La forma dipende dagli iperparametri: la stessa t-SNE con perplexity diversa disegna mappe di aspetto diverso (cluster che si spezzano o si fondono; guardale tracciando i grafici con le etichette a colori).

I tre errori hanno la stessa origine: la mappa t-SNE conserva soltanto i vicinati locali, e tutto il resto (distanze, dimensioni, forme) è un effetto collaterale dell'ottimizzazione e degli iperparametri. Una figura t-SNE si legge per chi sta con chi, mai per quanto è lontano. Dell'effetto degli iperparametri misuriamo ciò che si può misurare in una tabella:

for perp in (2, 5, 30, 100):
    t = TSNE(perplexity=perp, init="pca", random_state=0).fit(Xs)
    print("perplexity %3d  KL %.3f   trustworthiness %.3f" % (perp, t.kl_divergence_, trustworthiness(Xs, t.embedding_, n_neighbors=10)))

perplexity   2  KL 0.583   trustworthiness 0.970
perplexity   5  KL 0.597   trustworthiness 0.980
perplexity  30  KL 0.519   trustworthiness 0.983
perplexity 100  KL 0.405   trustworthiness 0.984

Sul piano dei vicinati locali la trustworthiness resta alta per tutte le perplexity (da 0.970 a 0.984): la qualità «locale» è simile, ma la disposizione globale delle isole può cambiare, e questo la tabella non lo mostra. Attenzione al valore di KL: quando cambia la perplexity cambia la distribuzione P stessa, quindi le divergenze non sono confrontabili tra una perplexity e l'altra. Una KL più bassa con perplexity 100 non significa che la mappa sia «migliore». Un'altra regola: la perplexity deve essere minore del numero di punti, e valori tipici stanno tra 5 e 50. La guida di scikit-learn definisce la perplexity come 2 elevato all'entropia di Shannon della distribuzione condizionale, cioè il numero effettivo di vicini considerati.

Il listato ripete la t-SNE di scikit-learn con quattro perplexity (sugli stessi 500 campioni) e stampa KL e trustworthiness: è la tabella appena commentata. Restano tre limiti di uso, non di lettura:

• Nessuna funzione di trasformazione: la t-SNE ottimizza direttamente le posizioni dei punti dati; non c'è un modo di proiettare un punto nuovo. La PCA ha una matrice W e proietta qualunque punto.
• Non usarla come passo prima di un clustering o un classificatore: è pensata per guardare, non come estrazione di feature.
• Varia con il seme: ripeti con più semi prima di fidarti di una struttura.

Questi tre limiti, uniti ai precedenti, spiegano perché t-SNE è uno strumento di esplorazione e non un passo di pipeline: la guida di scikit-learn lo dice in termini diversi, ricordando che il suo scopo principale è la visualizzazione di dati ad alta dimensione, che l'algoritmo è stocastico e che la struttura globale non è esplicitamente preservata.

PCA contro t-SNE: quale scegliere

• PCA: veloce, deterministica (a meno del segno), lineare, invertibile (puoi ricostruire), proietta nuovi dati, le componenti sono interpretabili tramite i pesi (loadings). Non cattura strutture curve.
• t-SNE: cattura strutture non lineari, ottima per esplorare; lenta (O(n²) nella versione esatta), stocastica, nessuna proiezione di punti nuovi, distanze globali non interpretabili.
• Prassi comune: PCA a 30-50 componenti per togliere il rumore e velocizzare, poi t-SNE su quelle per visualizzare.

La scelta si riduce a una domanda: ti serve una trasformazione che applichi anche a dati nuovi e che puoi invertire, oppure una mappa per guardare? Nel primo caso PCA, nel secondo t-SNE. Se vuoi entrambe le cose, la PCA riduce prima il rumore e la t-SNE disegna dopo, ma la mappa t-SNE non entra in produzione.

Trappole reali

• PCA senza scaling con unità diverse (vedi sopra).
• PCA fit su tutto il dataset prima della divisione train/test: leakage. Dentro la pipeline.
• Scambiare «varianza» per «informazione utile alla classificazione»: la direzione di massima varianza può essere irrilevante per la classe che vuoi prevedere. La PCA non guarda le etichette.
• Interpretare gli assi della t-SNE: gli assi non hanno significato.
• Segno delle componenti: due esecuzioni possono dare componenti con segno opposto. Se salvi un modello, salva anche la matrice W.

Le cinque trappole sono di due tipi. Quelle sulla PCA (scala, leakage, varianza contro utilità, segno) falsano il risultato senza segnalarlo; quella sugli assi della t-SNE è un errore di lettura. Per tutte vale la regola della lezione: misura prima di interpretare, e confronta la pipeline con quella senza riduzione.

Tre modi di calcolare la PCA: autovalori, SVD e scikit-learn

Finora la PCA è stata presentata come diagonalizzazione della matrice di covarianza. In pratica non si calcola quasi mai la covarianza esplicitamente: si usa la decomposizione ai valori singolari (SVD) della matrice dei dati centrati X_c = U Σ Vᵀ. Il legame è diretto: la covarianza è X_cᵀ X_c / (n − 1) = V Σ² Vᵀ / (n − 1), quindi le colonne di V (i «vettori singolari destri») sono gli autovettori della covarianza, e gli autovalori sono σ_i² / (n − 1). L'SVD è preferita perché è numericamente più stabile (non eleva al quadrato i numeri di condizionamento) e non richiede di costruire la matrice d × d, che con migliaia di feature sarebbe enorme. Verifichiamo che i tre percorsi diano la stessa cosa sulle cifre.

import warnings; warnings.filterwarnings("ignore")
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression

X, y = load_digits(return_X_y=True)
Xc = X - X.mean(axis=0)                                                  # dati centrati

cov = Xc.T @ Xc / (len(X) - 1)
eig_values, eig_vectors = np.linalg.eigh(cov)                            # percorso 1: autovalori della covarianza
order = np.argsort(eig_values)[::-1]
eig_values, eig_vectors = eig_values[order], eig_vectors[:, order]

U, S, Vt = np.linalg.svd(Xc, full_matrices=False)                        # percorso 2: SVD dei dati
svd_values = S ** 2 / (len(X) - 1)

pca = PCA().fit(X)                                                       # percorso 3: scikit-learn

print("differenza massima tra le prime 20 varianze, autovalori vs SVD: %.1e" % np.abs(eig_values[:20] - svd_values[:20]).max())
print("differenza massima, scikit-learn vs SVD: %.1e" % np.abs(pca.explained_variance_[:20] - svd_values[:20]).max())
sign = np.sign((eig_vectors[:, :10] * Vt[:10].T).sum(axis=0))             # gli autovettori sono definiti a meno del segno
print("differenza massima tra i primi 10 vettori (dopo aver allineato il segno): %.1e" % np.abs(eig_vectors[:, :10] * sign - Vt[:10].T).max())
print("rango della matrice centrata:", np.linalg.matrix_rank(Xc), "  pixel a varianza zero:", int((X.std(axis=0) == 0).sum()))

differenza massima tra le prime 20 varianze, autovalori vs SVD: 4.5e-13
differenza massima, scikit-learn vs SVD: 2.6e-13
differenza massima tra i primi 10 vettori (dopo aver allineato il segno): 1.2e-14
rango della matrice centrata: 61   pixel a varianza zero: 3

Il listato calcola la stessa decomposizione in tre modi: gli autovalori della covarianza esplicita (eigh, che sfrutta la simmetria della matrice), l'SVD dei dati centrati, e la classe PCA di scikit-learn (che internamente usa l'SVD). Le prime 20 varianze coincidono fino a 4,5·10⁻¹³ tra autovalori e SVD e fino a 2,6·10⁻¹³ tra scikit-learn e SVD: differenze di arrotondamento in virgola mobile, non di sostanza. Gli autovettori coincidono fino a 1,2·10⁻¹⁴ dopo aver allineato il segno: ogni autovettore v è valido quanto −v, e librerie diverse possono restituirlo con segni opposti, quindi un test che confronta le componenti senza allineare il segno è un test fragile.

L'ultima riga è un promemoria utile: la matrice centrata ha rango 61 e non 64, perché 3 pixel dei bordi sono sempre zero in tutto il dataset (varianza nulla) e le ultime tre componenti non portano informazione. Se non conoscessi questo fatto, la sbiancatura qui sotto ti sorprenderebbe.

Eliminare il rumore con la PCA e la trappola della sbiancatura

La PCA ha un secondo uso, oltre alla riduzione di dimensione per i modelli: la ricostruzione. Se proietti un'immagine sulle prime k componenti e poi torni indietro con inverse_transform, ottieni la migliore approssimazione di rango k (nel senso dei minimi quadrati). Il rumore, che è sparso in modo uniforme in tutte le direzioni, viene in gran parte scartato insieme alle componenti a bassa varianza, mentre il segnale, concentrato nelle prime componenti, resta. Misuriamolo: aggiungiamo rumore gaussiano di deviazione 4 (su pixel che vanno da 0 a 16) alle immagini di test e ricostruiamo con k componenti calcolate sul solo training set.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)
rng3 = np.random.default_rng(0)
noisy = X_test + rng3.normal(0, 4, X_test.shape)
rmse = lambda a, b: np.sqrt(((a - b) ** 2).mean())

print("errore (RMSE) rispetto alle immagini pulite; il rumore da solo vale %.2f" % (rmse(noisy, X_test)))
print("k    RMSE dopo la ricostruzione   RMSE se l'immagine fosse già pulita")
for k in (2, 5, 10, 20, 40, 64):
    pc = PCA(k).fit(X_train)
    denoised = pc.inverse_transform(pc.transform(noisy))
    clean_rec = pc.inverse_transform(pc.transform(X_test))
    print("%-4d  %-27.3f  %.3f" % (k, rmse(denoised, X_test), rmse(clean_rec, X_test)))

errore (RMSE) rispetto alle immagini pulite; il rumore da solo vale 4.00
k    RMSE dopo la ricostruzione   RMSE se l'immagine fosse già pulita
2     3.743                        3.676
5     3.160                        2.950
10    2.750                        2.247
20    2.668                        1.448
40    3.196                        0.485
64    3.996                        0.000

Il rumore da solo vale un errore quadratico medio di 4,00 (è la deviazione che abbiamo scelto). Poi per ogni k adatta la PCA sul training set, proietta le immagini rumorose, le ricostruisce e misura l'errore rispetto alle immagini pulite. L'errore scende da 3,743 con k = 2 a un minimo di 2,668 con k = 20, poi risale: 3,196 con k = 40 e 3,996 con k = 64, dove la ricostruzione è l'identità e il rumore torna intero.

Questo è il compromesso bias-varianza in veste di compressione: con k piccolo butti anche il segnale (la colonna di destra, l'errore su immagini già pulite, parte da 3,676 e scende fino a zero con 64 componenti), con k grande tieni anche il rumore. L'ottimo, k = 20, abbatte l'errore da 4,00 a 2,67 (circa un terzo in meno) usando solo 20 numeri per immagine. Nota che il k ottimo per l'eliminazione del rumore non è quello che «conserva il 95% della varianza»: dipende dalla quantità di rumore, ed è un altro parametro da scegliere con dati di validazione.

Sbiancare (whitening): utile e pericoloso

La sbiancatura divide ogni componente per la sua deviazione standard, così che le componenti abbiano tutte varianza 1 e siano incorrelate: la covarianza dei dati trasformati diventa la matrice identità. Alcuni modelli (quelli basati sulla distanza, o con regolarizzazione isotropa) ne traggono vantaggio perché nessuna direzione domina. Ma c'è un rovescio: dividere per una deviazione piccola amplifica proprio le direzioni a minore varianza, che sono quelle dominate dal rumore, e se una deviazione è praticamente zero la divisione esplode.

pw = PCA(20, whiten=True).fit(X_train)
cov_w = np.cov(pw.transform(X_train).T)
print("covarianza dopo la sbiancatura: diagonale tra %.3f e %.3f, massimo elemento fuori diagonale %.1e" % (cov_w.diagonal().min(), cov_w.diagonal().max(), np.abs(cov_w - np.diag(cov_w.diagonal())).max()))

print("\naccuratezza sul test, con e senza sbiancatura:")
print("k    KNN senza   KNN con   regr. logistica senza   regr. logistica con")
for k in (10, 20, 40, 64):
    row = []
    for model in (KNeighborsClassifier(5), LogisticRegression(max_iter=10000)):
        for whiten in (False, True):
            row.append(make_pipeline(PCA(k, whiten=whiten), model).fit(X_train, y_train).score(X_test, y_test))
    print("%-4d  %-10.3f  %-8.3f  %-22.3f  %.3f" % (k, row[0], row[1], row[2], row[3]))

covarianza dopo la sbiancatura: diagonale tra 1.000 e 1.000, massimo elemento fuori diagonale 1.6e-15

accuratezza sul test, con e senza sbiancatura:
k    KNN senza   KNN con   regr. logistica senza   regr. logistica con
10    0.969       0.956     0.941                   0.939
20    0.974       0.972     0.948                   0.948
40    0.980       0.957     0.961                   0.965
64    0.980       0.237     0.967                   0.961

La prima riga conferma la definizione: dopo la sbiancatura di 20 componenti la covarianza dei dati trasformati ha diagonale esattamente 1,000 e il più grande elemento fuori diagonale vale 1,6·10⁻¹⁵ (zero a meno dell'arrotondamento). Poi confronta quattro configurazioni con 10, 20, 40 e 64 componenti. Con poche componenti la sbiancatura non cambia quasi nulla: con 20 componenti il KNN passa da 0,974 a 0,972 e la regressione logistica resta a 0,948. Con 40 il KNN peggiora (0,980 → 0,957) mentre la regressione logistica migliora di poco (0,961 → 0,965): effetti piccoli e di segno diverso, che su un test di 540 immagini non sono distinguibili dal caso.

Con 64 componenti accade il disastro: il KNN con sbiancatura crolla a 0,237, mentre senza sbiancatura resta a 0,980. La causa è quella annunciata: le ultime componenti hanno varianza nulla o quasi (ricordi i 3 pixel sempre zero) e dividerle per la loro deviazione molto piccola le trasforma in numeri enormi, che dominano la distanza euclidea. Nessun errore, nessun avviso: il codice gira e restituisce un modello pessimo. La regola pratica è che la sbiancatura va usata solo insieme a una scelta di k che escluda le componenti trascurabili, o con un termine di regolarizzazione nella divisione; e va sempre confrontata con la stessa pipeline senza sbiancatura.

Quante componenti per davvero: varianza spiegata contro accuratezza

Il numero di componenti da tenere si sceglie di solito con una soglia di varianza spiegata. Confrontiamo quella soglia con l'accuratezza di un classificatore, per vedere se le due misure dicono la stessa cosa.

print("k    varianza spiegata   accuratezza CV (regressione logistica, 5 fold)")
for k in (2, 5, 10, 20, 30, 64):
    scores = cross_val_score(make_pipeline(PCA(k), LogisticRegression(max_iter=10000)), X, y, cv=5)
    print("%-4d  %-18.3f  %.3f +- %.3f" % (k, pca.explained_variance_ratio_[:k].sum(), scores.mean(), scores.std()))

k    varianza spiegata   accuratezza CV (regressione logistica, 5 fold)
2     0.285               0.574 +- 0.015
5     0.545               0.823 +- 0.040
10    0.738               0.888 +- 0.033
20    0.894               0.895 +- 0.029
30    0.959               0.910 +- 0.030
64    1.000               0.913 +- 0.024

Per ciascun k il listato mette la PCA dentro una pipeline con la regressione logistica e misura l'accuratezza con la validazione incrociata a 5 fold, che rifà la PCA su ogni fold di training (così non c'è fuga di informazione dal fold di validazione). Confronta la varianza spiegata cumulata (quanta informazione geometrica tieni) con l'accuratezza (quanta ne serve davvero al classificatore). Con k = 2 si tiene il 28,5% della varianza e l'accuratezza è 0,574; con k = 5 il 54,5% e 0,823; con k = 10 il 73,8% e 0,888; con k = 20 l'89,4% e 0,895; con k = 30 il 95,9% e 0,910; con tutte le 64 il 100% e 0,913.

L'accuratezza si appiattisce molto prima della varianza: tra 10 e 64 componenti guadagna 2,5 punti, mentre la deviazione standard tra i fold è 0,024–0,033, quindi le differenze tra k = 20, 30 e 64 sono dello stesso ordine del rumore di misura. La conseguenza è una regola operativa: invece di fissare «il 95% di varianza» a priori, scegli il k più piccolo la cui accuratezza non sia distinguibile dal massimo entro la variabilità tra i fold. Qui sarebbe un valore tra 20 e 30, che dimezza o più la dimensione dell'input. Tieni presente che questo risultato riguarda le cifre 8×8, un problema molto facile: su dati diversi la curva avrà un'altra forma.

Tre prove da fare

Esercizio 1. Ripeti l'esperimento di eliminazione del rumore con rumore di deviazione 2 e di deviazione 8: il k ottimale cambia? Perché ti aspetti che si sposti in quella direzione? Esercizio 2. Sbianca solo le componenti con varianza spiegata superiore a una soglia e lascia le altre invariate: recuperi i vantaggi della sbiancatura senza il crollo del KNN? Esercizio 3. Scegli k come il più piccolo valore per cui l'accuratezza CV resta entro una deviazione standard del massimo, e confrontalo con la regola del 95% di varianza spiegata: coincidono?

Per allenarti sulle misure antropometriche per le calotte dei caschi c'è uno scheletro di esercizi in Python con tre parti da completare: autovalori e autovettori di una matrice 2x2 in forma chiusa, varianza spiegata ed errore di ricostruzione, e le probabilità condizionate di t-SNE che non cambiano quando si triplicano tutte le distanze. Finché i punti da completare sono vuoti, stampa quali esercizi mancano.

Che cosa sai fare adesso

Sai derivare la PCA come problema agli autovalori e verificare che massima varianza e minimo errore di ricostruzione coincidono. Sai scegliere quante componenti tenere guardando l'accuratezza e non solo la varianza, usare la PCA per togliere rumore e riconoscere il crollo della sbiancatura. Sai scrivere t-SNE a mano e leggerla senza attribuire significato a distanze e dimensioni dei gruppi. Nella lezione successiva, la regressione lineare, ricominciamo da una retta.

Fonti

Lezione precedente del percorso: Clustering: K-Means e DBSCAN senza etichette. Nella stessa serie rimandano a questa lezione anche: Algebra lineare 2: autovalori, SVD e PCA, KNN: classificare guardando i vicini (e la maledizione della dimensionalità).