Overfitting, ridge regression e cross-validation

Nella Regressione lineareNell'apprendimento supervisionato si impara una funzione $F(x)$ dagli esempi $(x,y)$: regressione se $y$ è continua, classificazione se è categorica. Il modello lineare è $F_\beta(x)=\beta_0+\beta_1x_1+\dots+\beta_px_p=X\beta$ (con una colonna di uni per $\beta_0$) e i parametri si scelgono minimizzando l'errore quadratico medio $\mathrm{MSE}=\frac1n\sum_i(y_i-F_\beta(x_i))^2$, funzione convessa dei parametri. Annullando il gradiente di $J(\beta)=|y-X\beta|^2$ si ottengono le equazioni normali $X^TX\beta=X^Ty$ e la soluzione dei minimi quadrati ordinari $\beta=(X^TX)^{-1}X^Ty$. Il coefficiente di determinazione $R^2=1-SS_{res}/SS_{tot}$ misura la qualità del fit (0 = come la media, negativo = peggio della media). Un modello va valutato su un test set mai usato per addestrare: l'errore sul training è ottimistico e un polinomio di grado alto lo azzera senza generalizzare.Regressione lineare → si è visto che l'errore sul training peggiora la stima della qualità e che un polinomio di grado alto passa per tutti i punti senza generalizzare. Questa nota (Lezione 8 · Ridge regression, overfitting e cross-validation) spiega come misurare l'overfitting (cross-validation, bias e varianza) e come limitarlo (regolarizzazione, in particolare la ridge). La variante con penalità a valore assoluto è in LASSO e discesa del gradienteIl LASSO è la regressione regolarizzata con penalità $L_1$: minimizza $\sum_i(y_i-x_i^T\beta)^2+\lambda\sum_{j\ge1}|\beta_j|$. A differenza della ridge, porta alcuni coefficienti esattamente a zero (soluzione sparsa, selezione delle feature): geometricamente le curve di livello dell'errore toccano il vincolo $\sum|\beta_j|\le s$ (un rombo) in uno spigolo. Non ha formula chiusa, quindi si minimizza con la discesa del gradiente $W\leftarrow W-\eta,\nabla J(W)$, usando il subgradiente $\operatorname{sign}(\beta_j)$ per il valore assoluto (nel punto 0 qualunque valore in $[-1,1]$). Il passo $\eta$ è critico: per l'errore quadratico converge se $\eta<1/\mu_{\max}(X^TX)$; si può usare $\eta_t=\eta_0/(1+\gamma t)$. L'Elastic Net combina le penalità $L_1$ e $L_2$ con $\lambda_1=\alpha\lambda$, $\lambda_2=(1-\alpha)\lambda$.LASSO e discesa del gradiente →; il laboratorio è Lezione 9 · Laboratorio di regressione lineare e ridge con Esercizio - OLS, R quadro e aumento polinomiale delle feature, Esercizio - Ridge regression con standardizzazione e coefficienti in unità originali e Esercizio - Cross-validation k-fold e Monte Carlo per OLS e ridge sul dataset Advertising.

Perché non basta il training

Ogni modello visto ha parametri da addestrare (per esempio i $\beta$ della regressione lineare) e per addestrarli servono due ingredienti: dati e una procedura di addestramento con una funzione di costo. Per l'OLS l'addestramento è semplice (formula chiusa). Ma l'addestramento non basta: bisogna testare il modello per

• capire le sue prestazioni reali: generalizza? Quali prestazioni aspettarsi nel mondo reale?
• scegliere, tra le alternative di una pipeline, quella migliore.

Si divide il dataset in una parte di training e una di test (per esempio $80\%$ e $20\%$, a caso). Ma la scelta casuale è «sicura»? Con dataset piccoli il risultato può cambiare molto a seconda di quali punti finiscono nel test.

Cross-validation

Nella modellazione i dati si dividono in training e validation (costruzione del modello) e test (stima finale della prestazione). Per non dipendere da una singola divisione si ripete la valutazione più volte e si media: è la cross-validation (CV).

K-fold

Definizione (K-fold cross-validation). Si mescolano i dati e si dividono in $k$ parti (fold) di uguale dimensione. Per $i=1,\dots,k$: si usa il fold $i$ come insieme di valutazione e i restanti $k-1$ per addestrare il modello; si calcola la metrica (per esempio l'MSE) $\mathrm{MSE}_i$ sul fold $i$. La stima finale è la media $\overline{\mathrm{MSE}}=\frac1k\sum_{i=1}^k\mathrm{MSE}_i$.

Ogni dato è usato esattamente una volta per la valutazione e $k-1$ volte per l'addestramento. L'unica scelta di progetto è il numero di fold $k$. Nel caso estremo $k=n$ si parla di leave-one-out (si lascia fuori un solo dato per volta).

Esempio. Con $n=20$ dati e $k=5$ si formano 5 fold da 4 dati (dati 1-4, 5-8, 9-12, 13-16, 17-20 dopo il mescolamento): a ogni giro si addestra su 16 dati e si valuta su 4. Per il polinomio quadratico dell'esempio del laboratorio gli MSE dei 5 fold sono $0{,}0106;\ 0{,}0219;\ 0{,}0204;\ 0{,}0194;\ 0{,}0243$ e la media vale $\frac{0{,}0106+0{,}0219+0{,}0204+0{,}0194+0{,}0243}{5}=\frac{0{,}0966}{5}=0{,}0193$.

Monte Carlo cross-validation (MCCV)

Definizione (MCCV). Si ripete $k$ volte: divisione casuale del dataset in training e test, con una frazione fissata $q$ di dati nel test; addestramento; calcolo di $\mathrm{MSE}_i$. Si media: $\frac1k\sum_i\mathrm{MSE}_i$. Le scelte di progetto sono due: il numero di ripetizioni $k$ e la quota di test $q$.

Il nome viene dal metodo Monte Carlo, che stima una quantità ripetendo un campionamento casuale e facendo la media (Legge dei grandi numeri e metodo Monte CarloSe X₁, X₂, ... sono i.i.d. con media μ, la media campionaria X̄ₙ = (X₁ + ... + Xₙ)/n converge a μ: in probabilità (legge debole, dimostrata con Chebyshev se la varianza è finita: P(|X̄ₙ − μ| > ε) ≤ σ²/(nε²)) e quasi certamente (legge forte). Metodo Monte Carlo: ∫ g = E[g(U)] si stima con la media di g(U₁), ..., g(Uₙ) per uniformi indipendenti.Legge dei grandi numeri e metodo Monte Carlo →): nel laboratorio lo si mostra stimando $\pi$ come $4\times$ (frazione di punti casuali di un quadrato che cadono nel cerchio inscritto), con errore che scende come $1/\sqrt{\text{numero di punti}}$. Qui il «campione» è la divisione casuale dei dati, e mediare più divisioni riduce la varianza della stima.

Confronto

Guardare solo la media a volte basta, ma è meglio considerare la distribuzione degli errori sui fold (ad esempio con un box plot, Correlazione e visualizzazione dei datiLa correlazione di Pearson $r=\sum(X_i-\bar X)(Y_i-\bar Y)/\big(\sqrt{\sum(X_i-\bar X)^2}\sqrt{\sum(Y_i-\bar Y)^2}\big)\in[-1,1]$ misura la relazione lineare tra due variabili (covarianza divisa per le deviazioni standard); correlazione non implica causalità. Serve a capire quali variabili contano per il target e a eliminare quelle quasi duplicate (|r| molto alto). Gli indicatori di sintesi non bastano (quartetto di Anscombe, Datasaurus): vanno affiancati ai grafici: istogramma, KDE, box plot, violin plot, heatmap di correlazione, scatter plot e matrice di scatter plot.Correlazione e visualizzazione dei dati →).

Esempio (dalle slide, California Housing). K-fold con $k=10$ e MCCV con $k=100$, $q=0{,}2$, confrontando un modello a una feature (MedInc) con uno a quattro feature: con più feature la media dell'MSE è più bassa (meno bias) ma il box plot è più largo (più varianza).

Bias e varianza

Gli errori di un modello sono di tipi diversi.

Definizione (bias). L'incapacità di un metodo di cogliere la vera relazione tra ingresso e uscita. Un modello troppo semplice (la retta quando i dati seguono una parabola) ha bias alto e non la catturerà mai, per quanti dati si abbiano.

Definizione (varianza). La sensibilità del modello ai dati di addestramento: se, cambiando il campione di training, le previsioni cambiano molto, la varianza è alta. Un modello molto flessibile (polinomio con 20 coefficienti) ha bias basso sul training ma varianza alta.

Con la cross-validation i due difetti si distinguono: il polinomio di grado alto è ottimo sul training (basso bias) ma pessimo sui fold di valutazione, e le prestazioni cambiano molto da fold a fold (varianza alta); la retta ha errori alti ma simili in tutti i fold (varianza bassa, bias alto).

Teorema (decomposizione bias-varianza). Sia $y=f(x)+\varepsilon$ con rumore $E[\varepsilon]=0$, $\operatorname{Var}(\varepsilon)=\sigma^2$, indipendente dal training. Per un punto fissato $x$ e un modello $\hat f$ addestrato su un campione casuale: $$E\big[(y-\hat f(x))^2\big]=\underbrace{\big(f(x)-E[\hat f(x)]\big)^2}_{\text{bias}^2}+\underbrace{E\big[(\hat f(x)-E[\hat f(x)])^2\big]}_{\text{varianza}}+\sigma^2.$$

Dimostrazione, passo per passo. (I valori attesi sono quelli di Valore attesoIl valore atteso E[X] = Σ x p_X(x) è la media dei valori di X pesata con le loro probabilità (esiste se la serie converge assolutamente); per una funzione g vale E[g(X)] = Σ g(x) p_X(x) senza trovare la legge di g(X), ed E è lineare: E[aX + bY + c] = aE[X] + bE[Y] + c.Valore atteso →, la varianza è in Varianza e momentiI momenti E[X^k] e i momenti centrati E[(X − μ)^k] descrivono la forma di una legge; la varianza Var(X) = E[(X − μ)²] = E[X²] − E[X]² misura quanto X si disperde attorno alla media, vale Var(aX + b) = a² Var(X) e Var(X) = 0 solo se X è costante.Varianza e momenti →.)

1. Si sostituisce $y=f+\varepsilon$ e si sviluppa il quadrato: $E[(f+\varepsilon-\hat f)^2]=E[(f-\hat f)^2]+2E[\varepsilon(f-\hat f)]+E[\varepsilon^2]$.
2. Il termine misto è nullo: $\varepsilon$ è indipendente da $\hat f$ e $E[\varepsilon]=0$, quindi $E[\varepsilon(f-\hat f)]=E[\varepsilon]\,E[f-\hat f]=0$. Inoltre $E[\varepsilon^2]=\sigma^2$.
3. Resta $E[(f-\hat f)^2]$. Si somma e si sottrae $E[\hat f]$: $f-\hat f=(f-E\hat f)+(E\hat f-\hat f)$. Elevando al quadrato: $(f-E\hat f)^2+2(f-E\hat f)(E\hat f-\hat f)+(E\hat f-\hat f)^2$.
4. Si prende il valore atteso. Il primo termine è una costante. Nel secondo, $f-E\hat f$ è una costante e $E[E\hat f-\hat f]=0$, quindi si annulla. Il terzo è per definizione la varianza di $\hat f$. ∎

Esempio. Se per un punto $f=3$, la media delle previsioni su molti campioni è $E\hat f=2{,}5$ con varianza $0{,}04$ e $\sigma^2=0{,}01$, l'errore quadratico atteso è $(3-2{,}5)^2+0{,}04+0{,}01=0{,}25+0{,}04+0{,}01=0{,}30$.

Il rumore $\sigma^2$ è irriducibile; bias e varianza si bilanciano nella complessità del modello:

• Underfitting (sotto-adattamento): modello troppo semplice, i dati non sono sfruttati («ha imparato troppo poco»), associato a bias alto;
• Overfitting (sovra-adattamento): modello troppo complesso che cattura il rumore invece degli andamenti veri, associato a varianza alta: non generalizza.

Il punto fondamentale è che nessun modello è «ottimo» in sé: bisogna trovare l'equilibrio, e la cross-validation permette di farlo guardando dati non usati per addestrare.

Esempio (laboratorio, $n=20$ punti da $y=3x^2-2x+1+\text{rumore}$). Errore medio con 5-fold cross-validation (fold da 4 dati consecutivi) per polinomi OLS: grado 1 $\to0{,}121$ (underfitting: tutti i fold alti), grado 2 $\to0{,}0193$, grado 5 $\to0{,}0168$, grado 9 $\to$ circa $1{,}7\cdot10^{5}$ (un fold esplode: overfitting estremo). Con la ridge sul grado 9 e $\lambda=0{,}01$ la media scende a $0{,}0150$.

Regolarizzazione

Definizione (regolarizzazione). Tecnica per prevenire l'overfitting aggiungendo alla funzione di costo un termine di penalità sulla complessità del modello: si minimizza $$J=\sum_{i=1}^n\big[y_i-\hat y_i\big]^2+\gamma\,R,$$ dove $R$ misura la complessità e $\gamma\ge0$ è il parametro di regolarizzazione, un iperparametro (un valore che si sceglie prima dell'addestramento e che non è appreso dai dati). Se $\gamma=0$ non c'è regolarizzazione; per $\gamma$ grande conta quasi solo la penalità.

L'idea è quella del rasoio di Occam («tra le spiegazioni dei fenomeni si preferisce la più semplice possibile»): a parità di adattamento ai dati, meglio i parametri piccoli. Dato un test (o una validazione), si sceglie il valore dell'iperparametro che dà il miglior compromesso tra complessità e accuratezza.

Ridge regression

Definizione (ridge regression, penalità $L_2$). Con $R=\sum_{j=1}^p\beta_j^2$ (somma dei quadrati dei coefficienti, senza $\beta_0$) si minimizza $$J(\beta)=\sum_{i=1}^n\big(y_i-\beta_0-\sum_{j=1}^p\beta_jx_{ij}\big)^2+\lambda\sum_{j=1}^p\beta_j^2 .$$

Perché non si penalizza l'intercetta. Si vogliono ridurre i pesi delle feature rispetto alla risposta per diminuire la complessità; l'intercetta non è il peso di una feature: è il valore medio della risposta quando tutte le feature sono zero. Penalizzare $\beta_0$ spingerebbe le previsioni verso zero, cosa sbagliata in generale.

Standardizzazione necessaria. La penalità tratta tutti i $\beta_j$ allo stesso modo: se le feature hanno scale diverse, una feature con valori enormi ha coefficienti minuscoli e viene penalizzata poco, e una con valori piccoli ha coefficienti grandi e viene penalizzata molto, per un'unica ragione di unità di misura. Perciò le feature si standardizzano (media $0$, deviazione standard $1$, calcolate sul training, Statistica per il machine learningI dati di un problema ML si organizzano nella matrice di progetto $X$ ($n$ osservazioni, $p$ variabili). La statistica serve a capirli, ripulirli e prepararli: i momenti (media $\mu$, varianza $\sigma^2$, asimmetria, curtosi), i quartili con lo scarto interquartile $\mathrm{IQR}=Q_3-Q_1$ (all'esame senza interpolazione), la moda per i dati categorici. Con queste quantità si imputano i dati mancanti (media o mediana), si eliminano le variabili costanti e si standardizza con lo z-score $z=(x-\mu)/\sigma$, usando sempre media e deviazione standard del solo training set.Statistica per il machine learning →).

Derivazione della soluzione in forma chiusa. Si prendono le feature standardizzate (colonne a media zero) e si mette in $X$ la colonna di uni come primo elemento. Si indica con $\tilde I$ la matrice identità con zero nella posizione $(0,0)$ (così la penalità non tocca $\beta_0$): $$J(\beta)=(y-X\beta)^T(y-X\beta)+\lambda\,\beta^T\tilde I\beta.$$

1. Si sviluppa il primo termine come per l'OLS: $y^Ty-2y^TX\beta+\beta^TX^TX\beta$.
2. Gradiente del termine di penalità: $\beta^T\tilde I\beta$ è una forma quadratica con matrice simmetrica $\tilde I$, e la regola $\frac{\partial}{\partial\beta}\beta^TA\beta=2A\beta$ dà $2\lambda\tilde I\beta$ (Differenziabilità e gradientef è differenziabile in x0 se f(x) = f(x0) + ∇f(x0)·(x − x0) + o(‖x − x0‖): vicino a x0 il grafico si confonde con il piano tangente z = f(x0) + ∇f(x0)·(x − x0). Differenziabile ⇒ continua, derivabile e D_v f = ∇f·v; derivate parziali continue ⇒ differenziabile. Il gradiente indica la direzione di massima crescita (pendenza ‖∇f‖) ed è ortogonale alle curve di livello.Differenziabilità e gradiente →).
3. Gradiente totale: $-2X^Ty+2X^TX\beta+2\lambda\tilde I\beta$.
4. Lo si pone uguale a zero e si dividono i due membri per $2$: $X^TX\beta+\lambda\tilde I\beta=X^Ty$.
5. Si raccoglie $\beta$: $(X^TX+\lambda\tilde I)\beta=X^Ty$.
6. Se la matrice tra parentesi è invertibile (lo è sempre per $\lambda>0$, vedi sotto): $$\boxed{\hat\beta_{\text{ridge}}=(X^TX+\lambda\tilde I)^{-1}X^Ty}.$$

(Con $y$ centrata e le feature standardizzate si può togliere l'intercetta, perché $\hat\beta_0=\bar y$; la formula diventa $\hat\beta=(X^TX+\lambda I)^{-1}X^Ty$ con $I$ identità completa. Nelle slide si indica semplicemente $(X^TX+\lambda I)^{-1}X^Ty$.)

Formula (ridge regression). $\hat\beta_{\text{ridge}}=(X^TX+\lambda\tilde I)^{-1}X^Ty$. Per $\lambda=0$ coincide con l'OLS; al crescere di $\lambda$ i coefficienti si riducono (shrinkage) verso zero.

Esempio (una sola feature). Dati $x=[-1{,}5;-0{,}5;0{,}5;1{,}5]$ (centrata) e $y_c=[-1{,}5;-0{,}5;1{,}5;0{,}5]$ ($y=[2,3,5,4]$ meno la media $3{,}5$): $X^TX=2{,}25+0{,}25+0{,}25+2{,}25=5$ e $X^Ty_c=2{,}25+0{,}25+0{,}75+0{,}75=4$. Con $\lambda=0$: $\hat\beta=4/5=0{,}8$ (l'OLS). Con $\lambda=1$: $4/(5+1)=0{,}667$; con $\lambda=5$: $4/10=0{,}4$; con $\lambda=20$: $4/25=0{,}16$; per $\lambda\to\infty$: $\hat\beta\to0$. Il trace plot (coefficienti in funzione di $\lambda$) mostra questa discesa.

Collinearità e stabilità. Quando le feature sono molto correlate, $X^TX$ è mal condizionata: ha autovalori vicini a zero (Autovalori e autovettoriUn autovettore è un vettore non nullo che una funzione lineare manda in un suo multiplo; si trovano gli autovalori come radici del polinomio caratteristico det(A − λI) e gli autovettori come nucleo di A − λI. Matrici simili hanno gli stessi autovalori.Autovalori e autovettori →). Per una matrice simmetrica si può scrivere $X^TX=Q\Lambda Q^T$ con $Q$ ortogonale e $\Lambda$ diagonale degli autovalori $\mu_i\ge0$ (Teorema spettrale e forme quadraticheUna funzione lineare è simmetrica se f(v)·w = v·f(w); in una base ortonormale ha matrice simmetrica. Teorema spettrale: f è simmetrica se e solo se esiste una base ortonormale di autovettori, cioè A simmetrica ⇔ PᵀAP diagonale con P ortogonale. Applicato alle forme quadratiche, permette di scriverle come somma di quadrati con gli autovalori come coefficienti.Teorema spettrale e forme quadratiche →, DiagonalizzazioneUna matrice è diagonalizzabile se è simile a una diagonale, cioè se esiste una base di autovettori: allora A = S D S⁻¹ con gli autovettori nelle colonne di S e gli autovalori in D. Criterio: tutti gli autovalori nel campo e molteplicità geometrica uguale a quella algebrica. Le matrici simmetriche reali hanno autovalori reali.Diagonalizzazione →); l'inversa ha autovalori $1/\mu_i$, enormi se $\mu_i\approx0$: i coefficienti OLS variano moltissimo con piccoli cambiamenti dei dati (Matrice inversaL'inversa di una matrice quadrata A è la matrice A⁻¹ con A A⁻¹ = A⁻¹ A = I; esiste se e solo se rango(A) = n e si calcola con Gauss-Jordan riducendo (A | I) fino a (I | A⁻¹).Matrice inversa →). Aggiungendo $\lambda I$ (caso senza intercetta) si ha $X^TX+\lambda I=Q(\Lambda+\lambda I)Q^T$, con autovalori $\mu_i+\lambda\ge\lambda>0$: la matrice è sempre invertibile e l'inversa ha autovalori al più $1/\lambda$. In più, nella base degli autovettori il coefficiente ridge è il coefficiente OLS moltiplicato per $\dfrac{\mu_i}{\mu_i+\lambda}\in(0,1)$: le direzioni con piccola varianza ($\mu_i$ piccolo) sono ridotte quasi a zero, quelle con grande varianza quasi non cambiano.

Esempio (feature quasi identiche). $x_1=[-3,-1,1,3]$, $x_2=x_1+[0{,}1,-0{,}1,-0{,}1,0{,}1]$, $y=[-2{,}9;\,-1{,}3;\,1{,}1;\,3{,}1]$. Allora $X^TX=\begin{bmatrix}20&20\\20&20{,}04\end{bmatrix}$ con determinante $20\cdot20{,}04-20\cdot20=0{,}8$ e autovalori $0{,}02$ e $40{,}02$ (quasi singolare). L'OLS dà $\beta=(0{,}02;\ 1{,}00)$. Se si cambia un solo valore di $y$ di $+0{,}2$ (da $-1{,}3$ a $-1{,}1$), l'OLS diventa $(0{,}51;\ 0{,}50)$: i coefficienti si spostano di circa $0{,}5$ per una variazione minima. La ridge con $\lambda=1$ dà $(0{,}488;\ 0{,}508)$ prima e $(0{,}4925;\ 0{,}4928)$ dopo: quasi identici. Nelle slide si ricorda lo stesso fenomeno con $y=100x_1-99x_2+\dots$: coefficienti enormi e opposti che si compensano.

Un caso con $p=9$ variabili (laboratorio). Con il grado $9$ le feature $x,x^2,\dots,x^9$ sono quasi collineari. I coefficienti OLS (standardizzati) hanno valori come $78{,}5;\ -796{,}4;\ 3818{,}6;\ -10\,063;\ 15\,578;\ -14\,175;\ 7032;\ -1469$ mentre quelli ridge con $\lambda=0{,}1$ stanno tra $-0{,}28$ e $0{,}28$. Training: OLS ha MSE $0{,}0061$ ($R^2=0{,}964$), ridge $0{,}0120$ ($R^2=0{,}928$): sembra peggiore. Test su 20 punti nuovi: OLS $0{,}0251$ ($R^2=0{,}674$), ridge $0{,}0124$ ($R^2=0{,}838$): la ridge generalizza meglio. È l'errore classico confrontare i modelli sul training.

La relazione tra $\lambda$ e le prestazioni è a «U» sul test: poca regolarizzazione dà overfitting, troppa dà underfitting.

Coefficienti in unità originali

Se si standardizza, i coefficienti $\beta_j$ si riferiscono alle feature standardizzate. Per tornare alle unità originali, da $\hat Y=\beta_0+\sum_j\beta_j\dfrac{X_j-\mu_j}{\sigma_j}$ si separa la parte costante e quella che moltiplica $X_j$: $$\hat Y=\Big(\beta_0-\sum_j\beta_j\frac{\mu_j}{\sigma_j}\Big)+\sum_j\frac{\beta_j}{\sigma_j}X_j,$$ quindi $\beta_0^{\text{orig}}=\beta_0-\sum_j\beta_j\mu_j/\sigma_j$ e $\beta_j^{\text{orig}}=\beta_j/\sigma_j$. Esempio. Se $\beta_1=2$ per una feature con $\mu_1=10$, $\sigma_1=4$, e $\beta_0=7$ (una sola feature): $\beta_1^{\text{orig}}=2/4=0{,}5$ e $\beta_0^{\text{orig}}=7-2\cdot10/4=2$. Verifica: per $X=14$ il modello standardizzato dà $7+2\cdot(14-10)/4=9$ e quello originale $2+0{,}5\cdot14=9$ ✓.

Scegliere $\lambda$ senza barare: cross-validation annidata

Se si sceglie $\lambda$ con la cross-validation (grid search, random search), si seleziona il modello migliore sulla validazione e poi lo si valuta su un test separato; ma se il test ha già influenzato la scelta (anche solo indirettamente) la stima finale è ottimistica. La soluzione è la cross-validation annidata (nested CV), con due cicli:

1. Ciclo esterno: i dati sono divisi in training (con validation) e test; il test serve solo alla valutazione finale.
2. Ciclo interno: all'interno del training, per ciascuna combinazione di iperparametri si fa una cross-validation: si addestra sul training interno e si valuta sul validation interno, e si mediano le prestazioni.
3. Si sceglie la combinazione migliore, si riaddestra sul training esterno completo con quella combinazione e si valuta sul test esterno.
4. Si ripete per ogni fold esterno e si aggregano (media e deviazione standard). La combinazione finale è la media (o la moda) di quelle scelte; con essa si addestra un modello definitivo su tutti i dati.

Esempio. Con 5 fold esterni, 4 fold interni e 28 combinazioni di iperparametri (LASSO e discesa del gradienteIl LASSO è la regressione regolarizzata con penalità $L_1$: minimizza $\sum_i(y_i-x_i^T\beta)^2+\lambda\sum_{j\ge1}|\beta_j|$. A differenza della ridge, porta alcuni coefficienti esattamente a zero (soluzione sparsa, selezione delle feature): geometricamente le curve di livello dell'errore toccano il vincolo $\sum|\beta_j|\le s$ (un rombo) in uno spigolo. Non ha formula chiusa, quindi si minimizza con la discesa del gradiente $W\leftarrow W-\eta,\nabla J(W)$, usando il subgradiente $\operatorname{sign}(\beta_j)$ per il valore assoluto (nel punto 0 qualunque valore in $[-1,1]$). Il passo $\eta$ è critico: per l'errore quadratico converge se $\eta<1/\mu_{\max}(X^TX)$; si può usare $\eta_t=\eta_0/(1+\gamma t)$. L'Elastic Net combina le penalità $L_1$ e $L_2$ con $\lambda_1=\alpha\lambda$, $\lambda_2=(1-\alpha)\lambda$.LASSO e discesa del gradiente →) servono $5\cdot28\cdot4=560$ addestramenti per la ricerca, più $5$ per la valutazione.

Codice

import numpy as np
def fit_ridge(X, y, lam):                        # X già standardizzata
    X = np.hstack([np.ones((len(X), 1)), X])     # colonna di uni per beta_0
    I = np.eye(X.shape[1]); I[0, 0] = 0          # non si penalizza l'intercetta
    return np.linalg.solve(X.T @ X + lam * I, X.T @ y)

idx = np.random.permutation(len(X)); k = 5       # K-fold
for fold in np.array_split(idx, k):
    train = np.setdiff1d(idx, fold)
    mu, sd = X[train].mean(0), X[train].std(0)   # statistiche del solo training
    beta = fit_ridge((X[train]-mu)/sd, y[train], 10.0)
    pred = beta[0] + ((X[fold]-mu)/sd) @ beta[1:]
    mse = np.mean((y[fold] - pred) ** 2)

np.linalg.solve risolve il sistema $(X^TX+\lambda\tilde I)\beta=X^Ty$ senza calcolare l'inversa.

Errori tipici

• Scegliere il modello (o $\lambda$) guardando l'errore sul training: la ridge sembra peggiore dell'OLS sul training ma è migliore sul test.
• Penalizzare l'intercetta, o non standardizzare le feature prima della ridge.
• Calcolare media e deviazione standard su tutti i dati e non sul solo training di ogni fold (fuga di informazione).
• Usare lo stesso insieme per scegliere $\lambda$ e per stimare la prestazione finale (serve la cross-validation annidata).
• Pensare che $\lambda$ si «impari» dai dati dell'addestramento: è un iperparametro.
• Confondere bias e varianza: il bias non diminuisce con più dati se il modello è troppo semplice.