04 - Bias-Variance Tradeoff: il dilemma centrale del ML

Un modello sbaglia in due modi opposti, e il rimedio per l'uno peggiora l'altro. Se è troppo rigido, una retta per dati che disegnano una curva, sbaglia nello stesso modo su qualunque campione gli dai: errore sistematico, in gergo bias. Se è troppo flessibile, un polinomio di grado 9 su trenta punti, insegue il rumore del campione che ha visto e cambia forma ad ogni campione nuovo: errore di instabilità, in gergo varianza. Questo articolo ti dà il modo di capire quale dei due ti sta rovinando il modello e che cosa fare: ricava la scomposizione dell'errore, la misura con una simulazione in cui la verità è nota, e la usa per leggere le curve di apprendimento, scegliere la regolarizzazione e capire perché mediare più modelli aiuta. Servono il capitolo 02 (minimi quadrati e Ridge) e un'idea di media e varianza.

Tutto è stato eseguito con numpy 2.4.6 e scikit-learn 1.9.0, su dati sintetici con seme fisso e funzione vera dichiarata, quindi senza download. Ogni uscita riportata è identica in due esecuzioni.

Il problema in un esempio

La funzione vera è f(x) = sin(2πx) su [0, 1]; ogni misura è y = f(x) + ε con rumore gaussiano di deviazione standard 0,3, quindi varianza σ² = 0,09. Il rumore è la parte dell'errore che nessun modello può togliere: anche conoscendo f alla perfezione, una nuova misura sbaglierebbe in media di 0,09 al quadrato. Prendiamo trenta punti e adattiamo polinomi di grado 1, 3 e 9. Il grado 1 non riuscirà mai a seguire una sinusoide (sottoadattamento, underfitting); il grado 9 ha dieci coefficienti per trenta punti e li piega sul rumore (sovradattamento, overfitting): sul campione che ha visto sbaglia poco, su uno nuovo molto.

La scomposizione dell'errore

Immagina di ripetere l'esperimento: ogni volta peschi un campione di addestramento D diverso e ottieni un modello f̂_D diverso. In un punto fissato x₀ l'errore quadratico atteso su una nuova misura si scompone in tre termini: E[(y − f̂_D(x₀))²] = Var[f̂_D(x₀)] + (E[f̂_D(x₀)] − f(x₀))² + σ². Il primo è la varianza, quanto cambia la previsione da un campione all'altro; il secondo è il bias al quadrato, la distanza fra la previsione media e la verità; il terzo è il rumore irriducibile. Il manuale di James, Witten, Hastie e Tibshirani (An Introduction to Statistical Learning, paragrafo 2.2.2) riporta la stessa scomposizione. La dimostrazione si fa sommando e sottraendo f(x₀) e E[f̂_D(x₀)] dentro il quadrato: i prodotti misti si annullano perché il rumore della nuova misura è indipendente dal campione, e perché la media degli scarti dalla media è zero.

La scomposizione ha una conseguenza pratica. Con un solo dataset non puoi misurare bias e varianza separatamente, perché ti serve la verità e più campioni. Puoi però riprodurre la situazione con una simulazione dove la verità è dichiarata, e vedere se la scomposizione tiene.

Misurare bias e varianza con una simulazione

Il listato genera 500 campioni di 30 punti, adatta su ciascuno un polinomio, e confronta le 500 curve ottenute in 49 punti della griglia. Il bias al quadrato è la distanza media fra la curva media e f, la varianza è la dispersione media delle curve attorno alla loro media. Come controllo indipendente, l'ultimo calcolo valuta le stesse curve contro etichette fresche e rumorose: se la scomposizione è giusta, il risultato deve coincidere con la somma dei tre termini.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(0)
f = lambda x: np.sin(2 * np.pi * x)                 # funzione vera
sigma, n, M = 0.3, 30, 500                          # rumore, esempi per dataset, numero di dataset
xg = np.linspace(0.02, 0.98, 49)                    # punti in cui si misura l'errore

print(f"rumore irriducibile sigma^2 = {sigma**2:.4f}")
for grado in (1, 3, 9):
    P = np.empty((M, len(xg)))                      # previsioni di M modelli, ciascuno su un dataset diverso
    for m in range(M):
        x = rng.uniform(0, 1, n)
        y = f(x) + sigma * rng.standard_normal(n)
        P[m] = Polynomial.fit(x, y, grado)(xg)
    bias2 = np.mean((P.mean(axis=0) - f(xg)) ** 2)
    var = np.mean(P.var(axis=0))
    y_nuovo = f(xg) + sigma * rng.standard_normal(P.shape)    # etichette fresche, mai viste
    test = np.mean((y_nuovo - P) ** 2)
    print(f"grado {grado}: bias^2={bias2:.4f} varianza={var:.4f} somma con sigma^2={bias2 + var + sigma**2:.4f}  errore di test misurato={test:.4f}")

# Output:
# rumore irriducibile sigma^2 = 0.0900
# grado 1: bias^2=0.1820 varianza=0.0231 somma con sigma^2=0.2951  errore di test misurato=0.2951
# grado 3: bias^2=0.0037 varianza=0.0158 somma con sigma^2=0.1095  errore di test misurato=0.1096
# grado 9: bias^2=0.0220 varianza=8.7130 somma con sigma^2=8.8250  errore di test misurato=8.8381

Il grado 1 ha un bias al quadrato di 0,1820 e una varianza di 0,0231: l'errore è quasi tutto sistematico, e nessun campione in più lo toglie. Il grado 3 ha bias 0,0037 e varianza 0,0158, e l'errore atteso scende a 0,1095, vicino al rumore irriducibile di 0,09. Il grado 9 ha un bias piccolo (0,0220) ma una varianza di 8,7130, e l'errore diventa 8,825. La colonna di controllo conferma la scomposizione: 0,2951 contro 0,2951, 0,1095 contro 0,1096, 8,825 contro 8,838; la piccola differenza dell'ultimo caso è rumore di campionamento, perché con code così pesanti poche curve estreme dominano la media. La varianza così alta del grado 9 non è un difetto del metodo di misura: con 10 coefficienti su 30 punti rumorosi il polinomio è poco vincolato fra i punti.

Le curve di apprendimento

Nella pratica non conosci la verità, ma puoi guardare come cambia l'errore sul training e su dati nuovi al crescere del numero di esempi. La guida di scikit-learn alle curve di validazione e di apprendimento (learning_curve, validation_curve) descrive lo stesso strumento. Qui si calcola con la simulazione: per ogni dimensione n di campione si adattano 200 modelli e si riporta la mediana dell'errore sul training e dell'errore su una griglia di test, rumore incluso. Si usa la mediana perché per il grado 9 la media è dominata da pochi campioni estremi.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(1)
f = lambda x: np.sin(2 * np.pi * x)
sigma, M = 0.3, 200
xt = np.linspace(0, 1, 400)                          # griglia di test densa
print("mediana su 200 dataset; n | grado 1: train  test | grado 9: train  test")
for n in (12, 20, 40, 80, 160, 320):
    riga = []
    for grado in (1, 9):
        tr, te = [], []
        for _ in range(M):
            x = rng.uniform(0, 1, n)
            y = f(x) + sigma * rng.standard_normal(n)
            p = Polynomial.fit(x, y, grado)
            tr.append(np.mean((y - p(x)) ** 2))
            te.append(np.mean((f(xt) - p(xt)) ** 2) + sigma ** 2)
        riga += [np.median(tr), np.median(te)]
    print(f"{n:4d} | {riga[0]:12.4f} {riga[1]:6.4f} | {riga[2]:12.4f} {riga[3]:6.4f}")

# Output:
# mediana su 200 dataset; n | grado 1: train  test | grado 9: train  test
#   12 |       0.2344 0.3182 |       0.0103 998.4158
#   20 |       0.2388 0.3117 |       0.0397 1.2177
#   40 |       0.2614 0.3001 |       0.0659 0.1297
#   80 |       0.2719 0.2935 |       0.0778 0.1057
#  160 |       0.2797 0.2909 |       0.0845 0.0962
#  320 |       0.2834 0.2891 |       0.0871 0.0927

La lettura è la regola pratica della diagnosi. Con il grado 1 le due colonne si avvicinano e si fermano insieme a un livello alto: errore sul training 0,2834 ed errore sul test 0,2891 a n = 320, ben sopra il rumore irriducibile di 0,09. Quando entrambe le curve convergono in alto, il modello è troppo semplice e aggiungere dati non serve: serve un modello più flessibile o altre misure. Con il grado 9 l'errore sul training è bassissimo (0,0103 a n = 12, dove il polinomio quasi interpola i punti) mentre quello di test è enorme (mediana 998,4); il divario si chiude solo al crescere di n e a 320 i due errori sono 0,0871 e 0,0927, entrambi vicini al limite di 0,09. Un grande divario fra training e test indica varianza, e la cura passa dai dati (più esempi), da un modello più semplice, o dalla regolarizzazione.

Esercizio svolto: quanti dati servono per dimezzare la varianza?

Per un modello con un numero fisso di parametri la varianza delle previsioni scala come 1/n: raddoppiare gli esempi la dimezza, non servono quadrupli. Verifichiamolo sul grado 3, con 1000 campioni per ogni n.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(4)
f = lambda x: np.sin(2 * np.pi * x)
sigma, M = 0.3, 1000
xg = np.linspace(0.02, 0.98, 49)
prec = None
for n in (30, 60, 120, 240):
    P = np.empty((M, len(xg)))
    for m in range(M):
        x = rng.uniform(0, 1, n)
        P[m] = Polynomial.fit(x, f(x) + sigma * rng.standard_normal(n), 3)(xg)
    var = np.mean(P.var(axis=0))
    print(f"n={n:3d}  varianza (grado 3) = {var:.5f}" + (f"  rapporto col precedente = {prec / var:.2f}" if prec else ""))
    prec = var

# Output:
# n= 30  varianza (grado 3) = 0.01686
# n= 60  varianza (grado 3) = 0.00657  rapporto col precedente = 2.57
# n=120  varianza (grado 3) = 0.00310  rapporto col precedente = 2.12
# n=240  varianza (grado 3) = 0.00156  rapporto col precedente = 1.99

La varianza va da 0,01686 a 0,00156 quando n passa da 30 a 240, un fattore di circa 11 per un aumento di 8 volte dei dati; il rapporto fra due livelli successivi è 2,57, 2,12 e 1,99, quindi si avvicina a 2 al crescere di n, mentre con pochi dati la varianza scende più in fretta del previsto. Il bias invece non dipende da n: dimezzare la varianza non basta se l'errore è sistematico.

Regolarizzazione: scegliere dove stare sulla curva

La regolarizzazione aggiunge al problema dei minimi quadrati una penalità sulla grandezza dei pesi: la Ridge del capitolo 02 minimizza ||y − Xw||² + α||w||². Un α piccolo lascia il modello libero (poco bias, molta varianza); un α grande lo irrigidisce (più bias, meno varianza). Il listato usa un polinomio di grado 9 con le misure standardizzate, perché la penalità tratta tutti i pesi allo stesso modo, e misura bias e varianza per sei valori di α.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

rng = np.random.default_rng(2)
f = lambda x: np.sin(2 * np.pi * x)
sigma, n, M = 0.3, 30, 300
xg = np.linspace(0.02, 0.98, 49)[:, None]
print("alpha     bias^2  varianza  somma+sigma^2")
for alpha in (1e-6, 1e-3, 1e-2, 1e-1, 1.0, 10.0):
    P = np.empty((M, len(xg)))
    for m in range(M):
        x = rng.uniform(0, 1, (n, 1))
        y = f(x[:, 0]) + sigma * rng.standard_normal(n)
        mod = make_pipeline(PolynomialFeatures(9, include_bias=False), StandardScaler(), Ridge(alpha=alpha)).fit(x, y)
        P[m] = mod.predict(xg)
    b2 = np.mean((P.mean(axis=0) - f(xg[:, 0])) ** 2)
    v = np.mean(P.var(axis=0))
    print(f"{alpha:7.0e}  {b2:6.4f}  {v:8.4f}  {b2 + v + sigma ** 2:8.4f}")

# Output:
# alpha     bias^2  varianza  somma+sigma^2
#   1e-06  0.0002    0.1588    0.2490
#   1e-03  0.0007    0.1763    0.2671
#   1e-02  0.0031    0.0281    0.1211
#   1e-01  0.0318    0.0320    0.1538
#   1e+00  0.0694    0.0261    0.1855
#   1e+01  0.1306    0.0205    0.2411

La varianza scende di molto, da 0,1763 a 0,0281, passando da alpha=1e-3 a 1e-2, mentre il bias al quadrato cresce da 0,0007 a 0,0031 e poi a 0,1306 per alpha=10. L'errore totale ha un minimo a α = 0,01 (0,1211). Due osservazioni. Anche alpha=1e-6 è già una regolarizzazione: confrontata con il grado 9 senza penalità della simulazione precedente (varianza 8,713), la varianza qui è 0,1588, perché la penalità e la standardizzazione stabilizzano il sistema. Inoltre il minimo ottenuto (0,1211) è vicino ma non migliore del grado 3 non penalizzato (0,1095): regolarizzare un modello troppo grande non sostituisce la scelta di una complessità giusta. Come scegliere α senza conoscere la verità si fa con la validazione incrociata del capitolo 05.

Ridurre la varianza mediando

La seconda strada contro la varianza è mediare più modelli. Il bagging (bootstrap aggregating) addestra molti modelli, ciascuno su un campione ottenuto estraendo con reinserimento dal training, e ne media le previsioni. Un albero senza limiti di profondità è il caso estremo: bias quasi nullo, varianza alta. Il listato confronta un albero singolo e il bagging di 50 alberi sullo stesso problema, con 100 campioni di 60 punti.

import numpy as np
from sklearn.ensemble import BaggingRegressor
from sklearn.tree import DecisionTreeRegressor

rng = np.random.default_rng(3)
f = lambda x: np.sin(2 * np.pi * x)
sigma, n, M = 0.3, 60, 100
xg = np.linspace(0.02, 0.98, 49)[:, None]
modelli = {
    "un albero senza limiti": lambda: DecisionTreeRegressor(random_state=0),
    "bagging di 50 alberi": lambda: BaggingRegressor(DecisionTreeRegressor(), n_estimators=50, random_state=0),
}
for nome, fabbrica in modelli.items():
    P = np.empty((M, len(xg)))
    for m in range(M):
        x = rng.uniform(0, 1, (n, 1))
        y = f(x[:, 0]) + sigma * rng.standard_normal(n)
        P[m] = fabbrica().fit(x, y).predict(xg)
    b2 = np.mean((P.mean(axis=0) - f(xg[:, 0])) ** 2)
    print(f"{nome:24s} bias^2={b2:.4f} varianza={np.mean(P.var(axis=0)):.4f}")

# Output:
# un albero senza limiti   bias^2=0.0012 varianza=0.0936
# bagging di 50 alberi     bias^2=0.0007 varianza=0.0430

La varianza si dimezza (da 0,0936 a 0,0430) e il bias al quadrato non peggiora (da 0,0012 a 0,0007). Non si divide per 50, perché i 50 alberi sono addestrati su campioni rimescolati dello stesso training e quindi i loro errori sono correlati; la riduzione c'è, ma ha un limite. Questa è l'idea che i random forest spingono più in là decorrelando gli alberi, tema del capitolo 06.

Oltre il quadro classico

La rappresentazione a «U» dell'errore di test in funzione della complessità non è l'ultima parola. Belkin e colleghi, in «Reconciling modern machine learning practice and the bias-variance trade-off» (arXiv:1812.11118), discutono il comportamento di modelli con più parametri che esempi, in cui l'errore di test può calare di nuovo dopo il picco di interpolazione (double descent). In questo articolo non si riproduce l'effetto e non se ne riportano i risultati: gli esperimenti qui sopra restano nel regime classico, con modelli piccoli rispetto ai dati.

Che cosa sai fare adesso

Sai scrivere la scomposizione dell'errore quadratico e verificarla con una simulazione, distinguere con le curve di apprendimento un modello rigido da uno instabile, scegliere fra più dati, più regolarizzazione e un modello meno complesso in base a ciò che vedi, e sai perché mediare alberi riduce la varianza e non il bias. Resta da capire come stimare l'errore di test quando i dati sono pochi: 05 - Cross-validation e selezione del modello.

Nella stessa serie trovi anche: 01 - Introduzione al Machine Learning: tassonomia e workflow, 03 - Regressione Logistica e classificazione binaria e 12 - Deployment ML in produzione: serving e monitoring.