37 - Boosting: AdaBoost, Gradient Boosting e cenni su XGBoost

Prepari un esame con un quaderno di esercizi. Al primo giro svolgi tutto; poi non rifai gli esercizi che sai già, ma ti concentri su quelli che hai sbagliato. Al giro successivo, di nuovo: dedichi attenzione agli errori rimasti. Dopo dieci giri hai un metodo di studio che, esercizio dopo esercizio, ha corretto le proprie lacune. Il boosting fa lo stesso: costruisce una sequenza di modelli semplici, in cui ciascuno si concentra sugli errori dei precedenti.

La differenza con la Random Forest è di principio. La foresta addestra molti alberi profondi in parallelo e indipendenti e li media per ridurre la varianza. Il boosting addestra alberi piccoli (spesso di profondità 1-4) uno dopo l'altro e li somma per ridurre il bias. In questa lezione derivi e implementi a mano i due algoritmi fondamentali: AdaBoost e il Gradient Boosting.

Parole da chiarire subito

• Weak learner (apprendista debole): un modello appena meglio del caso. Lo stump (letteralmente «ceppo») è l'albero decisionale più piccolo possibile: una sola domanda, due foglie (vedi alberi decisionali).
• Peso di un esempio: un numero che dice quanto quell'esempio conta nella funzione di errore. Peso alto = «sbagliarlo costa di più».
• Residuo: la differenza y − previsione, cioè quanto il modello sbaglia ancora su un esempio.
• Gradiente: per una funzione di molte variabili, il vettore delle derivate parziali; punta nella direzione in cui la funzione cresce più rapidamente. Andare nella direzione opposta la fa calare più rapidamente: è l'idea della discesa del gradiente, che trovi in La discesa del gradiente.

Le quattro parole formano una catena: il weak learner è il mattone, il peso dice a ciascun mattone quali esempi guardare, il residuo e il gradiente dicono cosa resta da correggere. AdaBoost usa i pesi; il Gradient Boosting usa residui e gradiente. Se ricordi questa divisione, il resto della lezione è la stessa idea scritta in due modi.

Il disegno: a ogni round il modello (la linea tratteggiata) sbaglia alcuni esempi; nel round successivo quegli esempi sono disegnati più grandi (peso maggiore), e il nuovo stump è spinto a classificarli bene. Il risultato finale è una somma pesata di tutti gli stump.

AdaBoost: ripesare gli esempi

AdaBoost (Adaptive Boosting) lavora su classificazione binaria con etichette y ∈ {−1, +1}. Ogni esempio ha un peso w_i, all'inizio uguale per tutti (1/n). A ogni round t:

1. Addestra un weak learner h_t (tipicamente uno stump: un albero con un solo test) tenendo conto dei pesi.
2. Calcola l'errore pesato: ε_t = Σ_i w_i · [h_t(x_i) ≠ y_i].
3. Calcola il peso del learner: α_t = ½ · ln( (1 − ε_t) / ε_t ).
4. Aggiorna i pesi: w_i ← w_i · exp( −α_t · y_i · h_t(x_i) ), poi normalizza perché sommino a 1.

I quattro passi sono un ciclo, non un elenco: l'errore del passo 2 decide il peso del passo 3, e il peso decide come cambiano i dati che vedrà il learner del round successivo. Nessun learner conosce gli altri; l'unico canale di comunicazione fra i round sono i pesi degli esempi.

La previsione finale è un voto pesato: H(x) = segno( Σ_t α_t · h_t(x) ).

Perché proprio queste formule

Il fattore y_i·h_t(x_i) vale +1 se il learner ha ragione e −1 se sbaglia. Quindi l'aggiornamento moltiplica il peso degli esempi corretti per e^(−α) (minore di 1, scendono) e quello degli errori per e^(+α) (maggiore di 1, salgono). Con α_t = ½ ln((1−ε)/ε) si ha e^α = √((1−ε)/ε). Verifichiamo cosa succede alla massa di peso totale, prima di normalizzare:

massa degli esempi corretti dopo l'aggiornamento = (1 − ε) · e^(−α) = (1 − ε) · √(ε/(1−ε)) = √(ε(1−ε))
massa degli esempi sbagliati dopo l'aggiornamento = ε · e^(+α)      = ε · √((1−ε)/ε)     = √(ε(1−ε))

Le due masse sono uguali: dopo l'aggiornamento il peso totale degli errori precedenti è esattamente il 50%. Il learner successivo si trova di fronte a un problema in cui il learner precedente vale quanto una moneta, quindi è costretto a trovare qualcosa di nuovo. Il peso α_t è grande quando ε è piccolo (learner affidabile), nullo quando ε = 0,5 (come tirare a caso), negativo se ε > 0,5 (il suo voto viene invertito).

Implementazione a mano

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import AdaBoostClassifier

class AdaBoostDaZero:
    def __init__(self, n_round=50):
        self.n_round = n_round

    def fit(self, X, y):                       # y in {-1, +1}
        n = len(y)
        w = np.full(n, 1.0 / n)
        self.stump, self.alfa = [], []
        self.storia_pesi = []
        for t in range(self.n_round):
            h = DecisionTreeClassifier(max_depth=1, random_state=0).fit(X, y, sample_weight=w)
            pred = h.predict(X)
            eps = w[pred != y].sum()
            eps = np.clip(eps, 1e-10, 1 - 1e-10)
            alfa = 0.5 * np.log((1 - eps) / eps)
            w = w * np.exp(-alfa * y * pred)
            w = w / w.sum()
            self.stump.append(h); self.alfa.append(alfa); self.storia_pesi.append(w.copy())
        return self

    def decision_function(self, X):
        return sum(a * h.predict(X) for a, h in zip(self.alfa, self.stump))

    def predict(self, X):
        return np.where(self.decision_function(X) >= 0, 1, -1)

X, y01 = make_classification(n_samples=1500, n_features=10, n_informative=5, n_redundant=0,
                             class_sep=1.0, random_state=0)
y = 2 * y01 - 1
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)

stump = DecisionTreeClassifier(max_depth=1).fit(Xtr, ytr)
ada = AdaBoostDaZero(n_round=100).fit(Xtr, ytr)
lib = AdaBoostClassifier(n_estimators=100, random_state=0).fit(Xtr, ytr)
print("singolo stump:      test %.3f" % stump.score(Xte, yte))
print("AdaBoost a mano:    test %.3f" % (ada.predict(Xte) == yte).mean())
print("AdaBoost sklearn:   test %.3f" % lib.score(Xte, yte))
print("alfa dei primi 5 round:", np.round(ada.alfa[:5], 3))

singolo stump:      test 0.800
AdaBoost a mano:    test 0.849
AdaBoost sklearn:   test 0.849
alfa dei primi 5 round: [0.708 0.411 0.258 0.284 0.143]

Nel listato fit inizializza i pesi uguali (1/n) e per ogni round: addestra uno stump passando sample_weight=w (la libreria dà più importanza agli esempi con peso alto nel calcolo dell'impurità); calcola l'errore pesato eps sommando i pesi degli errori; limita eps con np.clip per evitare log(0) o divisioni per zero se lo stump fosse perfetto; ricava alfa; aggiorna e normalizza i pesi; salva stump, alfa e pesi. decision_function è la somma pesata dei voti; predict ne prende il segno. Notare y = 2*y01 - 1: AdaBoost richiede etichette −1/+1, perché la formula dell'aggiornamento usa il prodotto y·h(x).

Uno stump da solo fa 0.800 sul test: gli stump sono «deboli» per costruzione. Cento stump pesati insieme arrivano a 0.849, identico alla libreria. Nota gli alfa dei primi round (0.708, 0.411, 0.258, ...): i primi learner ricevono il voto più forte e il peso scende man mano che i problemi si fanno più difficili. I pesi della libreria sono il doppio dei nostri (misurato: estimator_weights_ dà 1.416, 0.822, 0.516, ... per i primi round, cioè 2 × 0.708, 2 × 0.411, 2 × 0.258): scikit-learn usa la variante SAMME, che nel caso binario scala il voto di un fattore costante. Il fattore non cambia il segno del voto pesato né i pesi normalizzati degli esempi, ed è per questo che le due accuratezze coincidono.

Il lato oscuro: AdaBoost e le etichette sbagliate

L'aggiornamento dei pesi ha un effetto collaterale importante: un esempio che viene sbagliato round dopo round vede il suo peso crescere esponenzialmente. Se quell'esempio è semplicemente etichettato male, l'algoritmo spende sempre più risorse per «imparare» un errore. Misuriamolo: invertiamo il 15% delle etichette di training e osserviamo dove va a finire il peso.

rng = np.random.default_rng(0)
sporche = rng.random(len(ytr)) < 0.15                 # 15% delle etichette di training invertite
ytr_sporco = np.where(sporche, -ytr, ytr)

ada_s = AdaBoostDaZero(n_round=100).fit(Xtr, ytr_sporco)
w_finale = ada_s.storia_pesi[-1]
print("quota di esempi con etichetta sbagliata: %.1f%%" % (100 * sporche.mean()))
print("quota del peso totale che portano dopo 100 round: %.1f%%" % (100 * w_finale[sporche].sum()))
print("test, addestrato su etichette pulite: %.3f   su etichette sporche: %.3f" % (
    (ada.predict(Xte) == yte).mean(), (ada_s.predict(Xte) == yte).mean()))

quota di esempi con etichetta sbagliata: 13.7%
quota del peso totale che portano dopo 100 round: 27.7%
test, addestrato su etichette pulite: 0.849   su etichette sporche: 0.833

Il listato: rng.random(len(ytr)) < 0.15 estrae per ogni esempio un numero casuale e lo marca «sporco» se minore di 0,15: quindi circa il 15% (qui 13,7%, essendo una scelta casuale). Si invertono le etichette sporche con np.where. Si riaddestra AdaBoost e si legge il vettore dei pesi dell'ultimo round: la somma dei pesi sugli esempi sporchi dice quanta attenzione l'algoritmo dedica agli errori di etichetta.

Il 13,7% degli esempi (quelli effettivamente invertiti) arriva a portare il 27,7% del peso totale: il doppio della loro quota. L'accuratezza sul test scende da 0.849 a 0.833. L'effetto qui è moderato; cresce con la quota di rumore e con il numero di round. AdaBoost è sensibile al rumore sulle etichette: se i tuoi dati sono sporchi, preferisci un boosting con perdite più robuste o con regolarizzazione forte, oppure pulisci prima i dati.

Gradient Boosting: la discesa del gradiente nello spazio delle funzioni

Il Gradient Boosting generalizza l'idea a qualsiasi funzione di perdita derivabile. Il modello è una somma di funzioni F_M(x) = F₀ + η · Σ_m h_m(x). A ogni passo si aggiunge un nuovo albero h_m che ci sposta nella direzione che fa calare più rapidamente la perdita: è una discesa del gradiente in cui i «parametri» sono i valori predetti per ciascun esempio.

Con la perdita quadratica L = ½ (y − F)² il gradiente rispetto alla previsione F è

∂L/∂F = −(y − F)     →     −∂L/∂F = y − F    (il residuo!)

Quindi, per la perdita quadratica, «seguire il gradiente negativo» significa semplicemente addestrare il nuovo albero sui residui del modello attuale. Per una perdita qualunque, si addestra l'albero sul gradiente negativo (detto pseudo-residuo). L'algoritmo:

1. Parti da una costante: F₀ = media di y (per la perdita quadratica è la costante ottima).
2. Per m = 1...M: calcola i residui r_i = y_i − F_{m−1}(x_i); addestra un albero piccolo h_m su (x_i, r_i); aggiorna F_m = F_{m−1} + η · h_m.

La costante iniziale è il modello più semplice possibile; ogni albero successivo sposta le previsioni un po' verso i residui, scalato da η. Il modello finale è la costante più la somma degli alberi, quindi per prevedere devi attraversarli tutti: M piccolo vuol dire modello snello, M grande modello accurato ma più lento e a rischio di overfitting, come mostra l'esperimento che segue.

Il learning rate η (tra 0 e 1, spesso 0,05-0,1) è lo shrinkage: ogni albero contribuisce solo una frazione di ciò che suggerisce. Passi piccoli e molti alberi generalizzano meglio di pochi passi grandi, al prezzo di più alberi da addestrare. Implementiamolo per la regressione e confrontiamo con scikit-learn:

from sklearn.datasets import make_friedman1
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error

Xf, yf = make_friedman1(n_samples=1000, n_features=10, noise=1.0, random_state=0)
Xa, Xb, ya, yb = train_test_split(Xf, yf, test_size=0.3, random_state=0)

class GBRegDaZero:
    def __init__(self, n_alberi=200, eta=0.1, prof=2):
        self.n_alberi, self.eta, self.prof = n_alberi, eta, prof

    def fit(self, X, y):
        self.f0 = y.mean()
        F = np.full(len(y), self.f0)
        self.alberi = []
        for m in range(self.n_alberi):
            residuo = y - F                                     # gradiente negativo della perdita quadratica
            h = DecisionTreeRegressor(max_depth=self.prof, random_state=0).fit(X, residuo)
            F = F + self.eta * h.predict(X)
            self.alberi.append(h)
        return self

    def predict(self, X, fino_a=None):
        F = np.full(len(X), self.f0)
        for h in self.alberi[:fino_a]:
            F = F + self.eta * h.predict(X)
        return F

gb = GBRegDaZero(n_alberi=200, eta=0.1, prof=2).fit(Xa, ya)
lib = GradientBoostingRegressor(n_estimators=200, learning_rate=0.1, max_depth=2, random_state=0).fit(Xa, ya)

print("MSE test a mano:  %.4f" % mean_squared_error(yb, gb.predict(Xb)))
print("MSE test sklearn: %.4f" % mean_squared_error(yb, lib.predict(Xb)))
d = np.abs(gb.predict(Xb) - lib.predict(Xb))
print("differenza tra le previsioni: media %.3f, massima %.3f (scala di y: deviazione std %.2f)" % (d.mean(), d.max(), yb.std()))

MSE test a mano:  2.3872
MSE test sklearn: 2.3882
differenza tra le previsioni: media 0.005, massima 0.223 (scala di y: deviazione std 5.16)

Il listato: make_friedman1 genera un problema di regressione classico con 10 feature (solo 5 realmente utili) e una relazione non lineare. GBRegDaZero.fit parte da f0, la media di y (la miglior costante), poi per ogni albero: calcola il residuo del modello corrente, ci addestra un DecisionTreeRegressor di profondità 2, e aggiorna F aggiungendo eta volte la previsione dell'albero. predict rifà la somma su nuovi dati; il parametro fino_a permette di usare solo i primi M alberi, comodo per tracciare l'errore durante la costruzione.

Il nostro ciclo e il GradientBoostingRegressor raggiungono lo stesso errore (2.387 contro 2.388). Le previsioni non sono identiche al decimale (differenza media 0.005 contro una deviazione standard di y pari a 5.16, massima 0.223; una causa plausibile sono gli spareggi tra divisioni equivalenti, non isolata), ma sono la stessa idea: per la perdita quadratica la libreria non fa altro che questo ciclo sui residui. Ora notiamo qualcosa che nella Random Forest non succedeva.

Il boosting può fare overfitting, e va fermato

Aggiungere alberi a una foresta è sempre innocuo. In un modello boosted, ogni albero aggiunto riduce l'errore di training, ma dopo un certo punto comincia ad adattarsi al rumore. Guardiamo l'errore di test al crescere degli alberi, con learning rate alto:

gb_alto = GBRegDaZero(n_alberi=600, eta=0.5, prof=4).fit(Xa, ya)
for M in (1, 5, 20, 50, 100, 300, 600):
    print("alberi %3d   MSE train %.3f   MSE test %.3f" % (
        M, mean_squared_error(ya, gb_alto.predict(Xa, M)), mean_squared_error(yb, gb_alto.predict(Xb, M))))

alberi   1   MSE train 11.938   MSE test 13.652
alberi   5   MSE train 2.336   MSE test 4.544
alberi  20   MSE train 0.587   MSE test 3.803
alberi  50   MSE train 0.091   MSE test 3.645
alberi 100   MSE train 0.006   MSE test 3.747
alberi 300   MSE train 0.000   MSE test 3.768
alberi 600   MSE train 0.000   MSE test 3.768

Il listato addestra un modello con learning rate alto (0,5) e alberi più profondi (4) per rendere visibile l'overfitting, e valuta l'errore a un numero crescente di alberi sfruttando predict(X, fino_a=M).

L'errore di training scende fino a zero, mentre quello di test tocca il minimo (3.645) intorno a 50 alberi e poi risale (3.747 a 100, 3.768 oltre): è l'overfitting nella sua forma più chiara, lo stesso che hai visto con gli alberi e con le SVM. La risalita qui è piccola perché la funzione è liscia, ma il segnale è chiaro: dopo il minimo ogni albero in più peggiora la generalizzazione. Si controlla con tre leve: learning rate più basso (con più alberi), alberi poco profondi, e soprattutto l'early stopping, cioè fermarsi quando l'errore su un set di validazione smette di migliorare (in scikit-learn: n_iter_no_change e validation_fraction). Il numero di alberi è un iperparametro da scegliere in validazione, non da impostare «abbastanza alto».

Da Gradient Boosting a XGBoost: cenni

XGBoost (Chen e Guestrin, 2016) è un'implementazione del gradient boosting con un obiettivo regolarizzato e molta ingegneria; lo sviluppo al secondo ordine che usa è, come dice lo stesso articolo, già presente nella letteratura (Friedman e coll.). Qui ne vedi l'idea; non lo implementiamo.

• Sviluppo al secondo ordine: usa non solo il gradiente g_i della perdita ma anche la derivata seconda h_i. Per un albero con foglie j, il peso ottimo di una foglia è w*_j = − G_j / (H_j + λ), con G_j e H_j somme di g_i e h_i degli esempi che finiscono nella foglia, e λ un termine di regolarizzazione L2.
• Regolarizzazione nell'obiettivo (la formula del guadagno è quella dell'articolo originale): penalizza il numero di foglie (γ per foglia) e la grandezza dei pesi (λ). Il guadagno di una divisione è ½ [ G_L²/(H_L+λ) + G_R²/(H_R+λ) − (G_L+G_R)²/(H_L+H_R+λ) ] − γ; se è negativo, la divisione non si fa.
• Ingegneria: ricerca delle divisioni su istogrammi, gestione nativa dei valori mancanti, parallelismo, sottocampionamento di righe e colonne.

Le prime due voci cambiano la matematica dell'albero: la foglia non ha più la media dei residui ma un valore che tiene conto anche della curvatura, e le divisioni inutili vengono scartate dalla penalità γ. La terza voce non cambia l'algoritmo, ma è ciò che lo rende veloce su milioni di righe. Per un confronto di principio basta il Gradient Boosting che abbiamo costruito.

LightGBM e CatBoost seguono la stessa idea con scelte diverse (la documentazione di LightGBM spiega che cresce l'albero foglia per foglia, scegliendo la foglia con la maggiore riduzione di perdita). Se non vuoi installare altro, scikit-learn include HistGradientBoostingClassifier, della stessa famiglia a istogrammi:

from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier

Xc, yc = make_classification(n_samples=3000, n_features=20, n_informative=6, flip_y=0.03, random_state=0)
Xc_tr, Xc_te, yc_tr, yc_te = train_test_split(Xc, yc, test_size=0.3, stratify=yc, random_state=0)

hgb = HistGradientBoostingClassifier(early_stopping=True, validation_fraction=0.2,
                                     n_iter_no_change=15, random_state=0).fit(Xc_tr, yc_tr)
rf = RandomForestClassifier(n_estimators=300, random_state=0).fit(Xc_tr, yc_tr)
print("alberi usati dal boosting con early stopping: %d (massimo consentito 100)" % hgb.n_iter_)
print("test boosting: %.3f    test random forest: %.3f" % (hgb.score(Xc_te, yc_te), rf.score(Xc_te, yc_te)))

alberi usati dal boosting con early stopping: 72 (massimo consentito 100)
test boosting: 0.917    test random forest: 0.908

Il listato: HistGradientBoostingClassifier è un boosting a istogrammi (le feature vengono discretizzate, per default in al più 255 intervalli, per accelerare la ricerca delle soglie). Con early_stopping=True mette da parte il 20% del training (validation_fraction) e si ferma quando il punteggio di validazione non migliora per 15 iterazioni consecutive (n_iter_no_change). n_iter_ riporta quanti alberi sono stati realmente usati.

Il boosting si è fermato da solo a 72 alberi e fa 0.917 contro 0.908 della foresta da 300 alberi: una differenza di 0.009 su 900 campioni, cioè circa 8 esempi. Non c'è un vincitore universale: sui dati tabellari le due famiglie sono spesso vicine e la differenza dipende dal dataset e dalla cura messa nella ricerca degli iperparametri.

Trade-off: boosting contro Random Forest

• Random Forest: più semplice da configurare, robusta ai valori di default, parallelizzabile per albero, tollera meglio il rumore sulle etichette, stima OOB gratuita.
• Boosting: spesso più accurato dopo la messa a punto; addestramento sequenziale (non parallelizzabile per albero); più iperparametri (learning rate, profondità, numero di alberi, regolarizzazione); rischio di overfitting se non si ferma.

Le due famiglie scambiano semplicità con accuratezza potenziale: la foresta funziona bene quasi senza messa a punto, il boosting chiede di scegliere con cura learning rate, profondità e numero di alberi, e in cambio può arrivare più in alto. Il confronto sopra (0,917 contro 0,908 su un solo dataset) non decide nulla: serve provare entrambi con la validazione.

Trappole reali

• Numero di alberi scelto sul test set: usa early stopping con un set di validazione dedicato.
• Etichette sporche con AdaBoost: i pesi esplodono sugli esempi sbagliati (vedi sopra).
• Learning rate e numero di alberi vanno tarati insieme: dimezzare η richiede più alberi per ottenere lo stesso risultato.
• Ordine temporale ignorato: nelle serie temporali la validazione per early stopping deve venire dopo il training nel tempo, altrimenti è leakage.
• Confrontare modelli con configurazioni non ottimizzate: una Random Forest tarata contro un boosting con valori di default non è un confronto equo (e viceversa).

Le cinque trappole hanno una radice comune: il boosting si adatta fino a dove glielo lasci fare. Il numero di alberi, il rumore sulle etichette, l'ordine del tempo e la taratura sono tutte forme dello stesso errore, cioè lasciare che il modello vinca sul training senza un controllo indipendente.

Learning rate e numero di alberi: due manopole che si compensano

Nel Gradient Boosting l'effetto complessivo di ciò che hai appreso è η · M, il learning rate per il numero di alberi: un η piccolo richiede più alberi per raggiungere lo stesso adattamento. Perché allora non usare sempre η molto piccolo? Perché costa più alberi (più tempo) e, soprattutto, perché la letteratura (J. Friedman, 2001) e la pratica mostrano che i passi piccoli generalizzano meglio: ogni albero apporta una correzione modesta, quindi il modello non si sbilancia su singoli alberi rumorosi. Lo verifichiamo con staged_predict, un metodo di scikit-learn che restituisce le previsioni dopo 1, 2, ..., M alberi, per ottenere l'intera curva dell'errore di test in un solo addestramento. Tre learning rate, fino a 800 alberi di profondità 2.

print("learning rate   alberi al minimo   MSE minimo   MSE con 50 alberi   con 200   con 800")
for eta in (0.5, 0.1, 0.02):
    g = GradientBoostingRegressor(n_estimators=800, learning_rate=eta, max_depth=2, random_state=0).fit(Xa, ya)
    mse = np.array([mean_squared_error(yb, p) for p in g.staged_predict(Xb)])
    print("%-14g  %14d          %10.3f        %8.3f          %6.3f       %6.3f" % (eta, mse.argmin() + 1, mse.min(), mse[49], mse[199], mse[799]))

learning rate   alberi al minimo   MSE minimo   MSE con 50 alberi   con 200   con 800
0.5                         81               2.346           2.598           2.401        2.510
0.1                        742               2.001           4.202           2.388        2.004
0.02                       800               2.445          14.327           5.288        2.445

Per ciascun learning rate il listato addestra 800 alberi una sola volta e, con staged_predict, calcola l'errore quadratico sul test dopo ogni albero; poi riporta il numero di alberi in cui l'errore è minimo e il minimo, più tre punti della curva. Con η = 0,5 il minimo (2,346) si raggiunge già con 81 alberi e dopo risale: il modello impara in fretta e inizia a fare overfitting. Con η = 0,1 servono 742 alberi e il minimo è migliore, 2,001. Con η = 0,02 la curva sta ancora scendendo a 800 alberi (2,445): servirebbero molti più alberi, e il minimo non è stato raggiunto.

Due letture. La prima: la combinazione 0,1 con quasi 750 alberi batte 0,5 con 81 di quasi il 15% di errore in meno: il passo piccolo con tanti alberi è meglio. La seconda, pratica: il numero di alberi giusto dipende dal learning rate, per cui non ha senso fissare M a priori; si fissa η (tipicamente 0,05-0,1), si usa un numero alto di alberi e si lascia all'early stopping il compito di fermarsi.

La profondità degli alberi e il sottocampionamento

Gli alberi del boosting sono «deboli» di proposito: troppo profondi, ciascuno si adatta troppo e il meccanismo correttivo perde senso. La profondità controlla anche l'ordine delle interazioni tra feature che il modello può rappresentare: alberi di profondità 1 (stump) sono modelli additivi, in cui ogni feature contribuisce da sola; con profondità 3 si catturano interazioni tra tre feature. Il sottocampionamento (subsample) addestra ogni albero su una frazione casuale degli esempi, introducendo il rumore benefico del bagging (Friedman 2002, stochastic gradient boosting). Misuriamo l'effetto di entrambi, con η = 0,1.

print("profondità   alberi al minimo   MSE minimo   MSE dopo 300 alberi")
for depth in (1, 2, 3, 5):
    g = GradientBoostingRegressor(n_estimators=300, learning_rate=0.1, max_depth=depth, random_state=0).fit(Xa, ya)
    mse = np.array([mean_squared_error(yb, p) for p in g.staged_predict(Xb)])
    print("%-10d  %15d   %11.3f   %17.3f" % (depth, mse.argmin() + 1, mse.min(), mse[-1]))

print("\nsubsample   MSE medio su 5 semi")
for subsample in (1.0, 0.5, 0.25):
    errors = [mean_squared_error(yb, GradientBoostingRegressor(n_estimators=300, learning_rate=0.1, max_depth=3, subsample=subsample, random_state=s).fit(Xa, ya).predict(Xb)) for s in range(5)]
    print("%-9g   %.3f +- %.3f" % (subsample, np.mean(errors), np.std(errors)))

profondità   alberi al minimo   MSE minimo   MSE dopo 300 alberi
1                       300         2.892               2.892
2                       297         2.241               2.243
3                       266         1.957               1.974
5                       161         2.525               2.534

subsample   MSE medio su 5 semi
1           1.974 +- 0.005
0.5         1.933 +- 0.087
0.25        1.926 +- 0.102

Il primo ciclo varia la profondità con 300 alberi e η = 0,1. Gli stump (profondità 1) arrivano a un errore di 2,89 e non hanno ancora smesso di migliorare con 300 alberi; la profondità 2 raggiunge 2,24, la 3 il minimo, 1,96, con 266 alberi; la 5 è peggiore, 2,53, e il suo minimo arriva prima (161 alberi) perché ogni albero impara molto. La funzione dei dati (Friedman 1) contiene termini di interazione tra coppie di feature e un prodotto, quindi gli stump additivi non bastano: è la ragione per cui un po' di profondità aiuta. Troppa, però, ritorna l'overfitting. La regola pratica: profondità 3-6 per i dati tabellari, da scegliere con la validazione.

Il secondo ciclo varia il sottocampionamento con profondità 3. Con 1,0 (tutti gli esempi) il risultato è quasi deterministico (1,974 ± 0,005: l'unica casualità è nei pareggi); con 0,5 e 0,25 l'errore medio è un po' più basso (1,933 e 1,926) ma la deviazione tra semi sale (0,087 e 0,102): la differenza non è significativa con cinque semi. Il sottocampionamento è, qui, un'opzione che non danneggia e che velocizza (si addestra su metà dei dati), non un miglioramento dimostrato. Non sempre un trucco della letteratura porta un beneficio misurabile su un dato problema.

Gradient Boosting per la classificazione: costruirlo a mano

Per la regressione con perdita quadratica il pseudo-residuo era il residuo y − F. Per la classificazione binaria, il modello somma i punteggi F(x) in scala logit e la probabilità è p = σ(F). La perdita è la log-loss di la regressione logistica, e il suo gradiente rispetto a F è, come hai derivato lì, p − y: quindi il pseudo-residuo (gradiente negativo) è r = y − p. Il passo di costruzione è lo stesso: si addestra un albero di regressione sui r e si aggiorna F. C'è una finezza (Friedman 2001): invece di usare come valore di ogni foglia la media dei residui, si fa un passo di Newton sulla perdita e il valore della foglia diventa γ = Σ r / Σ p(1 − p), il rapporto tra il gradiente e la curvatura della log-loss (p(1 − p) è la derivata seconda). Implementiamo tutto a mano e confrontiamo con GradientBoostingClassifier.

from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import log_loss, accuracy_score

X_cl, y_cl = make_classification(n_samples=1500, n_features=10, n_informative=4, flip_y=0.05, class_sep=0.9, random_state=0)
X_cl_tr, X_cl_te, y_cl_tr, y_cl_te = train_test_split(X_cl, y_cl, test_size=0.3, stratify=y_cl, random_state=0)
sigmoid = lambda z: 1 / (1 + np.exp(-z))

class ManualGBClassifier:
    def __init__(self, n_trees=100, eta=0.1, depth=2):
        self.n_trees, self.eta, self.depth = n_trees, eta, depth

    def fit(self, X, y):
        p0 = y.mean()
        self.f0 = np.log(p0 / (1 - p0))                                    # partenza: il logit della frequenza di base
        F = np.full(len(y), self.f0)
        self.trees = []
        for _ in range(self.n_trees):
            p = sigmoid(F)
            residual = y - p                                                # pseudo-residuo = gradiente negativo della log-loss
            tree = DecisionTreeRegressor(max_depth=self.depth, random_state=0).fit(X, residual)
            leaf = tree.apply(X)                                            # in quale foglia cade ogni esempio
            gamma = {l: residual[leaf == l].sum() / (p[leaf == l] * (1 - p[leaf == l])).sum() for l in np.unique(leaf)}   # passo di Newton per foglia
            F = F + self.eta * np.array([gamma[l] for l in leaf])
            self.trees.append((tree, gamma))
        return self

    def decision_function(self, X):
        F = np.full(len(X), self.f0)
        for tree, gamma in self.trees:
            F = F + self.eta * np.array([gamma.get(l, 0.0) for l in tree.apply(X)])
        return F

manual = ManualGBClassifier().fit(X_cl_tr, y_cl_tr)
p_manual = sigmoid(manual.decision_function(X_cl_te))
library = GradientBoostingClassifier(n_estimators=100, learning_rate=0.1, max_depth=2, random_state=0).fit(X_cl_tr, y_cl_tr)
p_library = library.predict_proba(X_cl_te)[:, 1]

print("                a mano    scikit-learn")
print("log-loss test    %.4f    %.4f" % (log_loss(y_cl_te, p_manual), log_loss(y_cl_te, p_library)))
print("accuratezza test %.4f    %.4f" % (accuracy_score(y_cl_te, p_manual > 0.5), accuracy_score(y_cl_te, p_library > 0.5)))
print("correlazione tra le probabilità previste: %.4f;  differenza massima: %.3f" % (np.corrcoef(p_manual, p_library)[0, 1], np.abs(p_manual - p_library).max()))

                a mano    scikit-learn
log-loss test    0.3471    0.3476
accuratezza test 0.8689    0.8689
correlazione tra le probabilità previste: 0.9997;  differenza massima: 0.150

La classe parte dal logit della frequenza di base (la miglior costante), poi per ognuno dei 100 round calcola le probabilità correnti, il pseudo-residuo y − p, addestra un albero di regressione di profondità 2 sui residui e, per ogni foglia, calcola il valore di Newton γ = Σ r / Σ p(1 − p); tree.apply(X) restituisce l'indice della foglia di ogni esempio, ed è ciò che permette di sostituire le medie di sklearn con i valori di Newton. La previsione somma le correzioni scalate di η e passa per la sigmoide.

Il confronto: log-loss sul test 0,3471 a mano e 0,3476 con scikit-learn, accuratezza identica (0,8689), correlazione delle probabilità 0,9997. Le due implementazioni non coincidono al decimale (differenza massima 0,15 su un singolo esempio, dovuta a dettagli di spareggio tra divisioni e di inizializzazione), ma sono lo stesso algoritmo: la differenza media è trascurabile e l'ordinamento degli esempi è praticamente identico. È la dimostrazione che il boosting per la classificazione non è magia: è la discesa del gradiente nello spazio delle funzioni, con pseudo-residuo y − p, la stessa quantità «uscita meno bersaglio» della regressione logistica.

Tre prove da fare

Esercizio 1. Ripeti la tabella di learning rate con η = 0,02 e 3 000 alberi: il minimo ora si raggiunge e vale meno di 2,001? Esercizio 2. Nel boosting a mano per la classificazione, sostituisci il passo di Newton con la media dei residui nelle foglie: quanto peggiora la log-loss a 100 round? Esercizio 3. Aggiungi alla regressione di Friedman cinque colonne di rumore e misura se gli stump o gli alberi di profondità 3 ne risentono di più.

Per allenarti sul controllo delle letture dei contatori c'è uno scheletro di esercizi in Python con tre parti da completare: un giro di AdaBoost con la verifica che l'errore pesato torna a 0,5, il Gradient Boosting a costanti confrontato con la formula chiusa, e i moltiplicatori di peso che spiegano perché una lettura trascritta male rovina l'ensemble. Finché i punti da completare sono vuoti, stampa quali esercizi mancano.

Che cosa sai fare adesso

Sai eseguire a mano un giro di AdaBoost e spiegare perché l'errore pesato torna a 0,5, sai costruire un Gradient Boosting per regressione e per classificazione sui pseudo-residui, e sai riconoscere quando fermarlo con la validazione. Sai scegliere learning rate, profondità e sottocampionamento guardando la curva dell'errore, e sai leggere XGBoost e LightGBM come varianti dello stesso schema. Dalla lezione successiva, clustering, togliamo le etichette.

Fonti