44 - Overfitting e regolarizzazione: L2, L1 e la trappola dell'accuracy alta

Due studenti preparano un esame. Il primo impara a memoria le risposte dei vecchi compiti, comprese le virgole. Il secondo capisce i principi. Sui vecchi compiti il primo prende il massimo dei voti, il secondo un po' meno; all'esame, con domande nuove, il secondo va molto meglio. Il primo ha fatto overfitting (sovra-adattamento): ha adattato la propria conoscenza al rumore e ai dettagli casuali dei dati visti, e non a ciò che vale in generale. Un modello può fare lo stesso. Il suo errore sui dati di addestramento scende verso zero, e questo sembra una vittoria; ma sull'esame vero (dati mai visti) peggiora.

Questa lezione chiude il blocco sulla regressione con la cura principale per l'overfitting: la regolarizzazione, cioè limitare la libertà del modello penalizzandone la complessità. Vedrai come si ricava, come si implementa a mano, in che cosa differiscono L2 e L1, e la trappola dell'accuratezza che sembra ottima e non lo è.

Le parole da fissare

Il termine che conta di più è generalizzazione: la capacità di un modello di funzionare su dati nuovi. È l'unico obiettivo reale, e l'errore sui dati di addestramento è solo un mezzo per arrivarci. Quanto un modello possa adattarsi si chiama capacità (o complessità): un polinomio di grado 15 ha una capacità enorme, una retta bassa. Se la capacità è troppo poca il modello sbaglia sia sul training sia sul test (underfitting, alto bias); se è troppa va benissimo sul training e male sul test (overfitting, alta varianza). La regolarizzazione è l'insieme delle tecniche che riducono l'overfitting limitando la complessità; qui è un termine aggiunto alla funzione di costo che penalizza i pesi grandi. La sua forza, λ, è un iperparametro: un valore scelto da te prima dell'addestramento, che si seleziona con la validazione e mai guardando il test.

L'idea: un costo che punisce i pesi grandi

Nella lezione sulla regressione polinomiale hai visto che un polinomio di grado 15 raggiunge un errore di training di 0,020 e un errore di test di 3,42, centosettanta volte peggiore. I suoi coefficienti erano nell'ordine dei miliardi: per far passare la curva vicino a ogni punto, il modello usa pesi giganteschi che si compensano a vicenda. Un peso grande significa che una piccola variazione dell'input sposta molto l'uscita: una funzione «nervosa», che insegue il rumore. L'idea è: se il costo include anche la grandezza dei pesi, il modello preferirà, a parità di errore sui dati, una soluzione con pesi piccoli, quindi più liscia. Nella regolarizzazione L2 (o ridge) si aggiunge la somma dei quadrati dei pesi:

J_L2(θ) = (1 / 2m) · Σ_i ( h(x_i) − y_i )²  +  (λ / 2m) · Σ_{j=1..n} θ_j²

Il primo termine è l'errore di sempre. Il secondo è la penalità, e λ ≥ 0 è la forza di regolarizzazione: con λ = 0 non c'è penalità (si torna al modello originale); con λ molto grande i pesi vengono schiacciati verso zero e il modello diventa quasi una costante. Si noti che la somma parte da j = 1: l'intercetta θ₀ non si penalizza, perché non misura la sensibilità a nessuna feature, ma solo il livello medio del target. Penalizzarla farebbe dipendere la soluzione da una costante sommata al target: spostare tutti i valori di y di 100 cambierebbe i pesi delle feature, che non c'entrano nulla.

Cosa cambia nel gradiente e nella soluzione chiusa

Derivare la penalità è facile: la derivata di (λ/2m)θ_j² rispetto a θ_j è (λ/m)θ_j. Il gradiente diventa:

∂J/∂θ_j = (1/m) Σ_i (h(x_i) − y_i) x_ij  +  (λ/m) θ_j         (per j ≥ 1; per j = 0 il secondo termine manca)

Nell'aggiornamento θ_j ← θ_j − α·∂J/∂θ_j il termine extra fa sì che a ogni passo il peso venga prima «ridotto» di un fattore (1 − αλ/m) e poi corretto dal gradiente dell'errore. Per questo L2 si chiama anche weight decay (decadimento dei pesi): in una rete neurale lo troverai con questo nome. Attenzione però: l'equivalenza tra penalità L2 e decadimento dei pesi vale per la discesa del gradiente semplice; per gli ottimizzatori adattivi come Adam le due cose non coincidono, e per questo esiste una variante «a decadimento disaccoppiato» (Loshchilov e Hutter, nelle fonti). La soluzione chiusa si ottiene azzerando il gradiente, come per la regressione lineare:

θ = ( XᵀX + λ·D )⁻¹ Xᵀ y          con D = matrice identità, tranne D₀₀ = 0 (non si penalizza l'intercetta)

Un effetto collaterale prezioso: aggiungere λ sulla diagonale rende XᵀX + λD sempre invertibile e ben condizionata, anche quando le feature sono collineari o più numerose degli esempi. Ridge cura quindi anche l'instabilità dei coefficienti della lezione precedente.

Ridge a mano: il polinomio di grado 15 domato

Ripartiamo dagli stessi dati della lezione sulla regressione polinomiale (sinusoide con rumore, 30 punti di training e 30 di test) e implementiamo la soluzione chiusa con la penalità:

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge

r = np.random.default_rng(1)
x = r.uniform(0, 6, 60)
y = np.sin(x) + r.normal(0, 0.25, 60)
x_tr, x_te, y_tr, y_te = train_test_split(x, y, test_size=0.5, random_state=1)

def potenze(x, grado):
    return np.column_stack([x ** k for k in range(1, grado + 1)])

P_tr = potenze(x_tr, 15)
mu, sd = P_tr.mean(axis=0), P_tr.std(axis=0)                              # standardizzazione dal training
A = np.column_stack([np.ones(30), (P_tr - mu) / sd])                      # (30, 16)
B = np.column_stack([np.ones(30), (potenze(x_te, 15) - mu) / sd])         # (30, 16)

def ridge(A, y, lam):
    D = np.eye(A.shape[1]) * lam
    D[0, 0] = 0.0                                                         # l'intercetta non si penalizza
    return np.linalg.solve(A.T @ A + D, A.T @ y)

print("lambda        MSE train   MSE test   |theta| massimo")
for lam in (0, 1e-4, 1e-2, 1e-1, 1, 10, 100):
    theta = ridge(A, y_tr, lam) if lam > 0 else np.linalg.lstsq(A, y_tr, rcond=None)[0]
    print("%-12g  %9.4f  %9.4f  %12.2f" % (lam, ((A @ theta - y_tr) ** 2).mean(), ((B @ theta - y_te) ** 2).mean(), np.abs(theta[1:]).max()))

sk = Ridge(alpha=1.0).fit(A[:, 1:], y_tr)
print("\nconfronto a lambda = 1: primi coefficienti a mano %s  sklearn %s" % (np.round(ridge(A, y_tr, 1.0)[:3], 5), np.round(np.r_[sk.intercept_, sk.coef_[:2]], 5)))

lambda        MSE train   MSE test   |theta| massimo
0                0.0196     3.4167  3403853961.78
0.0001           0.0259     0.0386         10.96
0.01             0.0328     0.0418          3.40
0.1              0.0488     0.0670          1.46
1                0.0920     0.1292          0.46
10               0.1246     0.1621          0.21
100              0.2562     0.3198          0.08

confronto a lambda = 1: primi coefficienti a mano [ 0.05268  0.01684 -0.45667]  sklearn [ 0.05268  0.01684 -0.45667]

Cosa fa il codice. potenze costruisce le colonne x, x², ..., x¹⁵; vengono standardizzate (senza questo passo i pesi di potenze diverse avrebbero scale incomparabili e la penalità colpirebbe in modo arbitrario, un punto che riprendiamo tra le trappole). ridge è la formula chiusa; D[0, 0] = 0 lascia libera l'intercetta. Il ciclo scorre sette valori di λ e stampa errore di training, errore di test e il coefficiente più grande in valore assoluto. Quattro letture:

• Con λ = 0 (nessuna regolarizzazione) si ritrova il disastro: test 3,42 e coefficienti dell'ordine dei miliardi.
• Già con λ = 10⁻⁴ il test crolla a 0,039 e il coefficiente massimo scende a 11: una penalità minuscola basta a sanare il modello. Gli errori di training e di test si avvicinano.
• Oltre il valore ottimo i due errori risalgono insieme: con λ = 100 il modello è troppo rigido (test 0,32, training 0,26). Questo è l'underfitting da regolarizzazione eccessiva.
• L'errore sul training cresce sempre con λ (da 0,020 a 0,256): la regolarizzazione sacrifica un po' di adattamento ai dati visti per guadagnare generalizzazione. L'errore di test ha forma a «U» (il disegno qui sopra): λ è un compromesso, e si sceglie guardando la validazione.

Le quattro letture, prese insieme, descrivono una curva e non un valore da imparare a memoria: la penalità sposta il modello lungo l'asse che va dall'overfitting (λ nullo, test 3,42) all'underfitting (λ = 100, test 0,32 ma con training 0,26 più alto di quello di ogni altro punto della griglia). Il punto utile sta in mezzo, e nei nostri numeri il test minimo (0,039) cade con λ = 10⁻⁴, il valore più piccolo della griglia dopo lo zero: significa che qui basta una penalità minuscola, ma anche che la griglia non è fine abbastanza per dire dove sia l'ottimo esatto. Il test usato come guida in questo esempio serve solo a mostrare la forma; nel lavoro vero lo stesso confronto si fa su un insieme di validazione.

La riga finale conferma che la formula chiusa coincide con sklearn.linear_model.Ridge (dove il parametro alpha è il nostro λ, nella convenzione in cui la penalità si somma alla somma degli errori quadratici, non alla media): primi coefficienti identici.

L1 (lasso): la regolarizzazione che spegne le feature

La regolarizzazione L1 (lasso) penalizza la somma dei valori assoluti dei pesi invece dei loro quadrati:

J_L1(θ) = (1 / 2m) · Σ_i ( h(x_i) − y_i )²  +  (λ / m) · Σ_j |θ_j|

La differenza sembra piccola ma ha un effetto qualitativo: L1 porta molti pesi esattamente a zero, cioè seleziona le feature. L'intuizione geometrica: con L2 il vincolo equivalente «‖θ‖₂ ≤ c» è un cerchio (una sfera), liscio, e la soluzione tocca il bordo in un punto generico, con tutte le componenti non nulle. Con L1 il vincolo «‖θ‖₁ ≤ c» è un rombo con gli spigoli sugli assi, e le curve di livello dell'errore tendono a toccarlo proprio negli spigoli, dove alcune componenti sono zero.

Un'altra lettura, algebrica: la derivata di |θ| è ±1 anche per θ minuscolo, quindi L1 spinge con forza costante verso lo zero fino ad annullare il peso; la derivata di θ² è 2θ, che sparisce quando θ si avvicina a zero, quindi L2 rimpicciolisce senza mai azzerare. Il valore assoluto non è derivabile in zero: si usano sottogradienti o algoritmi dedicati (coordinate descent in scikit-learn).

Mettiamo i tre approcci a confronto su un problema progettato apposta: 20 feature, ma solo le prime 5 contano davvero (coefficienti veri 3, −2, 1,5, 1, −1), le altre 15 sono rumore, e solo 50 esempi per addestrare.

from sklearn.linear_model import LinearRegression, Lasso

rng = np.random.default_rng(0)
n, p = 100, 20
X = rng.normal(size=(n, p))
beta = np.zeros(p); beta[:5] = [3, -2, 1.5, 1, -1]                        # solo 5 feature contano
yy = X @ beta + rng.normal(0, 1, n)
Xs, Xt, ys, yt = train_test_split(X, yy, test_size=0.5, random_state=0)

print("modello          MSE test   coefficienti esattamente 0 (su 15 di rumore)   |coef| massimo tra le 15 di rumore")
for nome, m in (("minimi quadrati", LinearRegression()), ("ridge  (alpha=1)", Ridge(alpha=1)),
                ("lasso  (alpha=0.05)", Lasso(alpha=0.05)), ("lasso  (alpha=0.1)", Lasso(alpha=0.1)), ("lasso  (alpha=0.5)", Lasso(alpha=0.5))):
    m.fit(Xs, ys)
    c = m.coef_
    print("%-20s %7.3f   %26d                       %10.3f" % (nome, ((m.predict(Xt) - yt) ** 2).mean(), (np.abs(c[5:]) < 1e-8).sum(), np.abs(c[5:]).max()))
print("coefficienti veri delle prime 5:  ", beta[:5])
print("lasso (alpha=0.1), prime 5:       ", np.round(Lasso(alpha=0.1).fit(Xs, ys).coef_[:5], 2))

modello          MSE test   coefficienti esattamente 0 (su 15 di rumore)   |coef| massimo tra le 15 di rumore
minimi quadrati        1.705                            0                            0.391
ridge  (alpha=1)       1.772                            0                            0.408
lasso  (alpha=0.05)    1.416                            5                            0.230
lasso  (alpha=0.1)     1.367                            8                            0.204
lasso  (alpha=0.5)     2.875                           14                            0.013
coefficienti veri delle prime 5:   [ 3.  -2.   1.5  1.  -1. ]
lasso (alpha=0.1), prime 5:        [ 2.81 -1.8   1.13  1.1  -0.82]

Cosa fa il codice. Genera un problema sparso (poche feature realmente rilevanti), divide metà e metà tra training e test e confronta cinque modelli. Il lasso con alpha = 0,1 azzera 8 delle 15 feature di rumore e ottiene il test migliore (1,37 contro 1,71 dei minimi quadrati). Aumentando alpha si azzerano più coefficienti (14 su 15 con alpha = 0,5), ma la penalità ridimensiona anche i coefficienti veri e introduce un errore sistematico (bias): l'errore di test risale a 2,88. Il ridge con alpha = 1 non migliora i minimi quadrati (1,77 contro 1,71): ridimensiona tutti i coefficienti, anche quelli veri, e non azzera nulla. È un risultato onesto che smentisce l'idea «la regolarizzazione migliora sempre»: dipende dal problema. Con pochi coefficienti veri e molto rumore, L1 aiuta di più; con molte feature che contribuiscono un poco ciascuna, L2 è di solito preferibile.

La regolarizzazione nella classificazione: C = 1/λ

La regressione logistica si regolarizza allo stesso modo, aggiungendo la penalità alla log-loss. In scikit-learn il parametro è C, l'inverso della forza: C piccolo = regolarizzazione forte, C enorme = nessuna. Qui la situazione tipica è «molte feature, pochi esempi»: 300 esempi, 100 feature di cui solo 5 informative, addestrando su 150.

from sklearn.datasets import make_classification
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

Xc, yc = make_classification(n_samples=300, n_features=100, n_informative=5, n_redundant=0, class_sep=1.2, random_state=0)
Xa, Xb, ya, yb = train_test_split(Xc, yc, test_size=0.5, stratify=yc, random_state=0)

print("C           accuratezza train   accuratezza test   |peso| massimo")
for C in (1e4, 1, 0.1, 0.01, 0.001):
    m = make_pipeline(StandardScaler(), LogisticRegression(C=C, max_iter=5000)).fit(Xa, ya)
    print("%-10g  %16.3f   %16.3f   %12.2f" % (C, m.score(Xa, ya), m.score(Xb, yb), np.abs(m[-1].coef_).max()))

griglia = GridSearchCV(make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000)),
                       {"logisticregression__C": np.logspace(-4, 2, 13)}, cv=5).fit(Xa, ya)
print("\nC scelto con 5-fold CV sul solo training: %g   accuratezza CV %.3f   accuratezza test %.3f" % (
    griglia.best_params_["logisticregression__C"], griglia.best_score_, griglia.score(Xb, yb)))

l1 = make_pipeline(StandardScaler(), LogisticRegression(l1_ratio=1, solver="liblinear", C=0.1)).fit(Xa, ya)
print("L1 con C=0.1: feature con peso non nullo = %d su 100, accuratezza test %.3f" % ((l1[-1].coef_ != 0).sum(), l1.score(Xb, yb)))

C           accuratezza train   accuratezza test   |peso| massimo
10000                  1.000              0.733           7.48
1                      1.000              0.760           2.23
0.1                    0.987              0.747           0.97
0.01                   0.947              0.767           0.29
0.001                  0.927              0.733           0.05

C scelto con 5-fold CV sul solo training: 0.001   accuratezza CV 0.773   accuratezza test 0.733
L1 con C=0.1: feature con peso non nullo = 6 su 100, accuratezza test 0.880

Cosa fa il codice. La Pipeline mette lo StandardScaler prima del classificatore (la penalità dipende dalla scala delle feature, quindi standardizzare è obbligatorio). Nel primo ciclo, con C enorme (nessuna regolarizzazione) il modello raggiunge accuratezza di training 1.000 e test 0,733: l'overfitting nella sua forma classica. Riducendo C il training scende (0,927 con C = 0,001) e il test sale leggermente (0,767 con C = 0,01): il guadagno è modesto, perché con 100 feature e 150 esempi il problema è difficile per qualunque modello L2. Il secondo blocco mostra l'uso corretto: GridSearchCV sceglie C con validazione incrociata sul solo training (qui C = 0,001) e solo dopo si guarda il test. Il risultato è istruttivo: la validazione incrociata sceglie la regolarizzazione più forte della griglia (C = 0,001) e il test dà 0,733, uguale al modello non regolarizzato.

Con 150 esempi di test, un solo esempio vale 0,0067, quindi le differenze tra 0,733 e 0,767 sono 5 esempi: rumore di campionamento. Non c'è un guadagno dimostrabile dalla sola penalità L2 su questo problema. L'ultima riga è la sorpresa: la penalità L1 con C = 0,1 tiene attive soltanto 6 feature su 100, quasi esattamente le 5 informative, e porta il test a 0,880. Quando la verità è sparsa, selezionare le feature conta più che ridimensionarle tutte. Una nota di versione: con scikit-learn 1.8.0 il parametro penalty risulta deprecato (nella firma vale 'deprecated') e la penalità si sceglie con l1_ratio: 0 è L2, 1 è L1. Con versioni precedenti si scrive penalty="l1"; il solver liblinear supporta L1.

La trappola dell'accuracy alta

In questa lezione un numero di accuratezza rassicurante ha nascosto un problema in tre modi diversi. Il primo è il confronto tra training (1,000) e test (0,733): il numero sul training non ha valore informativo, perché un modello abbastanza flessibile lo raggiunge sempre. Il secondo è usare il test per scegliere λ o C: se provi dieci valori e tieni il migliore sul test, quel test è diventato un secondo training set e il punteggio è ottimista; per questo c'è GridSearchCV sul training e il test resta intatto per la valutazione finale. Il terzo riguarda le classi sbilanciate, trattate nella lezione sul processo di un progetto di ML: un overfitting su pochi positivi cambia poco l'accuratezza complessiva, dominata dai negativi, quindi confronta sempre precision e recall tra training e validazione.

Altri modi di combattere l'overfitting

• Più dati: la cura più efficace, quando disponibile. A parità di modello, più esempi rendono più difficile memorizzare il rumore.
• Modelli più semplici: meno feature, grado più basso, alberi meno profondi.
• Early stopping: interrompere la discesa del gradiente quando l'errore di validazione smette di calare. È una forma implicita di regolarizzazione: pesi piccoli all'inizio e crescenti nel tempo.
• Data augmentation, dropout: tecniche specifiche delle reti neurali (le vedrai più avanti nel modulo).

Le quattro voci non sono alternative equivalenti. Le prime due agiscono sul problema (dati e modello), le altre sulla procedura di addestramento; la regolarizzazione esplicita di questa lezione sta in mezzo ed è la più economica da provare, perché è una riga di codice e un iperparametro. Nel dubbio si parte dai dati, che nessuna penalità sostituisce.

Trade-off e trappole

• L2 contro L1: L2 stabilizza e gestisce bene feature correlate (spartisce il peso fra loro); L1 seleziona e produce modelli compatti ma, con feature correlate, ne sceglie una quasi a caso e può cambiare scelta al variare dei dati. La combinazione si chiama elastic net.
• Dimenticare la standardizzazione: la penalità trattata sulle feature grezze colpisce di più le feature con valori piccoli (che richiedono pesi grandi) e di meno quelle con valori grandi: un effetto arbitrario dovuto alle unità di misura.
• Penalizzare l'intercetta (qui evitato).
• Scegliere λ sul test, o sul training: sempre validazione incrociata.
• Confondere alpha, λ e C: in scikit-learn Ridge(alpha) e Lasso(alpha) usano convenzioni di scala diverse tra loro, e C è l'inverso di λ. Controlla la documentazione dell'oggetto che usi.
• Leakage nella validazione: lo scaler va dentro la Pipeline passata a GridSearchCV, altrimenti i fold di validazione hanno già influenzato la standardizzazione.

Le sei voci si riducono a due famiglie. La prima è la scelta della penalità (L1, L2, elastic net) e ha una risposta che dipende dal problema, come mostrano gli esperimenti di questa lezione. La seconda è la disciplina della procedura: standardizzare dentro la pipeline, non toccare l'intercetta, scegliere la forza solo con la validazione. Gli errori più costosi vengono quasi sempre dalla seconda famiglia, perché producono punteggi buoni e falsi senza dare alcun avviso.

Bias e varianza misurati: perché un po' di distorsione conviene

Il motivo per cui la regolarizzazione funziona si può misurare con una simulazione. Fissiamo la verità (20 feature, di cui solo 5 con coefficienti non nulli: 3, 2, 1,5, 1, 0,5), generiamo 2 000 dataset indipendenti di 40 esempi con rumore di deviazione 2, e su ciascuno addestriamo il modello e prevediamo il valore in uno stesso punto fisso x. Poiché conosciamo il valore vero in quel punto, possiamo scomporre l'errore quadratico medio in bias² (quanto la media delle previsioni, sui vari dataset, si discosta dal vero) e varianza (quanto le previsioni oscillano da un dataset all'altro). Vale MSE = bias² + varianza (più il rumore irriducibile, che qui non contiamo perché confrontiamo con il valore senza rumore).

from sklearn.linear_model import Lasso, ElasticNet, lasso_path

rng_r = np.random.default_rng(0)
n_features, n_samples = 20, 40
true_beta = np.array([3, 2, 1.5, 1, 0.5] + [0] * 15, dtype=float)
x_query = rng_r.normal(size=(1, n_features))
true_value = float((x_query @ true_beta)[0])

alphas = (0, 1, 5, 20, 100)                                                 # 0 = minimi quadrati senza penalità
predictions = {a: [] for a in alphas}
for repetition in range(2000):
    gen = np.random.default_rng(100 + repetition)
    X_sim = gen.normal(size=(n_samples, n_features))
    y_sim = X_sim @ true_beta + gen.normal(0, 2, n_samples)
    for a in alphas:
        model = LinearRegression() if a == 0 else Ridge(alpha=a)
        predictions[a].append(float(model.fit(X_sim, y_sim).predict(x_query)[0]))

print("valore vero nel punto: %.3f" % true_value)
print("alpha (ridge)   bias²     varianza   MSE totale")
for a in alphas:
    v = np.array(predictions[a])
    print("%-14g   %.4f    %.4f     %.4f" % (a, (v.mean() - true_value) ** 2, v.var(), ((v - true_value) ** 2).mean()))

valore vero nel punto: 0.911
alpha (ridge)   bias²     varianza   MSE totale
0                0.0016    3.5930     3.5946
1                0.0057    2.9343     2.9399
5                0.0315    1.9987     2.0302
20               0.1411    1.3215     1.4626
100              0.4419    0.6957     1.1376

Cosa fa il codice. Per ciascun valore di alpha addestra 2 000 modelli su dataset diversi e misura, nel punto fisso (valore vero 0,911), il quadrato dello scarto tra la previsione media e il vero (bias²) e la varianza delle previsioni. Senza penalità (alpha = 0) il bias² è quasi nullo (0,0016), come ci si aspetta da uno stimatore non distorto, ma la varianza è 3,593: con 40 esempi e 20 feature le previsioni oscillano molto da un campione all'altro. Aumentando alpha il bias² cresce (0,0057; 0,0315; 0,1411; 0,4419) e la varianza cala molto più in fretta (2,934; 1,999; 1,322; 0,696).

L'MSE totale scende da 3,595 a 1,138: con alpha = 100 l'errore è circa un terzo di quello dei minimi quadrati, pagando un bias² di 0,44. Questo è il compromesso bias–varianza in forma numerica: un modello leggermente sbagliato in media (distorto verso zero) ma stabile batte il modello corretto in media ma instabile. Due avvertenze. Primo, in questo esperimento l'MSE non ha ancora smesso di scendere a alpha = 100, quindi l'ottimo è oltre la griglia testata: l'esercizio 1 propone di cercarlo. Secondo, è un solo punto di previsione: il bias cambia da punto a punto, e nella pratica alpha si sceglie con la validazione incrociata su molti punti, non con una simulazione che richiede di conoscere la verità.

La penalità dipende dalla scala delle feature: standardizza prima

La penalità ridge somma i quadrati dei coefficienti. Ma un coefficiente dipende dall'unità di misura della feature: se esprimi una lunghezza in millimetri invece che in metri, il suo coefficiente diventa mille volte più piccolo, e viene quindi penalizzato molto meno. Il risultato è che la regolarizzazione favorisce arbitrariamente le feature con scala grande. Lo dimostriamo con due feature identiche per importanza (coefficiente vero 1 per entrambe), dove moltiplichiamo la prima per un fattore s.

n_s = 300
f1, f2 = rng_r.normal(size=n_s), rng_r.normal(size=n_s)
target = f1 + f2 + rng_r.normal(0, 0.5, n_s)

print("scala s della prima feature   coefficienti ridge (alpha=100)   effetto in unità originali (coef × s)   dopo standardizzazione")
for s in (1, 10, 100):
    X_scaled = np.column_stack([f1 * s, f2])
    raw = Ridge(alpha=100).fit(X_scaled, target).coef_
    standardized = Ridge(alpha=100).fit(StandardScaler().fit_transform(X_scaled), target).coef_
    print("%-28d   %-31s   %-37s   %s" % (s, np.round(raw, 3), np.round(raw * np.array([s, 1]), 3), np.round(standardized, 3)))

scala s della prima feature   coefficienti ridge (alpha=100)   effetto in unità originali (coef × s)   dopo standardizzazione
1                              [0.781 0.767]                     [0.781 0.767]                           [0.793 0.758]
10                             [0.102 0.762]                     [1.023 0.762]                           [0.793 0.758]
100                            [0.01  0.762]                     [1.026 0.762]                           [0.793 0.758]

Cosa fa il codice. Due feature hanno lo stesso coefficiente vero (1). Con s = 1 ridge le tratta in modo simile (0,781 e 0,767). Moltiplicando la prima per 10 il suo coefficiente scende a 0,102: ha una scala dieci volte maggiore, quindi per produrre lo stesso effetto serve un coefficiente dieci volte più piccolo, e la penalità sul quadrato lo punisce molto meno. Convertito in unità originali (coefficiente × s) l'effetto della prima feature è 1,023, contro 0,762 della seconda: la feature con la scala maggiore è stata ritirata meno, solo per una questione di unità di misura. Con s = 100 il coefficiente è 0,01 e l'effetto 1,026.

Dopo la standardizzazione, invece, le due feature sono trattate in modo identico per ogni s (0,793 e 0,758, i valori non dipendono da s), e il risultato non dipende più dalle unità. La regola è semplice e non negoziabile: ogni modello con una penalità sui coefficienti (ridge, lasso, elastic net, regressione logistica con C) richiede feature standardizzate, e la standardizzazione va calcolata solo sul training set e dentro la pipeline di validazione incrociata (altrimenti fa fuga di informazione, come visto nella prima lezione del modulo). Variante utile: non penalizzare l'intercetta, cosa che scikit-learn già fa.

Feature correlate: lasso sceglie, elastic net divide, ridge distribuisce

Quando più feature portano la stessa informazione, i tre tipi di penalità si comportano in modo diverso. Lasso (L1) tende a sceglierne alcune e a pesarle in modo arbitrario e instabile (con una penalità più forte ne azzera altre). Ridge (L2) distribuisce il peso fra tutte. L'elastic net combina le due penalità, α·(l1_ratio·Σ|w| + (1 − l1_ratio)/2·Σw²), e tende a mantenere insieme i gruppi di feature correlate azzerando comunque quelle inutili. Creiamo quattro copie rumorose dello stesso segnale z (deviazione 0,1) più quattro feature di puro rumore, con bersaglio 4·z + rumore.

z = rng_r.normal(size=200)
X_corr = np.column_stack([z + rng_r.normal(0, 0.1, 200) for _ in range(4)] + [rng_r.normal(size=200) for _ in range(4)])
y_corr = 4 * z + rng_r.normal(0, 1, 200)

print("coefficienti delle 8 feature: le prime 4 sono copie del segnale, le ultime 4 sono rumore")
for name, model in (("lasso (alpha=0.1)", Lasso(alpha=0.1)), ("elastic net (alpha=0.1, l1_ratio=0.5)", ElasticNet(alpha=0.1, l1_ratio=0.5)), ("ridge (alpha=10)", Ridge(alpha=10))):
    coefficients = model.fit(X_corr, y_corr).coef_
    print("%-38s %s     somma delle prime 4: %.2f" % (name, np.round(coefficients, 2), coefficients[:4].sum()))

coefficienti delle 8 feature: le prime 4 sono copie del segnale, le ultime 4 sono rumore
lasso (alpha=0.1)                      [ 1.44  1.67  0.34  0.45 -0.   -0.02  0.   -0.  ]     somma delle prime 4: 3.89
elastic net (alpha=0.1, l1_ratio=0.5)  [ 1.05  1.09  0.87  0.89 -0.03 -0.06  0.   -0.  ]     somma delle prime 4: 3.90
ridge (alpha=10)                       [ 1.05  1.1   0.91  0.88 -0.09 -0.1   0.04 -0.03]     somma delle prime 4: 3.94

Cosa fa il codice. Con quattro copie quasi identiche del segnale e quattro feature di rumore, i tre modelli stimano coefficienti molto diversi sulle prime quattro, ma con somma quasi uguale (3,89; 3,90; 3,94 vs il valore 4 che la verità richiede). Il lasso con alpha = 0,1 pesa le quattro copie in modo molto diseguale: [1,44; 1,67; 0,34; 0,45]. Con un'altra estrazione dei dati le proporzioni sarebbero diverse, perché il modello non ha motivo di preferire una copia all'altra: è l'instabilità che si cita come difetto del lasso con feature correlate. In questo esperimento, con questa alpha, non ne azzera nessuna; con alpha più alta ne azzererebbe alcune.

L'elastic net con l1_ratio = 0,5 le tiene insieme con pesi molto più simili ([1,05; 1,09; 0,87; 0,89]) e ridge le divide quasi equamente ([1,05; 1,10; 0,91; 0,88]). Sulle quattro feature di puro rumore, il lasso e l'elastic net danno coefficienti quasi o esattamente nulli (da −0,06 a 0,00), mentre ridge ne lascia di piccoli ma non nulli (da −0,10 a 0,04): è la differenza fondamentale, perché ridge non seleziona mai, riduce soltanto. Quale scegliere? Se pensi che poche feature contino davvero e vuoi un modello leggibile, il lasso; se molte contribuiscono ciascuna un poco o le feature sono molto correlate, ridge o elastic net; nel dubbio, l'elastic net con l1_ratio scelto in validazione incrociata è il punto di partenza più prudente.

Il cammino di regolarizzazione del lasso: le feature entrano una alla volta

X_path = rng_r.normal(size=(100, 10))
true_path = np.array([3, -2, 1, 0.5, 0, 0, 0, 0, 0, 0.])
y_path = X_path @ true_path + rng_r.normal(0, 1, 100)
alpha_grid, coefficients, _ = lasso_path(X_path, y_path, alphas=[2, 1, 0.5, 0.2, 0.1, 0.02])

print("alpha   feature non nulle   coefficienti delle prime 4 (veri: 3, −2, 1, 0,5)")
for a, column in zip(alpha_grid, coefficients.T):
    print("%-6g  %-18d   %s" % (a, np.count_nonzero(column), np.round(column[:4], 2)))

alpha   feature non nulle   coefficienti delle prime 4 (veri: 3, −2, 1, 0,5)
2       2                    [ 0.63 -0.51  0.    0.  ]
1       2                    [ 1.72 -1.17  0.    0.  ]
0.5     3                    [ 2.3  -1.53  0.23  0.  ]
0.2     4                    [ 2.69 -1.77  0.61  0.24]
0.1     4                    [ 2.83 -1.85  0.74  0.36]
0.02    8                    [ 2.93 -1.91  0.85  0.44]

Cosa fa il codice. Con 100 esempi e 10 feature, di cui solo 4 con coefficiente non nullo (3, −2, 1, 0,5), calcola il lasso per sei valori decrescenti di alpha. Con alpha = 2 entrano solo due feature (le due più forti) con coefficienti molto ridotti (0,63 e −0,51: il lasso comprime tutto verso zero, anche ciò che tiene). Con alpha = 0,5 entra una terza (0,23), con alpha = 0,2 una quarta (0,24) e i coefficienti crescono verso i valori veri: 2,69, −1,77, 0,61, 0,24. Con alpha = 0,1 le quattro feature vere sono ancora le sole attive (2,83, −1,85, 0,74, 0,36): è l'intervallo in cui il lasso fa quello che promette, la selezione corretta.

Con alpha = 0,02 sono diventate otto: sono entrate anche quattro feature di puro rumore (con coefficienti piccoli), mentre i quattro coefficienti veri sono ormai vicini a quelli dei minimi quadrati (2,93; −1,91; 0,85; 0,44). Si nota anche che i coefficienti veri non vengono mai recuperati esattamente (il vero 3 è stimato 2,93 anche con la penalità minima): è il rumore del campione, non della regolarizzazione. Il messaggio operativo è doppio. Primo: la scelta di alpha non è un dettaglio, è l'atto di decidere quante feature tenere, e si fa con la validazione incrociata. Secondo: il lasso seleziona bene solo in un intervallo di alpha; fuori da quell'intervallo toglie segnale (alpha alti) o fa entrare rumore (alpha bassi).

Tre esercizi

Esercizio 1. Estendi la griglia degli alpha della simulazione bias–varianza fino a 1 000 e 10 000: dove si trova il minimo dell'MSE totale in quel punto? Che cosa succede al bias²? Perché non basta guardare un solo punto per scegliere alpha? Esercizio 2. Nell'esperimento delle scale, ripeti con il lasso invece di ridge: quale feature viene azzerata per prima quando s = 100, e perché? Esercizio 3. Nel cammino del lasso usa la validazione incrociata (LassoCV) per scegliere alpha e confronta il numero di feature selezionate con il numero vero (4): il lasso ne seleziona di più o di meno, e come si spiega?

Che cosa sai fare adesso

Sai scrivere il ridge a mano e riconoscere che coincide con quello di scikit-learn, scegliere λ o C con la validazione incrociata senza toccare il test, e decidere tra L2, L1 ed elastic net in base alla struttura del problema. Sai anche misurare il compromesso bias-varianza con una simulazione e sai perché la standardizzazione, dentro la pipeline, non è facoltativa. Con questa lezione hai chiuso il giro sui modelli lineari (la precedente è la regressione logistica); il passo successivo è il neurone, il percettrone e la ragione per cui servono le funzioni di attivazione non lineari.

Nella stessa serie rimandano a questa lezione anche: Ottimizzazione vincolata: moltiplicatori di Lagrange, dualità e condizioni KKT, Statistica bayesiana: prior, posterior, coniugati e perché la regolarizzazione è una credenza, Regressione multivariata, polinomiale e normalizzazione delle feature.