LASSO e discesa del gradiente

La ridge regressionUna buona prestazione sul training non basta: serve stimare quella su dati nuovi. La cross-validation (K-fold: $k$ parti, ciascuna a turno come test, errore medio; Monte Carlo: $k$ divisioni casuali con quota di test $q$; leave-one-out se $k=n$) evita di dipendere da una sola divisione casuale. L'errore atteso si scompone in $\text{bias}^2+\text{varianza}+\sigma^2$: i modelli semplici fanno underfitting (bias alto), quelli complessi overfitting (varianza alta). La regolarizzazione aggiunge alla perdita una penalità: la ridge regression minimizza $|y-X\beta|^2+\lambda\sum_{j\ge1}\beta_j^2$ e ha soluzione $\beta=(X^TX+\lambda\tilde I)^{-1}X^Ty$ (l'intercetta non si penalizza, le feature si standardizzano): riduce i coefficienti, rende l'inversa stabile con feature collineari, e $\lambda$ è un iperparametro scelto con la validazione (cross-validation annidata per non contaminare il test).Overfitting, ridge regression e cross-validation → riduce i coefficienti ma non li azzera mai. Cambiando la penalità si ottiene il LASSO, che produce modelli sparsi; ma per questo non c'è una formula chiusa e serve un algoritmo iterativo: la discesa del gradiente, che si usa anche per la regressione logistica e per le reti neurali (Lezione 10 · LASSO e discesa del gradiente, Lezione 12 · Laboratorio di regolarizzazione). Gli esercizi sono Esercizio - Discesa del gradiente per il LASSO e scelta del passo e Esercizio - Cross-validation annidata per il LASSO sul dataset prostate; la variante con funzione di costo asimmetrica è Esercizio - Regressione quantile con funzione pinball e discesa del gradiente.

Dalla ridge al LASSO

Si parte dalla funzione di costo regolarizzata $J=\sum_{i=1}^n[y_i-\hat y_i]^2+\lambda R$, con $R$ penalità sulla complessità. Per $R=\sum_{j=1}^p\beta_j^2$ si ha la ridge ($L_2$). E se si sceglie un'altra penalità?

Definizione (LASSO). Least Absolute Shrinkage and Selection Operator: penalità $L_1$, somma dei valori assoluti dei coefficienti, $$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|,$$ con $\beta_0$ non penalizzato (come nella ridge) e feature standardizzate (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 →). In forma matriciale, con $X$ comprensiva della colonna di uni: $J(\beta)=\lVert y-X\beta\rVert^2+\lambda\sum_{j\ge1}|\beta_j|$.

Esempio. Con $\beta=(\beta_0,\beta_1,\beta_2)=(3;\,2;\,-0{,}5)$ e $\lambda=4$ la penalità vale $4\,(|2|+|-0{,}5|)=4\cdot2{,}5=10$, mentre la ridge con lo stesso $\lambda$ varrebbe $4\,(4+0{,}25)=17$ (il valore assoluto penalizza di più i coefficienti piccoli, il quadrato quelli grandi).

Trace plot e sparsità

Se si riporta il valore di ciascun coefficiente in funzione di $\lambda$ (trace plot), sul dataset prostate (97 osservazioni, 8 feature cliniche, target: livello di log-PSA) si vede che i coefficienti della ridge si riducono gradualmente verso zero ma restano diversi da zero, mentre quelli del LASSO arrivano a zero uno dopo l'altro: a $\lambda$ elevato rimangono pochi coefficienti non nulli. Una soluzione con un sottoinsieme di coefficienti uguali a zero si dice sparsa.

Perché una soluzione sparsa è utile.

• Efficienza: previsioni più veloci, meno memoria, riaddestramento più rapido.
• Robustezza: meno sensibile ai dati rumorosi.
• Interpretabilità: meno variabili, modello più facile da capire.
• Gestione: sistema più semplice da mettere in produzione (anche su dispositivi piccoli).

Il LASSO fa la selezione delle feature, che è un passo di preprocessing, direttamente nella fase di modellazione.

Perché la penalità $L_1$ azzera i coefficienti: il caso con una variabile

Il caso più semplice si può risolvere a mano. Si abbia una sola feature con $x^Tx=1$ (colonna di norma 1 e centrata) e si indichi con $z$ il coefficiente OLS, $z=x^Ty$. Allora $\lVert y-x b\rVert^2=\text{costante}+(b-z)^2$ (perché $\lVert y-xb\rVert^2=y^Ty-2b\,x^Ty+b^2$ e $b^2-2bz=(b-z)^2-z^2$). Il costo da minimizzare è $$f(b)=(b-z)^2+\lambda|b|.$$

• Ridge (penalità $\lambda b^2$): $f'(b)=2(b-z)+2\lambda b=0\Rightarrow b=\dfrac{z}{1+\lambda}$. Il coefficiente è diviso per $1+\lambda$: non è mai zero se $z\neq0$.
• LASSO. Per $b>0$: $f'(b)=2(b-z)+\lambda=0\Rightarrow b=z-\lambda/2$, accettabile solo se $z>\lambda/2$. Per $b<0$: $f'(b)=2(b-z)-\lambda=0\Rightarrow b=z+\lambda/2$, accettabile solo se $z<-\lambda/2$. Se $|z|\le\lambda/2$ nessuno dei due vale e il minimo è in $b=0$ (nel punto angoloso il valore assoluto ha pendenza che va da $-1$ a $+1$ e riesce a bilanciare la pendenza $-2z$ del termine quadratico quando $|2z|\le\lambda$). In una formula (soglia morbida, soft-thresholding): $$\hat b_{\text{LASSO}}=\operatorname{sign}(z)\,\max\!\Big(|z|-\frac\lambda2,\,0\Big).$$

Esempio ($\lambda=1$). Per $z=0{,}8$: LASSO $0{,}8-0{,}5=0{,}3$, ridge $0{,}8/2=0{,}4$. Per $z=0{,}3$: LASSO $0$ (perché $0{,}3\le0{,}5$), ridge $0{,}15$. Per $z=-0{,}9$: LASSO $-0{,}4$, ridge $-0{,}45$. Il LASSO sottrae una quantità fissa e taglia a zero i coefficienti piccoli; la ridge li riduce in proporzione e non li annulla.

Interpretazione geometrica

Il problema regolarizzato è equivalente (per un opportuno $s$ legato a $\lambda$) a minimizzare l'errore quadratico con un vincolo: $\sum_j|\beta_j|\le s$ per il LASSO, $\sum_j\beta_j^2\le s$ per la ridge (relazione con i moltiplicatori di Lagrange: Massimi e minimi vincolati e moltiplicatori di LagrangeGli estremi di f sul vincolo g = c si cercano per sostituzione, per parametrizzazione o con i moltiplicatori di Lagrange: se f, g sono C¹, P0 è un estremo vincolato e ∇g(P0) ≠ 0, esiste λ con ∇f(P0) = λ∇g(P0) (curva di livello di f tangente al vincolo). Si risolve il sistema ∇f = λ∇g, g = c, si aggiungono i punti del vincolo con ∇g = 0 e si confrontano i valori; se il vincolo è compatto, Weierstrass garantisce massimo e minimo.Massimi e minimi vincolati e moltiplicatori di Lagrange →). Con due coefficienti il vincolo è un rombo per il LASSO (con gli spigoli sugli assi) e un cerchio per la ridge. Le curve di livello dell'errore sono ellissi centrate nella soluzione OLS; la soluzione regolarizzata è il primo punto in cui l'ellisse tocca il vincolo. Per il cerchio il contatto avviene in un punto qualsiasi del bordo (di solito con entrambe le coordinate non nulle); il rombo ha spigoli sugli assi, e l'ellisse tende a toccare uno spigolo, dove una coordinata è zero: soluzione sparsa.

(L'ellisse del disegno è stata scelta in modo da toccare lo spigolo; con altri dati il contatto cade altrove: la sparsità è una tendenza della penalità $L_1$, non una garanzia.)

Ridge e LASSO a confronto

Elastic Net. Combina le due penalità per superarne i difetti: $$J(\beta)=\sum_{i=1}^n(y_i-x_i^T\beta)^2+\lambda_1\sum_{j=1}^p|\beta_j|+\lambda_2\sum_{j=1}^p\beta_j^2,\qquad\lambda_1=\alpha\lambda,\ \ \lambda_2=(1-\alpha)\lambda.$$ Per $\alpha=1$ è il LASSO puro, per $\alpha=0$ la ridge pura. Il vincolo geometrico è un rombo con i lati incurvati, tra il rombo e il cerchio.

La discesa del gradiente

Idea

Per il LASSO non si può risolvere $\nabla J=0$ in forma chiusa. Come per ogni addestramento supervisionato (minimizzare una perdita $\mathcal L$, che può essere l'MSE o la cross-entropy della classificazione), si cerca il vettore di pesi che minimizza la perdita totale: $$W^\ast=\arg\min_WJ(W),\qquad J(W)=\frac1n\sum_{i=1}^n\mathcal L\big(f(x^{(i)};W),y^{(i)}\big).$$ Il gradiente $\nabla J(W)$, vettore delle derivate parziali (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 →, Gradiente e direzione di massima crescitaIl gradiente è il vettore delle derivate parziali ∇f(p) = (∂₁f(p), …, ∂ₙf(p)). Se f è C¹ (derivate parziali continue) vale la formula del gradiente D_u f(p) = ∇f(p)·u: tutte le derivate direzionali si ottengono dalle parziali e u ↦ D_u f(p) è lineare. Tra i versori, la crescita è massima lungo ∇f/|∇f| (pendenza |∇f|), minima lungo −∇f/|∇f| (pendenza −|∇f|), nulla lungo le direzioni ortogonali al gradiente. Utili: ∇|x| = x/|x|, ∇φ(|x|) = φ'(|x|) x/|x|.Gradiente e direzione di massima crescita →), punta nella direzione di massima crescita di $J$ ed è perpendicolare alle curve di livello. Per scendere si va nella direzione opposta.

Formula (discesa del gradiente, gradient descent). Partendo da un $W$ iniziale, si ripete $$W\leftarrow W-\eta\,\nabla J(W),$$ dove $\eta>0$ è il passo (learning rate), fino a convergenza (per esempio finché $\lVert W_{\text{nuovo}}-W\rVert<\text{tol}$ o dopo un numero massimo di iterazioni).

inizializza W (a caso, ~N(0, σ²), oppure con le scelte sotto)
ripeti fino a convergenza:
    g = gradiente di J in W
    W = W - eta * g
restituisci W

Esempio in una variabile. $J(w)=w^2$ ha gradiente $2w$ e minimo in $w=0$. L'aggiornamento è $w\leftarrow w-\eta\cdot2w=(1-2\eta)w$. Partendo da $w_0=4$:

Scelta del passo $\eta$

Il passo è critico: troppo grande, l'algoritmo diverge; troppo piccolo, la convergenza è lenta. Indicazioni:

1. Limite superiore. Per l'errore quadratico $J(\beta)=\lVert y-X\beta\rVert^2$ il gradiente è $\nabla J=2X^T(X\beta-y)$. Sia $\beta^\ast$ il minimo e $e=\beta-\beta^\ast$ l'errore. Poiché $X^T(X\beta^\ast-y)=0$, un passo dà $e\leftarrow e-2\eta X^TXe=(I-2\eta X^TX)e$. Si decompone $e$ nella base degli autovettori di $X^TX$ (simmetrica, autovalori $\mu_i\ge0$: 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 →, 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 →): la componente $i$ si moltiplica per $1-2\eta\mu_i$ a ogni passo, quindi tende a zero se e solo se $|1-2\eta\mu_i|<1$, cioè $0<\eta<1/\mu_i$. Per tutte le direzioni: $$\eta<\frac1L,\qquad L=\mu_{\max}(X^TX)$$ (il massimo autovalore di $X^TX$, la costante di Lipschitz delle slide). Con la perdita normalizzata per $n$ o con altri fattori il limite cambia per lo stesso fattore.
2. Passo adattivo. Si fa decrescere $\eta$ durante l'addestramento: passi grandi all'inizio (lontano dal minimo), piccoli alla fine: $$\eta_t=\frac{\eta_0}{1+\gamma t},$$ con $\eta_0$ passo iniziale, $t$ l'iterazione, $\gamma$ parametro di decadimento (per esempio $0{,}01$).

Esempio completo (stessi dati della regressione lineare). Per $x=[1,2,3,4]$, $y=[2,3,5,4]$ la matrice $X^TX=\begin{bmatrix}4&10\\10&30\end{bmatrix}$ ha autovalori $0{,}599$ e $33{,}40$; quindi $\eta<1/33{,}40=0{,}0299$. Con $\eta=0{,}02$ e $\beta_0=(0,0)$: il gradiente è $2X^T(X\beta-y)=-2X^Ty=-2(14,39)=(-28,-78)$ e il primo passo dà $\beta=(0{,}56;\ 1{,}56)$. Al secondo passo il gradiente è $(7{,}68;\ 26{,}8)$ e $\beta=(0{,}406;\ 1{,}024)$; al terzo $(-4{,}27;\ -8{,}43)$ e $\beta=(0{,}492;\ 1{,}193)$. Dopo 2000 iterazioni $\beta=(1{,}5;\ 0{,}8)$, la soluzione OLS (la direzione dell'autovalore piccolo converge lentamente). Con $\eta=0{,}035>0{,}0299$ il metodo diverge: i valori del gradiente passano da $(0{,}98;\ 2{,}73)$ a $(-0{,}23;\ -0{,}96)$, poi $(1{,}49;\ 3{,}94)$, e crescono in modulo.

La stessa convergenza si osserva nel laboratorio (MSE per $\eta=0{,}0001$, $0{,}001$, $0{,}009$: troppo piccolo, giusto, troppo grande). Con $\eta=0{,}009$ la perdita diverge, ma con il passo adattivo $\eta_t=0{,}009/(1+0{,}1t)$ si stabilizza.

Da dove partire

Il punto iniziale può essere casuale, ma alcune regole accelerano la ricerca:

1. se $\lambda$ è grande: partenza da zero (la soluzione ha molti coefficienti nulli);
2. se le feature sono molto correlate: partire dalla soluzione ridge;
3. se le feature sono indipendenti e $\lambda$ è piccolo: partire dalla soluzione OLS;
4. nel dubbio: ridge o OLS, che sono punti di partenza ragionevoli.

Il gradiente del LASSO: il subgradiente

Il valore assoluto non è derivabile in $0$ (Derivata - definizione e significatoLa derivata è il limite del rapporto incrementale; geometricamente è la pendenza della retta tangente. f è derivabile in x0 se e solo se f(x) = f(x0) + f'(x0)(x − x0) + o(x − x0); derivabile implica continua, non viceversa. Derivata destra e sinistra, punti angolosi, flessi a tangente verticale, cuspidi.Derivata - definizione e significato →): la derivata vale $+1$ per $w>0$, $-1$ per $w<0$ e per $w=0$ non è definita (il grafico ha uno spigolo). Si usa il subgradiente: un numero $g$ tale che la retta di pendenza $g$ per il punto sta sotto la funzione, $$f(w')\ge f(w)+g\,(w'-w)\quad\text{per ogni }w'.$$ Per $f(w)=|w|$ l'insieme dei subgradienti è $$\partial|w|=\begin{cases}+1&w>0\\-1&w<0\\ \text{qualsiasi valore in }[-1,1]&w=0.\end{cases}$$ Si può quindi usare $\operatorname{sign}(w)$ con la convenzione $\operatorname{sign}(0)=0$. Il termine dell'errore ha gradiente $2X^T(X\beta-y)$ (come per l'OLS), perciò con $\beta=(\beta_0,\dots,\beta_p)$:

Formula (subgradiente del LASSO). $$\nabla_\beta J=2X^T(X\beta-y)+\lambda\begin{bmatrix}0\\ \operatorname{sign}(\beta_1)\\ \vdots\\ \operatorname{sign}(\beta_p)\end{bmatrix},$$ con $0$ al posto di $\operatorname{sign}(\beta_0)$ perché l'intercetta non si penalizza.

Esempio. Con $\lambda=1$, $\beta=(0{,}5;\ -2;\ 0)$ e componente dell'errore $2X^T(X\beta-y)=(1;\ 3;\ -4)$: il subgradiente è $(1+0;\ 3+1\cdot(-1);\ -4+1\cdot0)=(1;\ 2;\ -4)$.

La funzione di costo è convessa in tutti e tre i casi (OLS, ridge, LASSO): somma di una quadratica convessa e di una penalità convessa (Funzioni convesse in più variabiliUn insieme C è convesso se contiene il segmento tra due suoi punti; f: C → R è convessa se f(tx + (1−t)y) ≤ t f(x) + (1−t) f(y) per t in [0,1], cioè il grafico sta sotto le corde. Se f è differenziabile, è convessa se e solo se f(y) ≥ f(x) + ∇f(x)·(y − x) (il grafico sta sopra ogni piano tangente). Se f è C² su un aperto convesso, è convessa se e solo se l'hessiana è semidefinita positiva in ogni punto; se è definita positiva ovunque f è strettamente convessa (non vale il viceversa: x⁴). Per una funzione convessa ogni punto critico è un minimo globale.Funzioni convesse in più variabili →); nelle curve di livello (cerchi per ridge e OLS standardizzati, forma «a diamante» smussato per il LASSO) c'è un solo minimo. Per questo la discesa del gradiente non resta intrappolata in minimi locali. Con il subgradiente, però, i coefficienti tendono a oscillare intorno a zero senza fermarsi esattamente su zero; per ottenere zeri esatti si usano metodi dedicati, come la soglia morbida vista sopra (coordinata per coordinata) o la Least Angle Regression (LARS), citati nelle slide con i «metodi per sottocoordinate».

import numpy as np

def lasso_gradient(X, y, beta, lam):
    g = 2 * X.T @ (X @ beta - y)          # parte dell'errore quadratico
    pen = lam * np.sign(beta)             # subgradiente di |beta_j|
    pen[0] = 0                            # niente penalità sull'intercetta
    return g + pen

def train_lasso_gd(X, y, lam, lr, max_iter=1000, tol=1e-6):
    X = np.hstack([np.ones((len(X), 1)), X])          # X standardizzata + colonna di uni
    beta = np.random.randn(X.shape[1], 1)
    for _ in range(max_iter):
        new = beta - lr * lasso_gradient(X, y, beta, lam)
        if np.linalg.norm(new - beta) < tol:           # criterio di arresto
            break
        beta = new
    return beta

y deve avere forma $(n,1)$ come beta, altrimenti la sottrazione X @ beta - y fa il broadcasting e i risultati sono sbagliati in silenzio.

Scegliere il passo e $\lambda$

Il passo $\eta$ e $\lambda$ sono iperparametri. Il laboratorio (Lezione 12 · Laboratorio di regolarizzazione) li sceglie con la cross-validation annidata sul dataset prostate (Overfitting, ridge regression e cross-validationUna buona prestazione sul training non basta: serve stimare quella su dati nuovi. La cross-validation (K-fold: $k$ parti, ciascuna a turno come test, errore medio; Monte Carlo: $k$ divisioni casuali con quota di test $q$; leave-one-out se $k=n$) evita di dipendere da una sola divisione casuale. L'errore atteso si scompone in $\text{bias}^2+\text{varianza}+\sigma^2$: i modelli semplici fanno underfitting (bias alto), quelli complessi overfitting (varianza alta). La regolarizzazione aggiunge alla perdita una penalità: la ridge regression minimizza $|y-X\beta|^2+\lambda\sum_{j\ge1}\beta_j^2$ e ha soluzione $\beta=(X^TX+\lambda\tilde I)^{-1}X^Ty$ (l'intercetta non si penalizza, le feature si standardizzano): riduce i coefficienti, rende l'inversa stabile con feature collineari, e $\lambda$ è un iperparametro scelto con la validazione (cross-validation annidata per non contaminare il test).Overfitting, ridge regression e cross-validation →): $K=8$ fold esterni, $D=4$ fold interni, griglia di 28 combinazioni (4 passi $\times$ 7 valori di $\lambda$, da $10^{-4}$ a $100$), con le feature standardizzate dentro ogni fold usando solo il training. Nelle combinazioni con passo troppo grande e $\lambda$ elevato l'algoritmo diverge (valori inf/nan): la selezione le scarta perché l'errore di validazione risulta il peggiore.

Errori tipici

• Penalizzare l'intercetta o non standardizzare le feature prima del LASSO.
• Scegliere $\eta$ troppo grande: l'errore cresce invece di scendere (controllare la curva di perdita).
• Usare lo stesso insieme per scegliere $\lambda$ e per stimare la prestazione.
• Aspettarsi dal subgradiente zeri esatti: servono soglia morbida o metodi dedicati.
• Dimenticare il verso: il gradiente punta verso l'alto, l'aggiornamento sottrae $\eta\nabla J$.
• Con feature molto correlate il LASSO ne sceglie una in modo instabile: meglio Elastic Net o ridge.