Analisi delle componenti principali (PCA)

Con più di tre variabili ($p>3$) i grafici 1D, 2D e 3D di Correlazione e visualizzazione dei datiLa correlazione di Pearson $r=\sum(X_i-\bar X)(Y_i-\bar Y)/\big(\sqrt{\sum(X_i-\bar X)^2}\sqrt{\sum(Y_i-\bar Y)^2}\big)\in[-1,1]$ misura la relazione lineare tra due variabili (covarianza divisa per le deviazioni standard); correlazione non implica causalità. Serve a capire quali variabili contano per il target e a eliminare quelle quasi duplicate (|r| molto alto). Gli indicatori di sintesi non bastano (quartetto di Anscombe, Datasaurus): vanno affiancati ai grafici: istogramma, KDE, box plot, violin plot, heatmap di correlazione, scatter plot e matrice di scatter plot.Correlazione e visualizzazione dei dati → non bastano. La PCA (Principal Component Analysis) trova pochi nuovi assi, combinazioni lineari delle variabili originali, che conservano la maggior parte della variabilità dei dati: con i primi due o tre assi si può disegnare il dataset (compressione con perdita, lossy), e le stesse componenti servono per ridurre il numero di variabili. L'implementazione da zero in NumPy è in Esercizio - PCA da zero sul dataset California Housing. La nota richiede 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 → (media, varianza, standardizzazione), la covarianza (Covarianza e coefficiente di correlazioneCov(X, Y) = E[(X − E X)(Y − E Y)] = E[XY] − E[X]E[Y] misura quanto X e Y variano insieme; è bilineare, Cov(X, X) = Var(X), Var(X + Y) = Var X + Var Y + 2Cov(X, Y); ρ = Cov / (σ_X σ_Y) sta in [−1, 1] e vale ±1 solo per legami lineari. Indipendenti ⇒ non correlate, ma non viceversa (tranne per i vettori gaussiani).Covarianza e coefficiente di correlazione →) e la nozione di autovalori e autovettori (Autovalori e autovettoriUn autovettore è un vettore non nullo che una funzione lineare manda in un suo multiplo; si trovano gli autovalori come radici del polinomio caratteristico det(A − λI) e gli autovettori come nucleo di A − λI. Matrici simili hanno gli stessi autovalori.Autovalori e autovettori →, DiagonalizzazioneUna matrice è diagonalizzabile se è simile a una diagonale, cioè se esiste una base di autovettori: allora A = S D S⁻¹ con gli autovettori nelle colonne di S e gli autovalori in D. Criterio: tutti gli autovalori nel campo e molteplicità geometrica uguale a quella algebrica. Le matrici simmetriche reali hanno autovalori reali.Diagonalizzazione →).

L'idea in due dimensioni

Si abbiano $n=6$ topi e di ciascuno l'espressione di $p=2$ geni (dati delle slide):

(Le slide usano questa tabella con un'inversione dei valori degli ultimi due topi nella prima versione: qui si usa quella delle slide successive.) Con un solo gene si può disegnare un asse, con due un piano $(x,y)$, con tre uno spazio 3D, ma con tre geni già non si «vede» bene.

Passo 1: centrare. Si calcola il centro dei dati con le medie: $\bar x_1=\frac{35}{6}\approx5{,}83$, $\bar x_2=\frac{21{,}8}{6}\approx3{,}63$. Si traslano i dati sottraendo la media di ciascuna variabile, così che il centro coincida con l'origine $(0,0)$. Le posizioni relative dei punti non cambiano. I punti centrati sono:

$(4{,}17;\,2{,}37),\ (5{,}17;\,0{,}37),\ (2{,}17;\,1{,}37),\ (-2{,}83;\,-0{,}63),\ (-3{,}83;\,-0{,}83),\ (-4{,}83;\,-2{,}63)$.

Passo 2: la retta che meglio si adatta. Si cerca la retta passante per l'origine che passa «più vicino possibile» ai punti. Ruotando la retta e proiettando i punti su di essa si possono seguire due strategie:

• minimizzare le distanze $b$ tra i punti e la loro proiezione sulla retta;
• massimizzare le distanze $c$ tra la proiezione di ogni punto e l'origine.

Sono la stessa cosa. Per ogni punto, il punto, l'origine e la proiezione formano un triangolo rettangolo (la proiezione ortogonale su una retta è in Complemento ortogonale e proiezioni ortogonaliL'ortogonale U⊥ di un sottospazio è un sottospazio di dimensione n − dim U, e R^n = U ⊕ U⊥; ogni vettore si scompone in proiezione su U più componente ortogonale; la proiezione è il punto di U più vicino e si calcola con un sistema o con la matrice di proiezione A(AᵀA)⁻¹Aᵀ.Complemento ortogonale e proiezioni ortogonali →) e per il teorema di Pitagora $$a^2=b^2+c^2,$$ dove $a$ è la distanza del punto dall'origine (non cambia, perché i dati sono fermi), $b$ la distanza dalla retta e $c$ la distanza della proiezione dall'origine. Se $a$ è fisso, aumentare $c$ significa diminuire $b$, e viceversa. Sommando su tutti i punti: $\sum a^2$ è costante, quindi minimizzare $\sum b^2$ equivale a massimizzare $\sum c^2$. Si sceglie di massimizzare $\sum c^2$, perché $\sum c^2/(n-1)$ è la varianza dei dati proiettati: la retta cercata è quella lungo cui i punti proiettati sono più sparpagliati.

La retta trovata è la prima componente principale (PC1). Si indica la somma dei quadrati delle distanze delle proiezioni dall'origine con $SS(\text{distances})=d_1^2+d_2^2+\dots+d_6^2$.

Esempio. Per i dati sopra, $\sum a^2=\sum(x_1^2+x_2^2)=110{,}47$; proiettando sulla migliore retta si ottiene $\sum c^2=106{,}42$ e quindi $\sum b^2=4{,}05$, con $106{,}42+4{,}05=110{,}47$. Qualunque altra retta per l'origine ha $\sum c^2<106{,}42$.

Passo 3: la «ricetta». La retta ha una pendenza, per esempio un vettore lungo $4$ in direzione del gene 1 e $1$ nella direzione del gene 2: la PC1 è «4 parti di gene 1 più 1 parte di gene 2», cioè una combinazione lineare delle variabili. Si normalizza il vettore a lunghezza 1: la lunghezza di $(4,1)$ è $\sqrt{4^2+1^2}=\sqrt{17}\approx4{,}12$, e dividendo ogni coefficiente per $4{,}12$ si ottiene $(0{,}97;\,0{,}242)$, lo stesso rapporto 4 a 1. I coefficienti del vettore normalizzato si chiamano loading (loading scores) e dicono quanto ciascuna variabile contribuisce alla componente; il vettore unitario è l'autovettore (o vettore singolare) della PC1. (Nelle slide i valori $0{,}97$ e $0{,}242$, e poi $15$ e $3$, sono numeri illustrativi della figura, non derivano dalla tabella: con i dati veri il loading della PC1 è $(0{,}942;\,0{,}336)$ e le varianze sono $21{,}28$ e $0{,}81$, vedi sotto.)

Passo 4: le altre componenti. In 2D, la seconda componente principale (PC2) è la retta per l'origine perpendicolare a PC1. Nell'esempio delle slide la ricetta è «$-1$ parti di gene 1 e $4$ parti di gene 2», cioè dopo la normalizzazione $(-0{,}242;\,0{,}97)$. In dimensione $p$ si ottengono $p$ componenti, tutte ortogonali tra loro, ciascuna scelta come retta che massimizza $SS$ tra quelle ortogonali alle precedenti. Il loro numero è al più il minore tra $p$ e $n$.

Autovalori, varianza spiegata e scree plot

Definizione (autovalore di una componente). $\lambda_k=\dfrac{SS(\text{distances per PC}_k)}{n-1}$, cioè la varianza dei dati proiettati sulla componente $k$. La radice quadrata di $SS(\text{distances})$ si chiama valore singolare. In forma matriciale $\lambda_k$ è l'autovalore della matrice di covarianza.

La variazione totale è la somma delle varianze delle componenti, $\sum_k\lambda_k$ (è anche la somma delle varianze delle variabili originali: la PCA ridistribuisce la variabilità, non la cambia).

Formula (varianza spiegata). $\text{frazione spiegata dalla PC}_k=\dfrac{\lambda_k}{\sum_{j=1}^p\lambda_j}.$ La frazione cumulata fino alla componente $k$ si calcola sommando.

Esempio (dati dei topi, due geni). La matrice di covarianza dei dati è $S=\begin{pmatrix}18{,}97&6{,}49\\6{,}49&3{,}13\end{pmatrix}$. Gli autovalori sono le radici di $\lambda^2-\operatorname{tr}(S)\lambda+\det S=0$, cioè $\lambda^2-22{,}09\lambda+17{,}23=0$, e valgono $\lambda_1=21{,}28$ e $\lambda_2=0{,}81$. PC1 spiega $21{,}28/22{,}09\approx96{,}3\%$ della varianza, PC2 il $3{,}7\%$. L'autovettore di $\lambda_1$ risolve $(18{,}97-21{,}28)v_1+6{,}49\,v_2=0$, cioè $v_2=0{,}357\,v_1$: normalizzato dà $(0{,}942;\,0{,}336)$.

Scree plot. Il grafico a barre delle frazioni di varianza spiegata da ciascuna componente. Serve per scegliere quante componenti tenere: se le prime due o tre barre sono alte e le altre basse, il disegno con due o tre componenti è una buona approssimazione; se le barre sono tutte simili, due componenti non rappresentano bene i dati. Una regola pratica è tenere il numero minimo $k$ di componenti per cui la frazione cumulata supera una soglia (per esempio il $95\%$).

Esempio (caso con tre geni, slide). Se PC1, PC2 e PC3 spiegano $79\%$, $15\%$ e $6\%$, un grafico 2D con PC1 e PC2 rappresenta il $94\%$ della variazione.

La PCA con il metodo della covarianza

La procedura usata in laboratorio, per una matrice $X$ con $n$ righe e $p$ colonne:

1. Standardizzare (z-score): $X_{\text{scaled}}=\dfrac{X-\mu}{\sigma}$ per ogni colonna (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 →).
2. Matrice di covarianza: $S=\dfrac{1}{n-1}X_c^{T}X_c$, con $X_c$ i dati centrati, di dimensione $p\times p$ (l'elemento $(j,k)$ è la covarianza tra le variabili $j$ e $k$, Covarianza e coefficiente di correlazioneCov(X, Y) = E[(X − E X)(Y − E Y)] = E[XY] − E[X]E[Y] misura quanto X e Y variano insieme; è bilineare, Cov(X, X) = Var(X), Var(X + Y) = Var X + Var Y + 2Cov(X, Y); ρ = Cov / (σ_X σ_Y) sta in [−1, 1] e vale ±1 solo per legami lineari. Indipendenti ⇒ non correlate, ma non viceversa (tranne per i vettori gaussiani).Covarianza e coefficiente di correlazione →; prodotti tra matrici: Operazioni tra matriciLe matrici m×n formano uno spazio vettoriale (somma e prodotto per scalare elemento per elemento); il prodotto righe per colonne corrisponde alla composizione di funzioni lineari, è associativo ma non commutativo; la trasposta scambia righe e colonne e (AB)^T = B^T A^T.Operazioni tra matrici →).
3. Decomposizione in autovalori e autovettori di $S$: autovettori (direzioni ortogonali) e autovalori (varianza lungo ciascuna).
4. Ordinare gli autovettori per autovalore decrescente.
5. Proiettare: le componenti principali dei dati (scores) sono $Z=X_c\,V_k$, dove $V_k$ ha per colonne i primi $k$ autovettori.

Teorema (perché gli autovettori). Tra tutti i vettori unitari $u$ ($u^Tu=1$), la varianza dei dati proiettati $Z=X_cu$ è $\operatorname{Var}(Z)=u^TSu$, ed è massima quando $u$ è l'autovettore di $S$ con autovalore più grande; il massimo vale $\lambda_1$.

Esempio. Per $S=\begin{pmatrix}2&0\\0&1\end{pmatrix}$ e $u=(\cos\theta,\sin\theta)$: $u^TSu=2\cos^2\theta+\sin^2\theta=1+\cos^2\theta$, massimo $2=\lambda_1$ per $\theta=0$ (l'autovettore $(1,0)$).

Perché. Con il vincolo $u^Tu=1$ si usano i moltiplicatori di Lagrange (Massimi e minimi vincolati e moltiplicatori di LagrangeGli estremi di f sul vincolo g = c si cercano per sostituzione, per parametrizzazione o con i moltiplicatori di Lagrange: se f, g sono C¹, P0 è un estremo vincolato e ∇g(P0) ≠ 0, esiste λ con ∇f(P0) = λ∇g(P0) (curva di livello di f tangente al vincolo). Si risolve il sistema ∇f = λ∇g, g = c, si aggiungono i punti del vincolo con ∇g = 0 e si confrontano i valori; se il vincolo è compatto, Weierstrass garantisce massimo e minimo.Massimi e minimi vincolati e moltiplicatori di Lagrange →): si annulla il gradiente di $u^TSu-\lambda(u^Tu-1)$, che dà $2Su-2\lambda u=0$, cioè $Su=\lambda u$ (autovettore); sostituendo, $u^TSu=\lambda\,u^Tu=\lambda$. Quindi il massimo è l'autovalore più grande. La componente successiva si cerca ortogonale alla precedente: gli autovettori di una matrice simmetrica sono ortogonali (teorema spettrale: Teorema spettrale e forme quadraticheUna funzione lineare è simmetrica se f(v)·w = v·f(w); in una base ortonormale ha matrice simmetrica. Teorema spettrale: f è simmetrica se e solo se esiste una base ortonormale di autovettori, cioè A simmetrica ⇔ PᵀAP diagonale con P ortogonale. Applicato alle forme quadratiche, permette di scriverle come somma di quadrati con gli autovalori come coefficienti.Teorema spettrale e forme quadratiche →), quindi la scelta è automatica.

Un'alternativa: la SVD. Se $X_c=U\Sigma V^T$ è la decomposizione ai valori singolari, le colonne di $V$ sono le componenti principali, i valori singolari sono $\sigma_k=\sqrt{SS_k}$ e $\lambda_k=\sigma_k^2/(n-1)$, gli scores sono $X_cV=U\Sigma$. È il modo in cui le librerie calcolano la PCA («sotto il cofano»). Gli autovettori sono definiti a meno del segno: due implementazioni possono dare grafici specchiati, ugualmente corretti.

Standardizzare: quando e perché

Se le variabili hanno scale molto diverse, la PCA sceglie la direzione di quella con varianza numericamente più grande, che dipende solo dall'unità di misura. Con il dataset California Housing (senza MedInc, 7 variabili):

• senza standardizzare, la prima componente spiega il $99{,}98\%$ della varianza e il suo autovettore ha $-1{,}0$ sulla variabile Population (che ha valori dell'ordine di migliaia) e quasi $0$ sulle altre: la PCA sta solo copiando Population;
• standardizzando, gli autovalori sono $2{,}022;\ 1{,}833;\ 1{,}257;\ 1{,}003;\ 0{,}670;\ 0{,}148;\ 0{,}066$ (somma $7$, pari a $p$ perché ogni variabile ha varianza $1$): la frazione spiegata è $28{,}9\%,\ 26{,}2\%,\ 18{,}0\%,\ 14{,}3\%,\ 9{,}6\%,\ 2{,}1\%,\ 0{,}9\%$, cumulata $28{,}9\%,\ 55{,}1\%,\ 73{,}0\%,\ 87{,}4\%,\ \mathbf{96{,}9\%},\ 99{,}1\%,\ 100\%$. Per spiegare almeno il $95\%$ servono $k=5$ componenti.

import numpy as np
Z = (X - X.mean(0)) / X.std(0)                  # 1. standardizza
S = np.cov(Z, rowvar=False)                     # 2. covarianza (p x p)
lam, V = np.linalg.eig(S)                       # 3. autovalori/autovettori
order = np.argsort(lam)[::-1]                   # 4. ordine decrescente
lam, V = lam[order], V[:, order]
k = np.argmax(np.cumsum(lam / lam.sum()) >= 0.95) + 1   # componenti per il 95 %
scores = Z @ V[:, :k]                           # 5. proiezione

np.cov vuole le variabili sulle righe: con le osservazioni sulle righe serve rowvar=False. Con scikit-learn basta PCA(n_components=2).fit_transform(Z).

Limiti della PCA e metodi non lineari

La PCA è lineare e fa in modo che oggetti lontani restino lontani (conserva la struttura globale), ma la struttura locale non è preservata e non tiene conto di eventuali varietà curve (manifold) su cui giacciono i dati. Per la sola visualizzazione si usano metodi non lineari.

t-SNE (t-distributed Stochastic Neighbor Embedding). Idea: punti simili nello spazio originale devono restare vicini nello spazio ridotto (2D). Concentra l'attenzione sulle piccole distanze, quindi conserva la struttura locale.

• Nello spazio originale la similarità tra $x_i$ e $x_j$ è una probabilità gaussiana centrata in $x_i$, $p_{j|i}=\dfrac{\exp(-\lVert x_i-x_j\rVert^2/2\sigma_i^2)}{\sum_{k\ne i}\exp(-\lVert x_i-x_k\rVert^2/2\sigma_i^2)}$, con $\sigma_i$ scelto in modo da avere una perplexity fissata; poi si simmetrizza, $p_{ij}=\dfrac{p_{j|i}+p_{i|j}}{2N}$.
• Nello spazio ridotto si usa una distribuzione t di Student (code pesanti): $q_{ij}=\dfrac{(1+\lVert y_i-y_j\rVert^2)^{-1}}{\sum_{k\ne l}(1+\lVert y_k-y_l\rVert^2)^{-1}}$.
• Le coordinate $y_i$ si ottengono minimizzando con la discesa del gradiente la divergenza di Kullback-Leibler $KL(p\Vert q)=\sum_{i\ne j}p_{ij}\log\dfrac{p_{ij}}{q_{ij}}$: se $p_{ij}$ è grande (vicini) e $q_{ij}$ piccolo la penalità è alta, quindi i vicini restano vicini; se $p_{ij}$ è piccolo e $q_{ij}$ grande la penalità è lieve, quindi punti lontani possono finire vicini.
• La perplexity è un iperparametro (impostato prima dell'addestramento, non appreso dai dati): è una misura continua del «numero effettivo di vicini» di ogni punto e bilancia struttura locale e globale. La scelta è delicata.
• Limiti: non scala bene, non preserva la struttura globale, le distanze tra cluster non hanno significato. Va bene per visualizzare, non per estrarre feature.

UMAP (Uniform Manifold Approximation and Projection). Costruisce un grafo «fuzzy» dei dati ad alta dimensione e ottimizza un grafo in bassa dimensione il più simile possibile. È molto più veloce di t-SNE, conserva meglio la struttura globale (categorie simili restano vicine) ed è usabile anche per il feature engineering. Parametri: n_neighbors (compromesso locale/globale: valori bassi struttura locale) e min_dist (distanza minima tra punti nello spazio ridotto: valori bassi danno agglomerati più compatti). Limiti: scelta degli iperparametri cruciale, dimensioni dei cluster e distanze tra cluster non significative, risultati diversi tra esecuzioni.

from sklearn.manifold import TSNE
X_2d = TSNE(n_components=2, random_state=42).fit_transform(X)   # es. digits: 64 -> 2 dimensioni

Nell'esempio di laboratorio (cifre scritte a mano, $8\times8$, quindi 64 variabili) due variabili qualsiasi non separano le cifre, mentre il t-SNE le raggruppa in cluster ben distinti.

Errori tipici

• Applicare la PCA a variabili con scale diverse senza standardizzare.
• Calcolare media e deviazione standard (e la PCA) sull'intero dataset prima della divisione train/test: i parametri vanno stimati sul training e applicati al test.
• Interpretare le distanze tra cluster in un grafico t-SNE o UMAP come distanze reali.
• Scambiare le componenti principali per variabili originali: sono combinazioni lineari, meno interpretabili.
• Confondere il segno degli autovettori: due grafici specchiati sono entrambi corretti.
• Tenere troppe poche componenti senza guardare lo scree plot.