24 - Algebra lineare 2: autovalori, SVD e PCA

Corso Machine Learning Engineer · Modulo 02 · Matematica e statistica · Lezione 9

Obiettivo. Trovare le direzioni speciali di una matrice (autovettori), capire che ogni matrice è «una rotazione, uno stiramento lungo assi, un'altra rotazione» (la decomposizione ai valori singolari, SVD), usarla per comprimere dati e immagini con la migliore approssimazione di rango basso, e arrivare alla PCA (analisi delle componenti principali), derivandola da zero e confrontandola con scikit-learn.

Falso mito da sfatare. «La PCA sceglie le variabili più importanti.» No: non seleziona nulla; costruisce nuove variabili, combinazioni lineari di tutte le originali, ordinate per quanta varianza portano. La differenza cambia come si interpreta il risultato.

Introduzione

Nella lezione precedente abbiamo visto che una matrice è una trasformazione dello spazio. Ora ci chiediamo: ci sono direzioni in cui la trasformazione si comporta in modo semplice? E come si descrive una matrice rettangolare qualunque (per esempio un dataset n × d) in modo da trovare la sua «struttura principale»? Le risposte sono gli autovalori e la SVD, e la loro applicazione più famosa nei dati è la PCA.

Il percorso definisce autovalori e autovettori, passa alle matrici simmetriche e al teorema spettrale, costruisce la SVD e l'approssimazione di rango basso, e arriva alla PCA derivandola da zero. Si chiude dicendo che cosa la PCA può e non può fare. I prerequisiti sono la lezione su vettori e matrici (rango, condizionamento, equazioni normali) e la covarianza; la riduzione di dimensionalità applicata ai modelli è ripresa in PCA da zero e t-SNE.

Autovalori e autovettori

L'analogia: il foglio di gomma

Immagina di stirare un foglio di gomma. La maggior parte dei punti si sposta in una direzione diversa da quella in cui si trovava rispetto al centro (ruota). Ma esistono direzioni speciali lungo le quali i punti si limitano ad allontanarsi o avvicinarsi dal centro, senza ruotare: restano sulla stessa retta. Quelle direzioni sono gli autovettori; il fattore di stiramento lungo ciascuna è il corrispondente autovalore.

Definizione. v ≠ 0 è un autovettore di A con autovalore λ se

A v = λ v.

L'effetto di A su v è soltanto moltiplicarlo per il numero λ (stiramento se |λ| > 1, compressione se |λ| < 1, ribaltamento se λ < 0).

Come si calcolano

A v = λ v equivale a (A − λI) v = 0, dove I è l'identità. Questo sistema omogeneo ha una soluzione v ≠ 0 solo se la matrice A − λI è singolare (lezione precedente: altrimenti l'unica soluzione sarebbe v = 0), cioè se

det(A − λI) = 0 (il polinomio caratteristico).

Per A = [[a, b], [c, d]]: det(A − λI) = (a − λ)(d − λ) − bc = λ² − (a + d)λ + (ad − bc) = 0, quindi

λ² − tr(A)·λ + det(A) = 0

con tr(A) = a + d la traccia (somma degli elementi sulla diagonale). Dalla formula del secondo grado si leggono due fatti: la somma degli autovalori è la traccia e il prodotto è il determinante.

Esempio. A = [[2, 1], [1, 2]]: tr = 4, det = 3, quindi λ² − 4λ + 3 = 0, λ = 3 o λ = 1. Per λ = 3: (A − 3I)v = [[−1, 1], [1, −1]]v = 0 dà v ∝ (1, 1). Per λ = 1: [[1, 1], [1, 1]]v = 0 dà v ∝ (1, −1).

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 2.0]])
autovalori, autovettori = np.linalg.eig(A)                # autovettori in COLONNA
print("autovalori:", autovalori)
print("autovettori (colonne):\n", np.round(autovettori, 4))

for i in range(2):
    v = autovettori[:, i]
    print(f"A v{i+1} = {np.round(A @ v, 4)} | λ{i+1} v{i+1} = {np.round(autovalori[i] * v, 4)}")

print("somma degli autovalori =", autovalori.sum(), "= traccia", np.trace(A))
print("prodotto degli autovalori =", round(float(autovalori.prod()), 6), "= det", round(float(np.linalg.det(A)), 6))

Output ottenuto:

autovalori: [3. 1.]
autovettori (colonne):
 [[ 0.7071 -0.7071]
 [ 0.7071  0.7071]]
A v1 = [2.1213 2.1213] | λ1 v1 = [2.1213 2.1213]
A v2 = [-0.7071  0.7071] | λ2 v2 = [-0.7071  0.7071]
somma degli autovalori = 4.0 = traccia 4.0
prodotto degli autovalori = 3.0 = det 3.0

Cosa fa il codice. np.linalg.eig restituisce due array: gli autovalori e una matrice i cui colonne sono gli autovettori (normalizzati a lunghezza 1). Il ciclo verifica la definizione Av = λv per ciascuna coppia: i due vettori coincidono. Le ultime righe confermano la somma (traccia) e il prodotto (determinante). Attenzione: eig può restituire i segni degli autovettori in un modo o nell'altro (se v è autovettore lo è anche −v).

A cosa servono

Potenze di una matrice. Applicare A molte volte: A·A·…·A v. Se scrivi v = c₁v₁ + c₂v₂ (combinazione di autovettori) allora Aᵏ v = c₁λ₁ᵏ v₁ + c₂λ₂ᵏ v₂: ogni componente cresce (o decade) come la potenza del suo autovalore, e dopo molte applicazioni domina l'autovettore con l'autovalore di modulo maggiore. È il metodo delle potenze.

Catene di Markov. Un sistema che passa tra stati con probabilità fisse (il tempo: sole o pioggia, con le probabilità di passaggio da un giorno al successivo). La distribuzione «a regime» è l'autovettore con autovalore 1 della matrice di transizione.

def metodo_delle_potenze(M, iterazioni=50, seme=0):
    v = np.random.default_rng(seme).normal(size=M.shape[0])
    for _ in range(iterazioni):
        v = M @ v
        v /= np.linalg.norm(v)                 # rinormalizzo per evitare che esploda
    return v, float(v @ M @ v)                 # autovettore e stima dell'autovalore (quoziente di Rayleigh)

v, lam = metodo_delle_potenze(A)
print("metodo delle potenze:", np.round(v, 4), "| λ stimato =", round(lam, 4), "(atteso 3, direzione (1,1)/√2 ≈ 0.7071)")

# tempo: P[i, j] = probabilità di passare dallo stato j allo stato i (colonne sommano a 1). Numeri inventati.
P = np.array([[0.8, 0.4],           # sole domani | sole oggi, pioggia oggi
              [0.2, 0.6]])          # pioggia domani
stato = np.array([0.0, 1.0])        # oggi piove con certezza
for giorno in range(1, 6):
    stato = P @ stato
    print(f"giorno {giorno}: P(sole) = {stato[0]:.4f}, P(pioggia) = {stato[1]:.4f}")

val, vec = np.linalg.eig(P)
regime = vec[:, np.argmax(val.real)].real
print("distribuzione a regime (autovettore di λ = 1, normalizzato a somma 1):", np.round(regime / regime.sum(), 4))

Output ottenuto:

metodo delle potenze: [-0.7071 -0.7071] | λ stimato = 3.0 (atteso 3, direzione (1,1)/√2 ≈ 0.7071)
giorno 1: P(sole) = 0.4000, P(pioggia) = 0.6000
giorno 2: P(sole) = 0.5600, P(pioggia) = 0.4400
giorno 3: P(sole) = 0.6240, P(pioggia) = 0.3760
giorno 4: P(sole) = 0.6496, P(pioggia) = 0.3504
giorno 5: P(sole) = 0.6598, P(pioggia) = 0.3402
distribuzione a regime (autovettore di λ = 1, normalizzato a somma 1): [0.6667 0.3333]

Cosa fa il codice. metodo_delle_potenze parte da un vettore casuale e applica ripetutamente M, rinormalizzando ogni volta (altrimenti i numeri crescerebbero come 3ᵏ); converge all'autovettore dominante (1, 1)/√2. Il rapporto vᵀMv (quoziente di Rayleigh) stima l'autovalore, vicino a 3. La seconda parte è una catena di Markov con probabilità inventate: partendo da «piove», la distribuzione si avvicina in pochi giorni a un valore stabile, che coincide con l'autovettore di λ = 1 della matrice di transizione (2/3 sole, 1/3 pioggia). Questa stessa logica, su una matrice enorme dei link tra pagine web, è alla base dell'algoritmo PageRank (Brin e Page, 1998).

Matrici simmetriche e il teorema spettrale

Una matrice è simmetrica se A = Aᵀ (uguale alla sua trasposta). Sono ovunque nei dati: la matrice di covarianza, XᵀX, le matrici delle distanze. Per loro vale un risultato molto forte, il teorema spettrale:

1. Gli autovalori sono numeri reali.
2. Esiste una base di autovettori a due a due perpendicolari.
3. Quindi A = Q Λ Qᵀ, con Q matrice ortogonale (le colonne sono gli autovettori di lunghezza 1, QᵀQ = I) e Λ diagonale con gli autovalori.

Perché gli autovettori sono perpendicolari (per due autovalori diversi). Siano Av₁ = λ₁v₁, Av₂ = λ₂v₂ con λ₁ ≠ λ₂, e A simmetrica. Calcoliamo (Av₁)·v₂ in due modi: = λ₁ (v₁·v₂) e, per la simmetria, (Av₁)·v₂ = v₁·(Aᵀv₂) = v₁·(Av₂) = λ₂ (v₁·v₂). Quindi (λ₁ − λ₂)(v₁·v₂) = 0 e, poiché λ₁ ≠ λ₂, v₁·v₂ = 0. ∎

Il significato geometrico: una matrice simmetrica stira lo spazio lungo assi perpendicolari (gli assi di un'ellisse, come nella figura), ciascuno con il suo fattore λ.

rng = np.random.default_rng(1)
M = rng.normal(size=(4, 4))
S = M + M.T                                   # una matrice simmetrica qualunque
lam, Q = np.linalg.eigh(S)                    # eigh: specializzata per matrici simmetriche
print("autovalori (reali, ordinati):", np.round(lam, 4))
print("Q^T Q = I ?                  ", np.allclose(Q.T @ Q, np.eye(4)))
print("S = Q Λ Q^T ?                ", np.allclose(Q @ np.diag(lam) @ Q.T, S))
print("traccia = somma autovalori   :", np.isclose(np.trace(S), lam.sum()))

Output ottenuto:

autovalori (reali, ordinati): [-2.1903  0.1975  1.4328  3.3985]
Q^T Q = I ?                   True
S = Q Λ Q^T ?                 True
traccia = somma autovalori   : True

Cosa fa il codice. Costruisce una matrice simmetrica casuale (M + Mᵀ lo è sempre) e ne calcola gli autovalori con eigh, la versione per matrici simmetriche (più veloce, più accurata, restituisce gli autovalori già ordinati in senso crescente e autovettori ortonormali). Verifica le tre affermazioni del teorema: QᵀQ = I, S = QΛQᵀ e la proprietà della traccia.

Perché è importante. La matrice di covarianza dei dati è simmetrica e semidefinita positiva (tutti gli autovalori sono ≥ 0, perché uᵀ C u è una varianza, che non può essere negativa). I suoi autovettori sono le direzioni principali dei dati, e gli autovalori le varianze lungo di esse: è la PCA, tra poco.

La decomposizione ai valori singolari (SVD)

Gli autovalori esistono per matrici quadrate (e non sempre sono reali). Per una matrice m × n qualunque (un dataset!) c'è una decomposizione che esiste sempre:

A = U Σ Vᵀ

dove U (m × m) e V (n × n) sono matrici ortogonali (rotazioni, eventualmente con un ribaltamento) e Σ è «diagonale» (m × n) con i valori singolari σ₁ ≥ σ₂ ≥ … ≥ 0 sulla diagonale. Nel significato geometrico: ogni trasformazione lineare è la composizione di una rotazione (Vᵀ), uno stiramento lungo gli assi con fattori σᵢ (Σ) e un'altra rotazione (U). La SVD è il cerchio che diventa ellisse, in generale.

Da dove viene (costruzione dagli autovalori)

La matrice AᵀA è simmetrica e semidefinita positiva (per ogni x, xᵀAᵀAx = ‖Ax‖² ≥ 0). Per il teorema spettrale ha autovettori ortonormali v₁, …, vₙ con autovalori λᵢ ≥ 0. Definiamo σᵢ = √λᵢ. Allora:

• ‖A vᵢ‖² = vᵢᵀAᵀA vᵢ = λᵢ = σᵢ²: il vettore vᵢ viene stirato di σᵢ.
• I vettori A vᵢ sono perpendicolari tra loro: (A vᵢ)·(A vⱼ) = vᵢᵀAᵀA vⱼ = λⱼ (vᵢ·vⱼ) = 0 per i ≠ j.
• Dunque, normalizzando, uᵢ = A vᵢ / σᵢ (per σᵢ > 0) sono vettori ortonormali, e A vᵢ = σᵢ uᵢ: è esattamente A V = U Σ, cioè A = U Σ Vᵀ. ∎

Gli autovalori di AᵀA sono quindi i quadrati dei valori singolari: σᵢ² = λᵢ(AᵀA). Le tre voci dicono nell'ordine che cosa fa A a ciascun vᵢ (lo allunga di σᵢ), che le immagini restano perpendicolari fra loro e che, normalizzate, sono le colonne di U: la costruzione non richiede nulla oltre il teorema spettrale, ed è quella che il listato seguente ripercorre in numeri.

rng = np.random.default_rng(2)
A = rng.normal(size=(6, 4))                                  # una matrice rettangolare qualunque

U, s, Vt = np.linalg.svd(A, full_matrices=False)             # forma «economica»
print("forme: U", U.shape, "| s", s.shape, "| Vt", Vt.shape)
print("valori singolari:", np.round(s, 4))
print("A = U Σ V^T ?", np.allclose(U @ np.diag(s) @ Vt, A))

# costruzione dagli autovalori di A^T A
lam, V = np.linalg.eigh(A.T @ A)
ordine = np.argsort(lam)[::-1]                               # ordine decrescente
lam, V = lam[ordine], V[:, ordine]
print("sqrt(autovalori di A^T A):", np.round(np.sqrt(np.clip(lam, 0, None)), 4), "= valori singolari")

U_costr = (A @ V) / np.sqrt(lam)                             # u_i = A v_i / σ_i
print("u_i ortonormali:", np.allclose(U_costr.T @ U_costr, np.eye(4)))
print("rango =", int((s > 1e-10).sum()), "| κ(A) = σ_max/σ_min =", round(float(s[0] / s[-1]), 4), "| np.linalg.cond:", round(float(np.linalg.cond(A)), 4))

Output ottenuto:

forme: U (6, 4) | s (4,) | Vt (4, 4)
valori singolari: [3.0945 2.1676 1.4749 0.365 ]
A = U Σ V^T ? True
sqrt(autovalori di A^T A): [3.0945 2.1676 1.4749 0.365 ] = valori singolari
u_i ortonormali: True
rango = 4 | κ(A) = σ_max/σ_min = 8.4786 | np.linalg.cond: 8.4786

Cosa fa il codice. np.linalg.svd(..., full_matrices=False) restituisce U, i valori singolari s (in un vettore, già ordinati in modo decrescente) e Vt (Vᵀ). Si verifica la ricostruzione A = UΣVᵀ. Poi ripercorre la costruzione della dimostrazione: autovalori di AᵀA, radice quadrata (ottenendo i valori singolari), e uᵢ = Avᵢ/σᵢ, che risultano ortonormali. Le ultime righe collegano la SVD ai concetti della lezione precedente: il rango è il numero di valori singolari non nulli, e il numero di condizionamento è σ_max/σ_min (la definizione che avevamo anticipato).

L'approssimazione di rango basso

Scrivendo il prodotto riga per colonna, A = Σᵢ σᵢ uᵢ vᵢᵀ: una somma di matrici di rango 1, ordinate per importanza (σ₁ ≥ σ₂ ≥ …). Tenendo solo i primi k addendi si ottiene

Aₖ = Σᵢ₌₁ᵏ σᵢ uᵢ vᵢᵀ.

Il teorema di Eckart-Young (1936) dice che Aₖ è la migliore approssimazione di rango k di A: tra tutte le matrici di rango al più k, è quella con l'errore minore (nella norma di Frobenius, la radice della somma dei quadrati di tutti gli elementi). L'errore vale

‖A − Aₖ‖_F = √(σ²ₖ₊₁ + σ²ₖ₊₂ + …).

(Lo enunciamo senza dimostrarlo; si veda Strang. Lo verificheremo numericamente.) Lo spazio occupato scende: invece di m·n numeri se ne memorizzano k·(m + n + 1).

Compressione di un'immagine sintetica

Costruiamo un'immagine 64 × 64 in scala di grigi con strutture semplici: un gradiente, un cerchio e delle strisce, più un po' di rumore.

m = n = 64
y, x = np.mgrid[0:m, 0:n]
immagine = 0.5 * (x / n) + 0.3 * (np.hypot(y - 32, x - 32) < 14) + 0.2 * (np.sin(x / 3.0) > 0.5)
immagine = immagine + np.random.default_rng(3).normal(0, 0.02, immagine.shape)

U, s, Vt = np.linalg.svd(immagine, full_matrices=False)
print("primi 8 valori singolari:", np.round(s[:8], 3))
print("energia cumulata dei primi k (somma dei σ²):")
energia = np.cumsum(s ** 2) / (s ** 2).sum()
for k in (1, 2, 3, 5, 10, 20):
    Ak = (U[:, :k] * s[:k]) @ Vt[:k]
    errore = np.linalg.norm(immagine - Ak)
    teoria = np.sqrt((s[k:] ** 2).sum())
    memoria = k * (m + n + 1) / (m * n)
    print(f"k = {k:>2}: energia {energia[k-1]:.4f} | errore {errore:.4f} = teoria {teoria:.4f} | memoria {memoria:6.1%} dell'originale")

Output ottenuto:

primi 8 valori singolari: [26.343  3.873  1.667  0.941  0.721  0.558  0.465  0.412]
energia cumulata dei primi k (somma dei σ²):
k =  1: energia 0.9703 | errore 4.6069 = teoria 4.6069 | memoria   3.1% dell'originale
k =  2: energia 0.9913 | errore 2.4950 = teoria 2.4950 | memoria   6.3% dell'originale
k =  3: energia 0.9952 | errore 1.8563 = teoria 1.8563 | memoria   9.4% dell'originale
k =  5: energia 0.9971 | errore 1.4281 = teoria 1.4281 | memoria  15.7% dell'originale
k = 10: energia 0.9984 | errore 1.0649 = teoria 1.0649 | memoria  31.5% dell'originale
k = 20: energia 0.9992 | errore 0.7609 = teoria 0.7609 | memoria  63.0% dell'originale

Cosa fa il codice. np.mgrid genera le coordinate di ogni pixel; l'immagine è la somma di un gradiente orizzontale, un disco e delle strisce, con rumore piccolo. La SVD ne dà i valori singolari. Il ciclo ricostruisce Aₖ con i primi k termini ((U[:, :k] * s[:k]) @ Vt[:k] è Σ σᵢ uᵢ vᵢᵀ scritto con un prodotto di matrici), misura l'errore e lo confronta con la formula √(Σ_{i>k} σᵢ²) del teorema: coincidono. L'energia cumulata Σᵢ≤ₖ σᵢ²/Σσ² dice quale frazione della «sostanza» dell'immagine è catturata. Con pochissimi termini (k = 5) l'errore è già piccolo e la memoria scende al 15,7% dell'originale (k·(m + n + 1)/(m·n)) e l'energia catturata è il 99,7%: l'immagine è «quasi di rango basso» perché ha strutture semplici. Un'immagine reale ricca di dettagli decade più lentamente.

Dove si usa. Compressione, riduzione del rumore (si scartano i σ piccoli, spesso rumore), sistemi di raccomandazione (la matrice utenti × prodotti approssimata con pochi fattori latenti), analisi semantica latente dei testi.

La PCA

L'analogia: l'ombra giusta

Hai un oggetto tridimensionale e lo proietti su un muro con una luce: l'ombra è una versione 2D. A seconda dell'orientamento del muro, l'ombra conserva più o meno informazione: una matita vista di punta è un puntino, vista di lato è un segmento lungo. La PCA cerca l'orientamento del muro che massimizza la dispersione dell'ombra, cioè conserva più «forma» possibile. Il muro 1D migliore è la direzione di massima varianza (la prima componente principale), il secondo asse è la direzione di massima varianza perpendicolare alla prima, e così via.

La derivazione

Dati n esempi con d caratteristiche, sia X la matrice n × d centrata (a ciascuna colonna è stata sottratta la sua media). La matrice di covarianza è C = XᵀX/(n − 1) (d × d, simmetrica, semidefinita positiva).

La varianza dei dati proiettati sulla direzione u (un versore, ‖u‖ = 1) è la varianza dei numeri Xu: (1/(n−1))·‖Xu‖² = uᵀ C u.

Vogliamo massimizzare uᵀCu con il vincolo uᵀu = 1. Nel punto di massimo il gradiente della funzione deve essere parallelo al gradiente del vincolo (altrimenti potremmo spostarci lungo il vincolo e aumentare la funzione: è il metodo dei moltiplicatori di Lagrange). Il gradiente di uᵀCu è 2Cu (C simmetrica), quello di uᵀu è 2u. Quindi

C u = λ u:

u deve essere un autovettore della matrice di covarianza, e la varianza ottenuta uᵀCu = λ uᵀu = λ è il suo autovalore. Per massimizzarla si prende l'autovettore con l'autovalore più grande. La seconda componente è l'autovettore del secondo autovalore, perpendicolare alla prima per il teorema spettrale; e così via. Le componenti principali sono dunque gli autovettori di C, e le varianze spiegate sono i suoi autovalori.

Il legame con la SVD. Se X = UΣVᵀ, allora XᵀX = VΣ²Vᵀ: gli autovettori di C sono le colonne di V e gli autovalori sono σᵢ²/(n−1). Perciò la PCA si calcola in modo numericamente più stabile con la SVD dei dati centrati, senza formare XᵀX (che raddoppierebbe il condizionamento: lezione precedente). È ciò che fa scikit-learn.

rng = np.random.default_rng(12)
n = 500
dati = rng.normal(size=(n, 2)) @ np.array([[2.2, 1.4], [0.0, 0.6]])         # nuvola allungata in diagonale
Xc = dati - dati.mean(axis=0)                                               # centraggio

# (a) PCA con gli autovalori della covarianza
C = Xc.T @ Xc / (n - 1)
lam, V = np.linalg.eigh(C)
ordine = np.argsort(lam)[::-1]
lam, V = lam[ordine], V[:, ordine]

# (b) PCA con la SVD dei dati centrati
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
var_svd = s ** 2 / (n - 1)

# (c) scikit-learn
from sklearn.decomposition import PCA
pca = PCA(n_components=2).fit(dati)

print("varianze (autovalori di C)  :", np.round(lam, 4))
print("varianze (σ²/(n-1) dalla SVD):", np.round(var_svd, 4))
print("varianze scikit-learn        :", np.round(pca.explained_variance_, 4))
print("quota di varianza spiegata   :", np.round(pca.explained_variance_ratio_, 4))
print("prima componente (eigh)  :", np.round(V[:, 0], 4))
print("prima componente (SVD)   :", np.round(Vt[0], 4))
print("prima componente (sklearn):", np.round(pca.components_[0], 4), "(il segno può essere opposto)")

punteggi = Xc @ V                                                           # dati nelle nuove coordinate
print("varianza dei punteggi per componente:", np.round(punteggi.var(axis=0, ddof=1), 4))
print("covarianza tra le componenti (≈ 0):", round(float(np.cov(punteggi.T)[0, 1]), 10))

Output ottenuto:

varianze (autovalori di C)  : [5.8321 0.2833]
varianze (σ²/(n-1) dalla SVD): [5.8321 0.2833]
varianze scikit-learn        : [5.8321 0.2833]
quota di varianza spiegata   : [0.9537 0.0463]
prima componente (eigh)  : [-0.8276 -0.5613]
prima componente (SVD)   : [-0.8276 -0.5613]
prima componente (sklearn): [0.8276 0.5613] (il segno può essere opposto)
varianza dei punteggi per componente: [5.8321 0.2833]
covarianza tra le componenti (≈ 0): 0.0

Cosa fa il codice. Genera una nuvola di 500 punti allungata in diagonale e la centra (passo indispensabile: la PCA lavora sulla dispersione attorno alla media). Calcola le componenti in tre modi: (a) autovalori e autovettori della covarianza con eigh (ordinati poi dal più grande), (b) SVD dei dati centrati (le varianze sono σ²/(n−1)), (c) sklearn.decomposition.PCA. Le tre varianze coincidono e le componenti coincidono a meno del segno (un autovettore e il suo opposto sono entrambi validi). Infine proietta i dati sulle componenti (Xc @ V): le varianze dei punteggi sono proprio gli autovalori, e la covarianza tra le componenti è praticamente zero: le nuove coordinate sono incorrelate.

Un esempio reale: i fiori di iris

Il dataset dei fiori (150 esempi, 4 misure) è incluso in scikit-learn. Vediamo quante dimensioni servono davvero, con e senza standardizzazione.

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

Xi, yi = load_iris(return_X_y=True)

for etichetta, dati_i in (("senza standardizzare", Xi), ("standardizzato", StandardScaler().fit_transform(Xi))):
    p = PCA().fit(dati_i)
    print(f"{etichetta:>21}: varianza spiegata per componente = {np.round(p.explained_variance_ratio_, 4)} | prime 2: {p.explained_variance_ratio_[:2].sum():.4f}")

p2 = PCA(n_components=2).fit(StandardScaler().fit_transform(Xi))
print("componenti (righe) sulle 4 misure, standardizzate:\n", np.round(p2.components_, 3))

# ricostruzione dalle prime 2 componenti: quanto si perde?
Z = StandardScaler().fit_transform(Xi)
rico = p2.inverse_transform(p2.transform(Z))
print("errore quadratico medio di ricostruzione con 2 componenti:", round(float(((Z - rico) ** 2).mean()), 4))

Output ottenuto:

 senza standardizzare: varianza spiegata per componente = [0.9246 0.0531 0.0171 0.0052] | prime 2: 0.9777
       standardizzato: varianza spiegata per componente = [0.7296 0.2285 0.0367 0.0052] | prime 2: 0.9581
componenti (righe) sulle 4 misure, standardizzate:
 [[ 0.521 -0.269  0.58   0.565]
 [ 0.377  0.923  0.024  0.067]]
errore quadratico medio di ricostruzione con 2 componenti: 0.0419

Cosa fa il codice. Per le quattro misure dei fiori calcola la quota di varianza spiegata da ciascuna componente, sia sui dati grezzi sia su quelli standardizzati (media 0, deviazione standard 1: lezione 5). Le due versioni danno risultati diversi: la PCA dipende dalla scala, perché la varianza di una variabile misurata in unità grandi pesa di più. Poi mostra i pesi (le «ricette») delle prime due componenti sulle quattro misure e misura l'errore di ricostruzione tornando indietro con inverse_transform: con due componenti (invece di quattro) si perde poca informazione.

Come scegliere quante componenti tenere

• Varianza cumulata: tieni le componenti che servono per raggiungere una soglia (per esempio il 95%), e dichiara la soglia.
• Grafico dello «scree»: la varianza spiegata per componente; si cerca il «gomito» in cui cala bruscamente.
• Scopo a valle: se la PCA precede un modello, scegli k dalla validazione incrociata della prestazione finale, dentro la Pipeline (lezione 8).

I tre criteri vanno dal più descrittivo al più operativo. La soglia cumulata e il gomito guardano solo i dati e sono arbitrari, ma comunicabili; la validazione incrociata guarda l'obiettivo e di solito è quella che conta, perché la quota di varianza non dice se l'informazione scartata serviva alla previsione.

Cosa può e non può fare la PCA

• È lineare: cattura strutture che si descrivono con direzioni rette. Su dati disposti lungo una spirale o un cerchio non trova la struttura.
• Dipende dalla scala: standardizza (a meno che le unità abbiano un senso comune).
• Massimizza la varianza, non la «rilevanza»: la direzione che separa le classi può avere poca varianza. Per la classificazione esistono metodi supervisionati.
• Le componenti sono combinazioni di tutte le variabili: sono più difficili da interpretare delle variabili originali.
• Si addestra sul solo training: se calcoli la PCA su tutto il dataset prima di dividere, hai la fuga di dati della lezione 8. Metti PCA in una Pipeline.
• I segni sono arbitrari: non confrontare i segni di componenti calcolate in momenti diversi.

I limiti si riducono a due domande da porsi prima di usare la PCA: i dati hanno davvero una struttura lineare, e la varianza è ciò che interessa a valle? Se la risposta a una delle due è no (spirale, classi separate da direzioni di poca varianza), la PCA comprime bene e informa male, e conviene un metodo non lineare o supervisionato.

Errori tipici (sintomo, causa, rimedio)

• Non centrare i dati. La prima componente punta verso la media invece che lungo la dispersione. Rimedio: Xc = X - X.mean(axis=0) (scikit-learn lo fa da sé).
• Non standardizzare variabili in scale diverse. Rimedio: StandardScaler.
• eig invece di eigh per matrici simmetriche: non sfrutta la simmetria, non ordina gli autovalori e, se la matrice è simmetrica solo a meno di arrotondamenti, può restituirne di complessi. Rimedio: eigh.
• Autovalori leggermente negativi (−1e-17) in una covarianza. Sono errori di arrotondamento: np.clip(lam, 0, None).
• PCA su tutto il dataset prima dello split. Rimedio: Pipeline.
• Interpretare le componenti come «le variabili importanti». Rimedio: guardare i pesi (components_), ricordando che mescolano tutto.
• Confrontare i segni degli autovettori. Rimedio: normalizzare la convenzione (per esempio componente più grande positiva).
• Usare la SVD completa quando basta quella economica su matrici grandi: spreco di memoria. Rimedio: full_matrices=False o algoritmi troncati.

Gli otto errori si ripartono in tre gruppi. Preparazione dei dati: centrare, standardizzare, fare la PCA dopo lo split. Scelta dello strumento: eigh invece di eig, SVD economica. Lettura dei risultati: autovalori negativi per arrotondamento, segni, componenti scambiate per variabili. Il primo gruppo cambia i numeri, gli altri due cambiano le conclusioni.

Come lavorare nella pratica

Parti dal grafico dei dati e della varianza spiegata. Standardizza, centra e solo dopo applica la PCA, dentro una pipeline, calcolandola con la SVD senza formare XᵀX. Verifica sempre la ricostruzione A ≈ UΣVᵀ e l'ortogonalità delle componenti, dichiara la soglia di varianza e il numero di componenti scelto, e per le matrici simmetriche usa eigh.

Il lessico. Un autovettore è una direzione che la matrice non ruota e l'autovalore il fattore di stiramento lungo di essa; la traccia è la somma degli elementi diagonali, uguale alla somma degli autovalori. Il teorema spettrale scrive una matrice simmetrica come QΛQᵀ; i valori singolari sono le radici quadrate degli autovalori di AᵀA e la SVD è A = UΣVᵀ, cioè rotazione, stiramento, rotazione. Una componente principale è un autovettore della matrice di covarianza, la varianza spiegata la quota di varianza totale che porta, il moltiplicatore di Lagrange la tecnica per massimizzare con un vincolo.

Mini-progetto

1. Applica la PCA a 3 nuvole sintetiche: (a) allungata in diagonale, (b) sferica, (c) due cluster lontani. Per ognuna stampa le varianze spiegate e commenta cosa significa per la dimensione effettiva.
2. Con il dataset delle cifre di scikit-learn (load_digits, immagini 8×8, 64 caratteristiche), trova il numero di componenti per spiegare il 90% della varianza, e confronta l'accuratezza di una regressione logistica con 64 caratteristiche e con quelle componenti, usando una Pipeline(StandardScaler, PCA, LogisticRegression) e la validazione incrociata.
3. Comprimi un'immagine a toni di grigio a scelta (o quella sintetica della lezione, ingrandita) con k = 1, 5, 20 e riporta errore e memoria.

Il primo punto allena a leggere la varianza spiegata come misura della dimensione effettiva: ti aspetti una sola componente dominante nel caso (a), tutte uguali nel (b) e una dominante che separa i gruppi nel (c). Il secondo è un confronto onesto, perché la PCA è dentro la validazione incrociata; il terzo riprende la compressione della lezione con un'immagine tua.

Challenge

• Easy. Calcola a mano gli autovalori e gli autovettori di [[4, 1], [2, 3]] e verificali con np.linalg.eig.
• Medium. Dimostra che gli autovalori di Aᵏ sono λᵏ e che una matrice è invertibile se e solo se nessun autovalore è zero (suggerimento: det A = Π λᵢ).
• Hard. Dimostra che la prima componente principale minimizza la somma delle distanze al quadrato dei punti dalla retta che passa per la media: «massimizzare la varianza proiettata» equivale a «minimizzare l'errore di ricostruzione» (suggerimento: Pitagora, per ogni punto ‖x‖² = (u·x)² + ‖x − (u·x)u‖²).
• Expert. Implementa la PCA con il metodo delle potenze e la deflazione (togliere la prima componente e ripetere) e confrontala con eigh.
• Enterprise. Prendi un dataset della tua organizzazione con almeno 15 caratteristiche numeriche. Applica la PCA dentro una pipeline, stima quante componenti servono per il 90% della varianza, e verifica se un modello addestrato sulle componenti ha prestazioni paragonabili a quello sulle variabili originali.

Le prime due sfide sono esercizi sugli autovalori che si verificano con eig; la Hard è la più istruttiva, perché dimostra che massimizzare la varianza proiettata e minimizzare l'errore di ricostruzione sono lo stesso problema. L'Expert ti fa scrivere a mano ciò che PCA nasconde.

Che cosa sai fare adesso

Sai calcolare autovalori e autovettori di una matrice 2×2 con il polinomio caratteristico, dimostrare che quelli di una matrice simmetrica sono perpendicolari, costruire la SVD dagli autovalori di AᵀA e usare il teorema di Eckart-Young per comprimere un'immagine misurando l'errore. Sai derivare la PCA con i moltiplicatori di Lagrange, calcolarla con la SVD dei dati centrati e confrontarla con scikit-learn, e sai quando non fidarti dei suoi risultati. Hai ora tutti i mattoni del modulo di matematica: il passo successivo è usarli in un modello, a partire da Il processo di un progetto di ML.

Nella stessa serie trovi anche: 25 - Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero (la lezione successiva), 26 - Ottimizzazione vincolata: moltiplicatori di Lagrange, dualità e condizioni KKT, 29 - Campionamento e metodi Monte Carlo: stimare quello che non si riesce a calcolare, 30 - Catene di Markov: la proprietà che dimentica il passato, e la distribuzione che non cambia più (che richiamano questa lezione).