09 - Principal Component Analysis (PCA) e riduzione dimensionalità

Un'immagine 8x8 di una cifra scritta a mano ha 64 numeri, un campione di vino ne ha 13, un profilo di espressione genica ne ha migliaia, e quasi sempre molti di quei numeri si muovono insieme: due pixel vicini sono accesi insieme, due misure chimiche crescono insieme. La Principal Component Analysis (PCA) cerca nuovi assi, combinazioni lineari delle variabili originali, ordinati per quanta varianza dei dati (la dispersione attorno alla media) catturano e incorrelati fra loro.

Tenere i primi e scartare gli altri comprime i dati perdendo, per costruzione, la minor varianza possibile. Due errori la rendono inutile o ingannevole: usarla su variabili con unità di misura diverse senza standardizzare, e scambiare la varianza spiegata per l'informazione utile al tuo problema. Questo articolo, il nono della serie ML Fundamentals, la calcola a mano su cinque punti, la riproduce con la SVD e con scikit-learn, e misura i due errori su dati veri.

Ti servono i concetti del clustering (la documentazione di scikit-learn suggerisce proprio la PCA prima del K-Means quando le dimensioni sono molte, perché le distanze euclidee in alta dimensione si gonfiano) e un po' di algebra lineare: autovalori e autovettori di una matrice simmetrica. Tutti i dati sono sintetici con seme fisso o sono dataset inclusi in scikit-learn; ogni numero riportato viene da un listato eseguito.

La PCA a mano su cinque punti

Prendi cinque punti nel piano già centrati sulla media, $(-2,-1)$, $(-1,-1)$, $(0,0)$, $(1,1)$, $(2,1)$. La matrice di covarianza mette sulla diagonale le varianze delle due variabili e fuori diagonale la loro covarianza: $S = X^\top X/(n-1)$. Qui $S_{11} = 10/4 = 2{,}5$, $S_{22} = 4/4 = 1$ e $S_{12} = 6/4 = 1{,}5$.

Gli autovettori di $S$ sono le direzioni delle componenti principali e i corrispondenti autovalori sono le varianze lungo quelle direzioni. Per una matrice 2x2 si ricavano dall'equazione $\lambda^2 - \mathrm{tr}(S)\,\lambda + \det(S) = 0$, cioè $\lambda^2 - 3{,}5\,\lambda + 0{,}25 = 0$, con soluzioni $\lambda = (3{,}5 \pm \sqrt{11{,}25})/2 = 3{,}4271$ e $0{,}0729$.

Il listato rifà il conto in tre modi: con gli autovalori della covarianza, con la SVD (decomposizione ai valori singolari) dei dati centrati, dove gli autovalori sono $\sigma_i^2/(n-1)$, e con PCA di scikit-learn. Poi verifica le tre proprietà che definiscono la PCA: la varianza lungo la prima componente è il primo autovalore, le coordinate nelle nuove assi sono incorrelate, e la varianza totale non cambia, perché una rotazione non crea né distrugge dispersione.

# PCA su cinque punti: covarianza, autovalori, SVD e sklearn danno la stessa risposta
import numpy as np
from sklearn.decomposition import PCA

X = np.array([[-2, -1], [-1, -1], [0, 0], [1, 1], [2, 1]], dtype=float)
n = len(X)
print("medie:", X.mean(axis=0).tolist(), "(i dati sono gia' centrati)")

S = X.T @ X / (n - 1)                                   # matrice di covarianza, divisa per n-1
print("covarianza:\n", S)
print("traccia =", np.trace(S), " determinante =", round(float(np.linalg.det(S)), 4))

lam, V = np.linalg.eigh(S)                              # autovalori in ordine crescente
lam, V = lam[::-1], V[:, ::-1]
print("autovalori:", np.round(lam, 4).tolist())
print("autovettore 1:", np.round(V[:, 0], 4).tolist(), " autovettore 2:", np.round(V[:, 1], 4).tolist())

# stessa cosa con la SVD dei dati centrati: autovalori = sigma^2 / (n-1)
U, s, Vt = np.linalg.svd(X, full_matrices=False)
print("sigma^2/(n-1):", np.round(s ** 2 / (n - 1), 4).tolist())

pca = PCA(n_components=2).fit(X)
print("sklearn explained_variance_:", np.round(pca.explained_variance_, 4).tolist())
print("sklearn explained_variance_ratio_:", np.round(pca.explained_variance_ratio_, 4).tolist())
print("componenti sklearn:", np.round(pca.components_, 4).tolist(), " (il segno puo' cambiare)")

Z = X @ V                                               # coordinate nelle nuove assi
print("varianza lungo la 1a componente:", round(float(Z[:, 0].var(ddof=1)), 4))
print("covarianza fra le due coordinate:", f"{float(np.cov(Z.T)[0, 1]):.1e}")
print("somma delle varianze prima e dopo:", float(np.trace(S)), round(float(Z.var(axis=0, ddof=1).sum()), 4))

# Output:
# medie: [0.0, 0.0] (i dati sono gia' centrati)
# covarianza:
#  [[2.5 1.5]
#  [1.5 1. ]]
# traccia = 3.5  determinante = 0.25
# autovalori: [3.4271, 0.0729]
# autovettore 1: [-0.8507, -0.5257]  autovettore 2: [0.5257, -0.8507]
# sigma^2/(n-1): [3.4271, 0.0729]
# sklearn explained_variance_: [3.4271, 0.0729]
# sklearn explained_variance_ratio_: [0.9792, 0.0208]
# componenti sklearn: [[0.8507, 0.5257], [-0.5257, 0.8507]]  (il segno puo' cambiare)
# varianza lungo la 1a componente: 3.4271
# covarianza fra le due coordinate: 2.6e-16
# somma delle varianze prima e dopo: 3.5 3.5

I tre metodi concordano: autovalori 3,4271 e 0,0729, uguali al conto a mano. La prima componente spiega 0,9792 della varianza totale (3,4271 su 3,5, la traccia della covarianza); la seconda solo 0,0208. La covarianza fra le coordinate nuove vale 2,6e-16, cioè zero a meno dell'arrotondamento, e la somma delle varianze è 3,5 prima e dopo.

Nota il segno dell'autovettore: eigh restituisce $(-0{,}8507, -0{,}5257)$ e scikit-learn $(0{,}8507, 0{,}5257)$. Sono la stessa direzione con verso opposto, e Jolliffe e Cadima avvertono che le coordinate possono cambiare segno con software diversi: non confrontare mai i segni fra due esecuzioni.

Perché la SVD e non la covarianza? La covarianza si forma moltiplicando $X^\top X$, un passo che può perdere precisione, mentre la SVD lavora direttamente sui dati. La documentazione di scikit-learn dice che PCA centra i dati ma non li scala, e che sceglie il solver in base alla forma: la covariance_eigh se ci sono meno di 1000 feature e più di dieci volte i campioni, il randomized (SVD troncata di Halko et al.) se la matrice supera 500x500 e chiedi meno dell'80% delle componenti, altrimenti la SVD completa. Quel «non li scala» è il punto della sezione seguente.

Le unità di misura decidono il risultato

La PCA cerca la varianza massima, e la varianza dipende dall'unità: la stessa lunghezza in metri o in millimetri ha varianza che differisce di un milione di volte. Il dataset del vino incluso in scikit-learn (178 campioni, 13 misure chimiche, 3 cantine) lo mostra bene: la varianza delle 13 feature va da 0,015 a 98.610.

Il listato calcola la PCA sui dati grezzi e su quelli standardizzati, e controlla anche che la PCA sui dati standardizzati sia la PCA sulla matrice di correlazione, la covarianza delle variabili con media 0 e varianza 1. Jolliffe e Cadima dedicano una sezione proprio a questa scelta fra covarianza e correlazione.

# PCA sul dataset vino (178 campioni, 13 misure chimiche): senza e con standardizzazione
import numpy as np
from sklearn.datasets import load_wine
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

wine = load_wine()
X, nomi = wine.data, np.array(wine.feature_names)
print("varianza delle 13 feature: min = %.3f, max = %.1f" % (X.var(axis=0).min(), X.var(axis=0).max()))

for etichetta, M in [("grezzi", X), ("standardizzati", StandardScaler().fit_transform(X))]:
    p = PCA().fit(M)
    r = p.explained_variance_ratio_
    top = nomi[np.argsort(-np.abs(p.components_[0]))[:2]]
    print(f"{etichetta:15s} varianza spiegata PC1 = {r[0]:.3f}, PC2 = {r[1]:.3f}, "
          f"cumulata a 2 = {r[:2].sum():.3f}; PC1 pesa soprattutto: {top.tolist()}")

# PCA sui dati standardizzati = PCA sulla matrice di correlazione
Z = StandardScaler().fit_transform(X)
corr = np.corrcoef(X, rowvar=False)
print("autovalori della correlazione, primi 3:", np.round(np.linalg.eigvalsh(corr)[::-1][:3], 3).tolist())
# StandardScaler divide per n, PCA per n-1: gli autovalori differiscono del fattore (n-1)/n = 177/178
print("explained_variance_ PCA su dati standardizzati (corretti), primi 3:", np.round(PCA().fit(Z).explained_variance_[:3] * 177 / 178, 3).tolist())

# Output:
# varianza delle 13 feature: min = 0.015, max = 98609.6
# grezzi          varianza spiegata PC1 = 0.998, PC2 = 0.002, cumulata a 2 = 1.000; PC1 pesa soprattutto: ['proline', 'magnesium']
# standardizzati  varianza spiegata PC1 = 0.362, PC2 = 0.192, cumulata a 2 = 0.554; PC1 pesa soprattutto: ['flavanoids', 'total_phenols']
# autovalori della correlazione, primi 3: [4.706, 2.497, 1.446]
# explained_variance_ PCA su dati standardizzati (corretti), primi 3: [4.706, 2.497, 1.446]

Sui dati grezzi la prima componente spiega il 99,8% della varianza, ma i pesi più grandi sono quelli di proline e magnesium: la PCA ha scoperto che una feature ha numeri grandi, non che esista una direzione chimica importante. Sui dati standardizzati la prima componente scende al 36,2%, la seconda al 19,2%, e il peso va a flavanoids e total_phenols, due misure legate fra loro.

Gli autovalori della matrice di correlazione (4,706; 2,497; 1,446) coincidono con quelli della PCA sui dati standardizzati, dopo aver corretto il fattore 177/178 che nasce perché StandardScaler divide per $n$ mentre PCA divide per $n-1$. La regola operativa: se le variabili hanno unità diverse, standardizza dentro la Pipeline; se hanno tutte la stessa unità (i pixel di un'immagine), puoi non farlo.

Quante componenti tenere

Gli autovalori ordinati dicono quanta varianza porta ciascuna componente. Il rapporto di varianza spiegata è $\lambda_i / \sum_j \lambda_j$, e il suo cumulo risponde alla domanda «tenendo $k$ componenti, quanta varianza conservo?». La regola popolare del 90% o 95% è una convenzione, non un teorema, e la sua giustificazione è solo la varianza.

Se i dati servono a un classificatore la domanda giusta è un'altra: quanto peggiora il classificatore? Il listato usa le 1.797 immagini 8x8 di load_digits, con i pixel riportati da 0..16 a 0..1 (hanno tutti la stessa unità), e confronta la varianza cumulata con l'accuratezza di una regressione logistica in cross-validation a 5 fold, con la PCA dentro la pipeline e quindi ricalcolata su ogni fold di training.

# Quante componenti tenere: varianza cumulata e accuratezza a valle, su digits (1797 immagini 8x8 = 64 pixel)
import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline

X, y = load_digits(return_X_y=True)
X = X / 16.0                       # pixel da 0..16 a 0..1: stessa unita' per tutti, convergenza piu' rapida
pca = PCA().fit(X)
cum = np.cumsum(pca.explained_variance_ratio_)
print("pixel con varianza zero:", int((X.var(axis=0) == 0).sum()), " -> rango di fatto al massimo", 64 - int((X.var(axis=0) == 0).sum()))
for q in (0.80, 0.90, 0.95, 0.99):
    print(f"componenti per spiegare il {q:.0%} della varianza: {int(np.searchsorted(cum, q) + 1)}")
print("PCA(n_components=0.95) sceglie:", PCA(n_components=0.95).fit(X).n_components_)
print("primi 5 autovalori:", np.round(pca.explained_variance_[:5], 3).tolist(), " autovalori numericamente nulli (< 1e-10):", int((pca.explained_variance_ < 1e-10).sum()))

cv = StratifiedKFold(5, shuffle=True, random_state=0)
print("\n   k  varianza cumulata  accuratezza (CV 5 fold, regressione logistica)")
for k in (2, 5, 10, 20, 29, 40, 64):
    mod = make_pipeline(PCA(n_components=k), LogisticRegression(max_iter=1000))
    acc = cross_val_score(mod, X, y, cv=cv)
    print(f"{k:4d}  {cum[k - 1]:16.3f}  {acc.mean():.3f} (dev. {acc.std():.3f})")
senza = cross_val_score(LogisticRegression(max_iter=1000), X, y, cv=cv)
print(f"senza PCA (64 pixel): {senza.mean():.3f} (dev. {senza.std():.3f})")

# Output:
# pixel con varianza zero: 3  -> rango di fatto al massimo 61
# componenti per spiegare il 80% della varianza: 13
# componenti per spiegare il 90% della varianza: 21
# componenti per spiegare il 95% della varianza: 29
# componenti per spiegare il 99% della varianza: 41
# PCA(n_components=0.95) sceglie: 29
# primi 5 autovalori: [0.699, 0.64, 0.554, 0.395, 0.272]  autovalori numericamente nulli (< 1e-10): 3
#
#    k  varianza cumulata  accuratezza (CV 5 fold, regressione logistica)
#    2             0.285  0.605 (dev. 0.020)
#    5             0.545  0.850 (dev. 0.009)
#   10             0.738  0.928 (dev. 0.010)
#   20             0.894  0.953 (dev. 0.007)
#   29             0.955  0.962 (dev. 0.007)
#   40             0.988  0.968 (dev. 0.008)
#   64             1.000  0.969 (dev. 0.009)
# senza PCA (64 pixel): 0.969 (dev. 0.009)

Servono 21 componenti per il 90% della varianza, 29 per il 95%, 41 per il 99%; tre pixel hanno varianza zero (sono sempre spenti), perciò i tre ultimi autovalori sono nulli e il rango è al massimo 61. PCA(n_components=0.95) sceglie da sola 29. L'accuratezza non segue la varianza: con 2 componenti (28,5% della varianza) è 0,605, con 10 (73,8%) è 0,928, con 29 (95,5%) è 0,962, contro 0,969 senza PCA.

Il guadagno cade fra 10 e 29 componenti e dopo i 40 si ferma. Qui il 95% di varianza costa meno di un punto di accuratezza, ma è una proprietà di questi dati e di questo classificatore, non di tutti: si misura caso per caso, trattando $k$ come iperparametro da scegliere in cross-validation.

Ricostruzione e anomalie

Proiettare su $k$ componenti e tornare indietro con inverse_transform restituisce un'immagine approssimata, e l'errore ha una forma esatta: la somma dei quadrati delle differenze, divisa per la varianza totale, vale $1 - \text{(varianza spiegata)}$. Questo giustifica la lettura della PCA come compressione: la varianza non spiegata è esattamente l'errore di ricostruzione.

Lo stesso errore, calcolato per immagine, può servire da punteggio di anomalia: le immagini che il sottospazio delle cifre non sa descrivere hanno errore alto. Il listato verifica l'identità per tre valori di $k$ e confronta l'errore di cifre vere, non viste in addestramento, con quello di immagini di rumore casuale.

# Ricostruzione: l'errore e' esattamente la varianza che hai scartato; e serve da punteggio di anomalia
import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split

X, y = load_digits(return_X_y=True)
Xtr, Xte = train_test_split(X, test_size=0.3, random_state=0)

print("   k  errore relativo = 1 - varianza spiegata   (misurato sul training)")
for k in (2, 10, 30):
    p = PCA(n_components=k).fit(Xtr)
    R = p.inverse_transform(p.transform(Xtr))
    err = ((Xtr - R) ** 2).sum() / ((Xtr - Xtr.mean(axis=0)) ** 2).sum()
    print(f"{k:4d}  errore {err:.4f}   1 - ratio {1 - p.explained_variance_ratio_.sum():.4f}")

# Punteggio di anomalia: errore di ricostruzione di ogni immagine con k = 20
p = PCA(n_components=20).fit(Xtr)
def errore(M):
    return ((M - p.inverse_transform(p.transform(M))) ** 2).sum(axis=1)

rng = np.random.default_rng(0)
rumore = rng.integers(0, 17, size=(len(Xte), 64)).astype(float)    # pixel casuali 0..16, come le cifre
e_ok, e_rum = errore(Xte), errore(rumore)
auc = roc_auc_score(np.r_[np.zeros(len(e_ok)), np.ones(len(e_rum))], np.r_[e_ok, e_rum])
print(f"errore mediano cifre vere = {np.median(e_ok):.1f}, immagini di rumore = {np.median(e_rum):.1f}")
print(f"AUC dell'errore di ricostruzione come punteggio di anomalia: {auc:.3f}")

# Output:
#    k  errore relativo = 1 - varianza spiegata   (misurato sul training)
#    2  errore 0.7073   1 - ratio 0.7073
#   10  errore 0.2573   1 - ratio 0.2573
#   30  errore 0.0403   1 - ratio 0.0403
# errore mediano cifre vere = 124.3, immagini di rumore = 2322.9
# AUC dell'errore di ricostruzione come punteggio di anomalia: 1.000

L'identità tiene alla quarta cifra decimale: 0,7073 con due componenti, 0,2573 con dieci, 0,0403 con trenta, uguali a $1 - \text{rapporto di varianza}$. Con 20 componenti l'errore mediano delle cifre vere è 124 e quello del rumore 2.323, e l'AUC è 1,000.

Non farti ingannare dalla perfezione: il rumore uniforme è un'anomalia facilissima, perché non ha niente della struttura dei tratti di una cifra. Per anomalie più sottili (una cifra scritta male, una cifra di un'altra scrittura) l'AUC sarebbe più bassa, e va misurata su esempi reali del tuo problema.

Quando la PCA lineare non basta

La PCA trova solo assi dritti. Per due cerchi concentrici, che nessuna retta separa, i due autovalori sono quasi uguali (0,503 e 0,497): non esiste una direzione privilegiata. La Kernel PCA applica lo stesso metodo a una trasformazione implicita delle variabili, con il trucco del kernel visto nella SVM, e secondo la documentazione di scikit-learn ottiene una riduzione non lineare.

Il listato usa il kernel RBF con tre valori di gamma e per ciascuno misura quanto bene una retta (regressione logistica in cross-validation) separa i due cerchi usando le prime $m$ componenti, da 1 a 4.

# Quando la PCA lineare non basta: due cerchi concentrici (cerchio interno a 0.4 del raggio esterno)
import numpy as np
from sklearn.datasets import make_circles
from sklearn.decomposition import PCA, KernelPCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X, y = make_circles(n_samples=400, noise=0.05, factor=0.4, random_state=0)
print("PCA lineare, varianza spiegata dalle due componenti:", np.round(PCA(2).fit(X).explained_variance_ratio_, 3).tolist())

def acc(trasf):
    """Accuratezza di una regressione logistica (cioe' una RETTA) sulle componenti, in CV a 5 fold."""
    return cross_val_score(make_pipeline(trasf, StandardScaler(), LogisticRegression()), X, y, cv=5).mean()

print(f"retta sui dati originali (PCA con 2 componenti): {acc(PCA(2)):.3f}")
print("\nKernelPCA RBF: accuratezza della retta sulle prime m componenti")
print("gamma      m=1    m=2    m=3    m=4")
for g in (1.0, 3.0, 10.0):
    r = [acc(KernelPCA(m, kernel="rbf", gamma=g)) for m in (1, 2, 3, 4)]
    print(f"{g:5.1f}  " + "  ".join(f"{v:6.3f}" for v in r))

# Output:
# PCA lineare, varianza spiegata dalle due componenti: [0.503, 0.497]
# retta sui dati originali (PCA con 2 componenti): 0.435
#
# KernelPCA RBF: accuratezza della retta sulle prime m componenti
# gamma      m=1    m=2    m=3    m=4
#   1.0   0.458   0.450   1.000   1.000
#   3.0   0.492   0.505   1.000   1.000
#  10.0   0.665   0.700   1.000   1.000

La retta sui dati originali dà 0,435, sotto il 50% del caso: non c'è nulla di lineare da imparare. Con la Kernel PCA le prime due componenti non bastano (0,450-0,700 secondo gamma), e con tre componenti l'accuratezza è 1,000 per tutti e tre i valori di gamma. La lezione è doppia: nei metodi a kernel il numero di componenti da tenere non si può dedurre dall'idea che la prima separi tutto, e gamma cambia la qualità delle prime due. Si prova, come per ogni iperparametro.

t-SNE: per disegnare, non per misurare

Per visualizzare dati ad alta dimensione in due dimensioni si usa spesso il t-SNE (van der Maaten e Hinton, 2008), che costruisce una mappa in cui i punti vicini nello spazio originale restano vicini e che, secondo gli autori, riduce la tendenza ad ammassare i punti al centro. A differenza della PCA non è una proiezione lineare, la sua funzione costo non è convessa (con inizializzazioni diverse si ottengono mappe diverse, dice la documentazione) e non ha un metodo transform per punti nuovi.

Per misurare la fedeltà del disegno uso la trustworthiness di scikit-learn, che vale 1 se i dieci vicini di ogni punto nel disegno sono anche vicini nello spazio originale, e la calcolo per la PCA a 2 componenti e per tre esecuzioni del t-SNE con semi diversi su load_digits.

# t-SNE contro PCA per disegnare digits in 2D: trustworthiness (1 = i vicini nel disegno sono vicini davvero)
import time
import numpy as np
from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE, trustworthiness

X, y = load_digits(return_X_y=True)
print("PCA a 2 componenti: trustworthiness =", round(trustworthiness(X, PCA(2).fit_transform(X), n_neighbors=10), 3))
T = []
for seed in (0, 1, 2):
    Z = TSNE(n_components=2, perplexity=30, init="pca", random_state=seed).fit_transform(X)
    T.append(trustworthiness(X, Z, n_neighbors=10))
print("t-SNE a 2 dimensioni, tre semi: trustworthiness =", np.round(T, 3).tolist())
print("TSNE ha un metodo transform per punti nuovi?", hasattr(TSNE(), "transform"))
print("PCA ha un metodo transform per punti nuovi?", hasattr(PCA(), "transform"))

# Output:
# PCA a 2 componenti: trustworthiness = 0.83
# t-SNE a 2 dimensioni, tre semi: trustworthiness = [0.993, 0.993, 0.993]
# TSNE ha un metodo transform per punti nuovi? False
# PCA ha un metodo transform per punti nuovi? True

La PCA a due componenti ha trustworthiness 0,830, il t-SNE 0,993 con tutti e tre i semi: per la sola vicinanza locale, il t-SNE disegna molto meglio. Il prezzo è che le distanze fra gruppi lontani nella mappa non hanno un significato quantitativo, e che i punti nuovi non si possono proiettare, quindi il t-SNE non è un passaggio da mettere in una pipeline di previsione: serve a guardare i dati. La PCA, con transform, resta lo strumento per decorrelare, comprimere e alimentare un modello.

Esercizio svolto

Tre punti sono allineati su una retta: $(-2,-1)$, $(0,0)$, $(2,1)$. Prima di eseguire qualsiasi cosa, calcola la matrice di covarianza, gli autovalori e la direzione della prima componente. Che cosa ti aspetti dalla seconda componente? Poi confronta con il listato.

# Esercizio: tre punti allineati
import numpy as np
from sklearn.decomposition import PCA

X = np.array([[-2.0, -1.0], [0.0, 0.0], [2.0, 1.0]])
S = X.T @ X / (len(X) - 1)
print("covarianza:", S.tolist())
lam, V = np.linalg.eigh(S)
print("autovalori:", np.round(lam[::-1], 6).tolist())
print("direzione principale:", np.round(np.abs(V[:, -1]), 4).tolist(), " 2/sqrt(5), 1/sqrt(5) =", round(2 / 5 ** 0.5, 4), round(1 / 5 ** 0.5, 4))
print("varianza spiegata dalla prima componente:", round(float(PCA(1).fit(X).explained_variance_ratio_[0]), 6))

# Output:
# covarianza: [[4.0, 2.0], [2.0, 1.0]]
# autovalori: [5.0, 0.0]
# direzione principale: [0.8944, 0.4472]  2/sqrt(5), 1/sqrt(5) = 0.8944 0.4472
# varianza spiegata dalla prima componente: 1.0

Che cosa sai fare adesso

Dopo questo articolo sai fare cinque cose, ciascuna con un controllo numerico pronto:

• calcolare la PCA a mano da covarianza, autovalori e autovettori, e riprodurla con la SVD;
• decidere se standardizzare, guardando le unità di misura e non l'abitudine;
• scegliere $k$ confrontando varianza cumulata e accuratezza a valle in cross-validation;
• leggere l'errore di ricostruzione come varianza scartata e usarlo come punteggio di anomalia, con prudenza;
• riconoscere quando serve un metodo non lineare, e quando un grafico t-SNE è solo un disegno.

Le cinque capacità hanno un filo comune: la PCA è una rotazione e una scelta di quali assi buttare, e ogni dubbio su di lei si risolve calcolando quanta varianza, e di quale variabile, stai tenendo. Dopo questo articolo restano tre tappe della serie: le reti neurali feedforward, le metriche di valutazione e il deployment in produzione. Per la PCA scritta da zero e il t-SNE, nella serie Machine Learning Engineer c'è un articolo dedicato.