5.6 Analisi delle componenti principali (PCA)

Consideriamo in questa sezione un problema tipico del caso vettoriale, che si può affrontare con tecniche simili a quelle introdotte sopra (ossia usando medie e varianza campionarie e stima di massima verosimiglianza). Il problema è di “ridurre la dimensionalità” (in inglese dimension reduction) di una variabile \(Y \in \R^d\) (o similmente di un certo numero \(n\) osservazioni di tale variabile), con \(d\) molto grande, ossia introdurre una variabile \(X \in \R^k\), con \(k\) molto più piccolo di \(d\) in modo da “riassumere” l’informazione di \(Y\) in modo efficace. Questo può essere utile per rappresentare graficamente \(Y\) (ad esempio se \(k=2\)) ma anche per velocizzare l’esecuzione di algoritmi che in dimensione alta possono risultare particolarmente lenti.

Esempio 5.5 Per fare un’esempio concreto, la foto di un volto di una persona incontrata per caso può essere presentata come una variabile \(Y\) a valori in uno spazio molto grande (una dimensione per ogni pixel nell’immagine). Tuttavia quando noi osserviamo un volto ne facciamo automaticamente un “riassunto” tramite caratteristiche quali il colore degli occhi, della pelle, dei capelli ecc. La variabile \(X\) contiene un “riassunto” efficace di \(Y\) (per la nostra memoria, però, mentre per una stampante certamente è più utile direttamente la \(Y\)).

Il problema è presente in tantissimi ambiti scientifici e le tecniche per affrontarlo sono molteplici. Una delle tecniche più semplici, ma comunque efficace, è l’analisi delle componenti principali (in inglese principal component analysis, abbreviato PCA). Dal punto di vista astratto la PCA si può spiegare in modo semplice ricordando il procedimento di standardizzazione di un vettore aleatorio. Data \(Y \in \R^d\), la matrice delle covarianze \(\Sigma_Y\) può essere diagonalizzata tramite il teorema spettrale, ossia esiste \(U_Y \in \R^{d\times d}\) ortogonale (\(U^T_Y = U^{-1}_Y\)) tale che \[ U_Y \Sigma_Y U^T_Y = D_Y\] è diagonale (e contiene gli autovalori di \(\Sigma_Y\)). Dal punto di vista delle variabili, questo significa che tramite un cambio di coordinate, ossia definendo \(Y' = U_Y Y\), la matrice delle covarianze di \(Y'\) è diagonale, e quindi le componenti non sono correlate. A questo punto, volendo “riassumere” \(Y\), si definisce \(X \in \R^k\) come la variabile congiunta delle \(k\) coordinate di \(Y'\) che hanno varianza maggiore. Si riassume quindi \(Y\) catturandone il sottospazio di dimensione \(k\) che presenta maggiore “variabilità” dal punto di vista della covarianza. La variabile \(X\) si ottiene proiettando \(Y\) su tale sottospazio, ossia algebricamente mediante una matrice di proiezione ortogonale \(\Pi_Y \in \R^{k \times d}\), e vale quindi \[ X = \Pi_Y Y.\]

L’idea teorica descritta sopra solleva almeno due problemi, uno pratico e uno teorico: 1. spesso si dispone solamente di un certo numero \(n\) di osservazioni \((y_1, \ldots, y_n)\) associate a variabili aleatorie \((Y_1, \ldots, Y_n)\), tutte indipendenti tra loro e con la stessa legge di \(Y\). Come stimare \(\Pi_Y\)? 2. la PCA è una procedura ad-hoc per questo problema oppure si può giustificare mediante le regole del calcolo della probabilità?

Una soluzione per il primo problema è immediata: invece di considerare la matrice delle covarianze teorica \(\Sigma_Y\) (che non è nota), partendo dalle osservazioni \(y = (y_1, \ldots, y_n)\), si procede allo stesso modo partendo però dalla matrice delle covarianze campionarie, \[\Sigma_{y} = \frac 1 n \sum_{i=1}(y_i-\bar{y})(y_i-\bar{y})^T.\] Ricordiamo infatti che è una matrice simmetrica e semi-definita positiva. Il teorema spettrale applicato \(\Sigma_y\) determina allora una matrice ortogonale \(U_y \in \R^{d\times d}\) e una matrice diagonale \(D_y \in \R^{d\times d}\) (contenente gli autovalori di \(\Sigma_y\)) tali che \[ U_y \Sigma_y U_y^T = D_y.\] A questo punto, si definisce \(\Pi_y \in \R^{k \times d}\) come la matrice di proiezione nel sottospazio associato alle direzioni dei \(k\) vettori di \(U_y\) per cui le componenti nella diagonale (le varianze) sono il più grande possibile. Il “riassunto” in questo caso non è una variabile \(X\) ma il vettore delle osservazioni proiettate \(x_i = \Pi_y y_i\) (anche se tipicamente ciò che interessa sono il sottospazio su cui si proietta e gli autovalori associati).

Esempio 5.6 Applichiamo la PCA in R usando la funzione specifica prcomp() (in alternativa, la decomposizione spettrale di una matrice si può ottenere in generale tramite il comando eigen()). Vediamo ad esempio sul dataset Iris (applicandolo solo ai dati numerici, escludendo la colonna della specie) :

iris_PCA <- prcomp(iris[, 1:4])

L’oggetto risultate contiene diverse informazioni sulla PCA, come ad esempio la base di vettori \(U\) (che è una matrice \(4\times 4\) in questo caso),

iris_PCA$rotation

##                      PC1         PC2         PC3        PC4
## Sepal.Length  0.36138659 -0.65658877  0.58202985  0.3154872
## Sepal.Width  -0.08452251 -0.73016143 -0.59791083 -0.3197231
## Petal.Length  0.85667061  0.17337266 -0.07623608 -0.4798390
## Petal.Width   0.35828920  0.07548102 -0.54583143  0.7536574

e le deviazioni standard delle varie componenti (ossia la radice quadrata dei vari autovalori della matrice delle covarianze empirica \(\Sigma_y\), o della sua diagonalizzata \(D_y\)).

iris_PCA$sdev

## [1] 2.0562689 0.4926162 0.2796596 0.1543862

Inoltre vi sono già calcolate tutte le osservazioni trasformate tramite la matrice \(U\), ossia una matrice che contiene i vettori \(Uy_i\) (come righe). Di questa possiamo ad esempio plottare la prima colonna (corrispondente alla direzione con varianza massima), o le prime due. La PCA ha l’effetto in questo caso di ben separare le tre specie, o quanto meno la prima dalle altre due, come evidenziamo con la diversa colorazione.

plot(iris_PCA$x[, 1], col = miei_colori[as.numeric(iris$Species)],
  pch = 16, xlab = "esemplare", ylab = "valore")

plot(iris_PCA$x[, 1:2], col = miei_colori[as.numeric(iris$Species)],
  pch = 16)

Veniamo alla seconda domanda, e mostriamo che è possibile giustificare il metodo di PCA in termini di una stima di massima verosimiglianza per un opportuno modello gaussiano. L’idea è che con la PCA stiamo recuperando un “segnale” (la \(X\)) osservandone una versione “rumorosa” e disposta su un sottospazio non noto.

Per introdurre un modello probabilistico, supponiamo fissata (e nota) la dimensione \(k\) introduciamo una variabile standardizzata \(Z \in \R^k\) e imponiamo che valga \[Y = AZ+W,\] dove \(A \in \R^{d \times k}\) è una matrice non nota (rispetto all’informazione priori). Il “segnale” da ricostruire è quindi \(AZ\) (quello che nella PCA abbiamo chiamato \(X\)) e \(W\) è una variabile che rappresenta il “rumore” aggiunto. Supponiamo che \(Z\), \(W\) siano indipendenti con densità gaussiane centrate e, oltre a \(\Sigma_Z = Id\), supponiamo che \(\Sigma_W = \sigma_0^2 Id\), per una costante opportuna (nota a priori e sufficientemente piccola). Supponendo nota la matrice \(A\), allora la densità di \(Y\), è anch’essa una gaussiana centrata, con covarianza \(\Sigma_Y = AA^T + \sigma_0^2 Id\), chè è una funzione di \(A\). Pertanto la verosimiglianza di \(A\) associata ad \(Y\) si scrive \[L(A; y) = p(Y=y | A) \propto \exp\bra{ - \frac 1 2 \bra{ y^T \sigma_Y^{-1} y + \log(\det(\Sigma_Y)) }}.\] Ovviamente, se invece di una singola osservazione \(Y=y\), supponiamo di avere \(n\) osservazioni indipendenti \(Y_i= y_i\), tutte gaussiane con gli stessi parametri – in particolare con la stessa matrice \(A\), la verosimiglianza si ottiene come prodotto della funzione sopra (cambiando i valori osservati) \[ L(A; y_1, \ldots, y_n) \propto \exp\bra{ - \frac n {2} \bra{ \frac 1 n\sum_{i=1}^n y_i^T \Sigma_Y^{-1} y_i + \log(\det(\Sigma_Y)) }}.\] Per stimare \(A\) si tratta quindi di determinare \(A \in \R^{d \times k}\) tale che la quantità sopra sia massima, ossia (passando al solito all’opposto del logaritmo e minimizzando) \[ A \mapsto \frac 1 n\sum_{i=1}^n y_i ^T \Sigma_Y^{-1} y_i + \log(\det(\Sigma_Y)\] sia minima (ricordiamo che \(\Sigma_Y\) dipende da \(A\)). Riconosciamo quindi una variante del problema della stima dei parametri di una densità gaussiana vettoriale a partire da osservazioni indipendenti, dove stavolta la dipendenza dei parametri (ossia \(A\)) è più complessa. Tuttavia, con qualche passaggio di algebra lineare e calcolo in più variabili si può dedurre che la stima di massima verosimiglianza per \(A\) è data da \[ A_{{\mle}} = U_{y|k} (D_{y|k} - \sigma^2_0Id)^{1/2},\] dove \(U_{y|k} \in \R^{d \times k}\) indica la matrice corrispondente ai \(k\) autovettori della covarianza campionaria \(\Sigma_y = \sum_{i=1}^n y_i y_i^T\) con autovalori più grandi, e \(D_{y|k} \in \R^{k\times k}\) indica la matrice diagonale contenente tali autovalori nell’ordine corrispondente. Tutto questo purché \(\sigma_0^2\) sia sufficientemente piccolo, affinché la diagonale di \(D_{y|k} - \sigma^2_0 Id\) sia positiva, altrimenti la radice quadrata (intesa solamente sulla diagonale) non avrebbe senso. Nel limite \(\sigma_0 \to 0\) si ottiene che \(A_{{\mle}}= U_{y|k} D_{y|k}^{1/2}\) e la variabile \(A_{{\mle}}X\) si identifica con \(\Pi_y Y\) (identificando \(\R^k\) con il sottospazio generato dai \(k\) vettori delle colonne di \(U_{y|k}\)).

5.6.1 Esercizi

Esercizio 5.10 Si consideri il dataset ‘mtcars’ contenente dati relativi ad alcuni modelli di auto (ormai d’epoca). Si usi opportuni comandi R per calcolare e visualizzare la matrice delle correlazioni campionarie e successivamente si applichi la PCA e si plotti un diagramma a nuvola di punti relativo alle prime due componenti principali.

Esercizio 5.11 Si consideri un dataset ricavato dalle visualizzazioni di pagine di Wikipedia (https://pageviews.toolforge.org/) costruito come nell’esercizio della sezione precedente. Si applichi la PCA e si plotti un diagramma a nuvola di punti relativo alle prime due componenti principali.