41 - La discesa del gradiente: derivata, learning rate e tre varianti

Sei su una montagna con la nebbia fitta e devi scendere a valle. Non vedi il fondovalle, ma senti sotto i piedi da che parte il terreno scende più ripido. Così fai un passo in quella direzione, poi ti fermi, senti di nuovo la pendenza, fai un altro passo. Passo dopo passo arrivi in basso, senza aver mai visto la mappa. La discesa del gradiente (gradient descent) è esattamente questo, applicato alla funzione di costo che hai incontrato nella lezione su ipotesi, errore quadratico medio e funzione di costo: la «montagna» è J(θ), la tua «posizione» sono i parametri θ, e «sentire la pendenza» significa calcolare una derivata.

Nella lezione precedente il minimo di J si trovava con una formula chiusa. Quella strada si chiude quando il modello ha troppi parametri o non è lineare, ed è il caso di ogni rete neurale. La discesa del gradiente è l'algoritmo che resta: funziona per qualsiasi funzione derivabile, ed è il motore dell'addestramento di tutto ciò che vedrai da qui in avanti. Ogni listato è stato eseguito due volte con seme fisso e l'uscita è identica; il blocco PyTorch è stato eseguito con PyTorch 2.12.0 su CPU (senza GPU, con 2 thread). La lezione costruisce l'algoritmo da zero: cos'è una derivata, perché «contro il gradiente» è la direzione giusta, come si sceglie il passo, e le tre varianti (batch, stocastica, mini-batch) con cui si usa nella pratica.

La derivata: la pendenza in un punto

Prendi una funzione f(θ) con un solo parametro. Se ti sposti di un piccolo tratto h, f cambia di f(θ + h) − f(θ). Il rapporto tra i due è la pendenza media in quel tratto; se fai tendere h a zero, ottieni la derivata:

f'(θ) = lim_{h→0} [ f(θ + h) − f(θ) ] / h

La derivata dice «di quanto cambia f per un piccolo spostamento di θ, per unità di spostamento». Se f'(θ) > 0, f cresce spostandosi a destra; se f'(θ) < 0, f scende spostandosi a destra; se f'(θ) = 0, sei in un punto piatto (un minimo, un massimo o una sella). Un esempio che userai spesso: per f(θ) = θ², la derivata è f'(θ) = 2θ. In θ = 3 la pendenza vale 6; in θ = −3 vale −6; in θ = 0 è zero, ed è il minimo.

Non devi fidarti della formula: la puoi controllare numericamente con una differenza finita centrale, che approssima la derivata scegliendo un h minuscolo e guardando da entrambi i lati: [f(θ + h) − f(θ − h)] / 2h.

import numpy as np

f = lambda t: t ** 2
h = 1e-5
t = 3.0
print("derivata numerica in 3: %.6f   (formula 2*theta = %.1f)" % ((f(t + h) - f(t - h)) / (2 * h), 2 * t))

derivata numerica in 3: 6.000000   (formula 2*theta = 6.0)

Il codice definisce f come una funzione anonima (lambda), sceglie un h piccolo ma non troppo (con h troppo vicino alla precisione della macchina, la differenza di due numeri quasi uguali perde cifre significative) e calcola il rapporto simmetrico. Il risultato coincide con 2θ = 6 fino alla sesta cifra. Questa tecnica tornerà utile più avanti come strumento di verifica dei gradienti calcolati a mano.

Con più parametri: derivate parziali e gradiente

La funzione di costo della regressione lineare dipende da due parametri, J(θ₀, θ₁). La derivata parziale ∂J/∂θ₁ è la derivata rispetto a θ₁ tenendo θ₀ fermo, come se l'altra variabile fosse una costante. Il gradiente è il vettore che le raccoglie tutte:

∇J(θ) = ( ∂J/∂θ₀ , ∂J/∂θ₁ )

Il gradiente ha una proprietà geometrica fondamentale: punta nella direzione in cui J cresce più rapidamente (la «salita più ripida»), e la sua lunghezza è la ripidezza. Quindi andare nella direzione opposta, −∇J, è la discesa più ripida. Da qui la regola di aggiornamento:

θ ← θ − α · ∇J(θ)

dove α (alfa) è il learning rate (tasso di apprendimento): la lunghezza del passo, in unità del gradiente. Si ripete finché J non smette di calare in modo significativo. Il passo è automaticamente piccolo vicino al minimo, dove il gradiente tende a zero, e grande lontano da esso.

Il gradiente per la regressione lineare, derivato con la regola della catena

Con h(x) = θ₀ + θ₁x e J = (1/2m) Σ (h(x_i) − y_i)², la regola della catena dice che la derivata di una funzione composta f(g(θ)) è f'(g(θ)) · g'(θ): si deriva «da fuori verso dentro» e si moltiplica. Chiamiamo e_i = h(x_i) − y_i l'errore. Il quadrato è la funzione esterna, e_i quella interna:

∂J/∂θ₀ = (1/2m) Σ 2·e_i · ∂e_i/∂θ₀ = (1/m) Σ e_i · 1
∂J/∂θ₁ = (1/2m) Σ 2·e_i · ∂e_i/∂θ₁ = (1/m) Σ e_i · x_i

Perché ∂e_i/∂θ₀ = 1 e ∂e_i/∂θ₁ = x_i? Perché e_i = θ₀ + θ₁x_i − y_i è lineare nei parametri: θ₀ compare con coefficiente 1 e θ₁ con coefficiente x_i. In forma matriciale, con la matrice X di forma (m, 2) che ha una colonna di uni e una con gli x, il gradiente è un'unica espressione ∇J = (1/m) Xᵀ(Xθ − y), la stessa che avevi azzerato per ottenere l'equazione normale. Ora, invece di risolverla, la usi per fare un passo.

Il learning rate: il parametro che decide tutto

Il disegno mostra tre discese sulla stessa conca. Con un α troppo piccolo si arriva al fondo ma ci vogliono moltissimi passi. Con un α adatto si converge in pochi passi. Con un α troppo grande ogni passo supera il fondo e atterra sul lato opposto, più in alto di prima: la discesa diverge. Possiamo capirlo con precisione sul caso più semplice, J(θ) = θ², con gradiente 2θ. L'aggiornamento diventa:

θ_nuovo = θ − α · 2θ = (1 − 2α) · θ        →        θ_k = (1 − 2α)^k · θ₀

A ogni passo θ viene moltiplicato per (1 − 2α). Se |1 − 2α| < 1, cioè 0 < α < 1, θ tende a zero (converge). Se 0 < α < 0,5 il fattore è positivo e θ scende «da un lato solo»; se 0,5 < α < 1 il fattore è negativo e θ oscilla da un lato all'altro ma con ampiezza decrescente; α = 1 oscilla senza mai ridursi; α > 1 oscilla con ampiezza crescente: diverge. Lo verifichiamo:

for alpha in (0.1, 0.5, 0.9, 1.0, 1.1):
    theta, traiettoria = 1.0, [1.0]
    for _ in range(10):
        theta = theta - alpha * 2 * theta        # un passo di discesa su J = theta^2
        traiettoria.append(theta)
    print("alpha = %.1f  primi passi: %s   dopo 10 passi: %.4f" % (alpha, np.round(traiettoria[:5], 3), theta))

alpha = 0.1  primi passi: [1.    0.8   0.64  0.512 0.41 ]   dopo 10 passi: 0.1074
alpha = 0.5  primi passi: [1. 0. 0. 0. 0.]   dopo 10 passi: 0.0000
alpha = 0.9  primi passi: [ 1.    -0.8    0.64  -0.512  0.41 ]   dopo 10 passi: 0.1074
alpha = 1.0  primi passi: [ 1. -1.  1. -1.  1.]   dopo 10 passi: 1.0000
alpha = 1.1  primi passi: [ 1.    -1.2    1.44  -1.728  2.074]   dopo 10 passi: 6.1917

Con α = 0,5 il fattore (1 − 2α) è zero e si arriva al minimo in un passo (caso fortunato dovuto alla perfetta simmetria della parabola). Con 0,9 si oscilla ma si converge; con 1,0 si rimbalza tra +1 e −1 per sempre; con 1,1 l'ampiezza cresce e la discesa esplode. Nella regressione lineare la soglia di stabilità dipende dalla curvatura di J, cioè dal più grande autovalore λ_max della matrice Hessiana H = XᵀX/m: la discesa converge se α < 2/λ_max. Più le feature hanno valori grandi, più λ_max è grande e più il passo ammissibile è piccolo. Questo ci porta alla trappola più comune.

Discesa del gradiente a mano per la regressione lineare

Implementiamo la regola θ ← θ − α∇J sugli stessi dati sintetici della lezione precedente (y = 3 + 2x + rumore) e la confrontiamo con la soluzione esatta:

rng = np.random.default_rng(0)
m = 100
x = rng.uniform(0, 10, m)
y = 3 + 2 * x + rng.normal(0, 2, m)
X = np.column_stack([np.ones(m), x])                       # (m, 2): colonna di uni + x
theta_esatto = np.linalg.solve(X.T @ X, X.T @ y)           # equazione normale, per confronto

def discesa(X, y, alpha, iterazioni):
    theta = np.zeros(X.shape[1])                           # partenza: tutti i parametri a zero
    costi = []
    for _ in range(iterazioni):
        gradiente = X.T @ (X @ theta - y) / len(y)         # (1/m) X^T (X theta - y)
        theta = theta - alpha * gradiente                  # passo contro il gradiente
        costi.append(((X @ theta - y) ** 2).sum() / (2 * len(y)))
    return theta, costi

H = X.T @ X / m
ev = np.linalg.eigvalsh(H)                                 # autovalori (H è simmetrica), in ordine crescente
print("autovalori dell'Hessiana: %s  ->  soglia 2/lambda_max = %.4f, condizionamento %.0f" % (np.round(ev, 4), 2 / ev.max(), ev.max() / ev.min()))
print("soluzione esatta:", np.round(theta_esatto, 4))
for alpha, it in ((0.01, 1000), (0.02, 1000)):
    th, costi = discesa(X, y, alpha, it)
    print("alpha=%.2f, %4d iterazioni -> theta = %s   J: %.2f (inizio) %.4f (fine)" % (alpha, it, np.round(th, 4), costi[0], costi[-1]))
with np.errstate(all="ignore"):
    th, costi = discesa(X, y, 0.06, 200)
print("alpha=0.06 (sopra la soglia): theta = %s   J finale = %.3g" % (np.round(th, 3), costi[-1]))

autovalori dell'Hessiana: [ 0.2296 40.0196]  ->  soglia 2/lambda_max = 0.0500, condizionamento 174
soluzione esatta: [2.9851 1.9739]
alpha=0.01, 1000 iterazioni -> theta = [2.7184 2.0113]   J: 43.14 (inizio) 1.8884 (fine)
alpha=0.02, 1000 iterazioni -> theta = [2.9584 1.9776]   J: 7.18 (inizio) 1.8802 (fine)
alpha=0.06 (sopra la soglia): theta = [-6.55996369e+28 -4.66845879e+29]   J finale = 4.45e+60

Il codice fa un'unica cosa ripetuta: calcola il gradiente con la formula matriciale (una moltiplicazione di Xᵀ per il vettore degli errori, diviso m), sposta θ di α volte il gradiente, e registra il costo per poter vedere la curva di apprendimento. Prima stampa gli autovalori dell'Hessiana (eigvalsh di NumPy, pensata per matrici simmetriche, documentazione) e la soglia 2/λ_max: con le feature x tra 0 e 10, λ_max è grande, e un α sopra la soglia fa esplodere il costo. np.errstate sopprime gli avvisi di overflow dell'ultimo esperimento, che diverge di proposito.

Osservazione importante: con α sotto soglia la discesa si avvicina alla soluzione esatta, ma anche dopo 1000 iterazioni non coincide del tutto (con α = 0,01 l'intercetta vale 2,72 contro 2,99). Il motivo sta nei due autovalori stampati: la curvatura è molto diversa nelle due direzioni, il passo è limitato dalla più ripida (λ_max) e la più piatta (λ_min, circa 0,23) avanza lenta. Il rapporto tra i due, il condizionamento, è il sintomo che la sezione sulla scala delle feature affronta e che più sotto prende un nome e una misura.

spessore x (decimi di mm)    1     2     3     4
ore       y                  1     3     7     8

La trappola della scala delle feature

Immagina di cambiare l'unità di misura di x, per esempio dai metri ai centimetri: i valori diventano cento volte più grandi. La forma di J cambia drasticamente: diventa una conca lunghissima e stretta, ripida in una direzione e quasi piatta nell'altra. Un unico α non può andare bene per entrambe: abbastanza piccolo da non divergere sulla direzione ripida significa lentissimo sulla direzione piatta. Verifichiamolo moltiplicando x per 100:

xr = x * 100                                               # stessa informazione, scala x100
Xr = np.column_stack([np.ones(m), xr])
for alpha in (1e-7, 1e-6, 1e-5):
    with np.errstate(all="ignore"):
        th, costi = discesa(Xr, y, alpha, 1000)
    print("scala originale x100, alpha=%.0e -> theta = %s   J finale = %s" % (alpha, np.round(th, 4), "nan (divergenza)" if np.isnan(costi[-1]) else "%.3f" % costi[-1]))

mu, sd = xr.mean(), xr.std()
Xs = np.column_stack([np.ones(m), (xr - mu) / sd])         # feature standardizzata: media 0, deviazione 1
th_s, costi_s = discesa(Xs, y, 0.1, 100)
print("standardizzata, alpha=0.1, 100 iterazioni -> J finale = %.4f   (minimo teorico %.4f)" % (costi_s[-1], ((Xs @ np.linalg.solve(Xs.T @ Xs, Xs.T @ y) - y) ** 2).sum() / (2 * m)))
# tornare alle unità originali: y = t0 + t1 * (x - mu)/sd  =>  pendenza t1/sd, intercetta t0 - t1*mu/sd
print("in unità originali: pendenza = %.5f (vera 0.02),  intercetta = %.4f (esatta 2.985)" % (th_s[1] / sd, th_s[0] - th_s[1] * mu / sd))

scala originale x100, alpha=1e-07 -> theta = [0.0001 0.0239]   J finale = 2.923
scala originale x100, alpha=1e-06 -> theta = [0.0007 0.0239]   J finale = 2.922
scala originale x100, alpha=1e-05 -> theta = [nan nan]   J finale = nan (divergenza)
standardizzata, alpha=0.1, 100 iterazioni -> J finale = 1.8801   (minimo teorico 1.8801)
in unità originali: pendenza = 0.01974 (vera 0.02),  intercetta = 2.9850 (esatta 2.985)

Cosa fa il codice. Per valori di α piccoli (1e-7, 1e-6) la discesa non diverge, ma dopo 1000 iterazioni il costo è fermo a 2.92, contro il minimo di 1.88: la direzione dell'intercetta non si è mossa quasi per niente. Con α = 1e-5 il passo è già sopra la soglia di stabilità e i valori esplodono fino a nan. Non esiste un α che vada bene. Dopo la standardizzazione (sottrarre la media e dividere per la deviazione standard, per ogni feature) la conca diventa tondeggiante: un α = 0,1 converge al minimo teorico in meno di 100 iterazioni. Le ultime righe mostrano come tornare ai coefficienti nelle unità originali: dato che y = θ₀ + θ₁(x − μ)/σ, la pendenza originale è θ₁/σ e l'intercetta θ₀ − θ₁μ/σ.

Regola operativa: standardizza sempre le feature prima di una discesa del gradiente (con media e deviazione calcolate solo sul training set), e ricorda di applicare la stessa trasformazione ai dati di test.

Tre varianti: batch, stocastica, mini-batch

Il gradiente (1/m) Σ ... somma su tutti gli m esempi: costa un passaggio sull'intero dataset per fare un solo passo. Con milioni di esempi è proibitivo. Il gradiente calcolato su un sottoinsieme è una stima rumorosa ma non distorta di quello completo, quindi si può fare un passo anche con pochi esempi. Le tre varianti differiscono per quanti:

• Batch: tutti gli m esempi per ogni passo. Traiettoria liscia e deterministica, ma un solo passo per passaggio sui dati e costo alto per passo.
• Stocastica (SGD): un solo esempio per passo. Moltissimi passi economici, ma rumorosi: il costo oscilla attorno al minimo senza fermarsi.
• Mini-batch: un piccolo gruppo (tipicamente 16-512 esempi). Compromesso: stima del gradiente meno rumorosa e calcolo efficiente su GPU. È quella usata in pratica nel deep learning.

Un'epoca è un passaggio completo su tutto il training set. Con batch di dimensione b si fanno m/b passi per epoca. Confrontiamo le tre varianti a parità di epoche (10):

def addestra(X, y, dim_batch, alpha, epoche, seed=0):
    rng = np.random.default_rng(seed)
    n = len(y); theta = np.zeros(X.shape[1]); passi = 0; costi = []
    for _ in range(epoche):
        ordine = rng.permutation(n)                        # rimescola ad ogni epoca
        for inizio in range(0, n, dim_batch):
            b = ordine[inizio:inizio + dim_batch]
            g = X[b].T @ (X[b] @ theta - y[b]) / len(b)    # gradiente stimato sul mini-batch
            theta = theta - alpha * g; passi += 1
        costi.append(((X @ theta - y) ** 2).sum() / (2 * n))
    return theta, passi, costi

Xz = np.column_stack([np.ones(m), (x - x.mean()) / x.std()])
minimo = ((Xz @ np.linalg.solve(Xz.T @ Xz, Xz.T @ y) - y) ** 2).sum() / (2 * m)
print("costo minimo teorico: %.4f" % minimo)
for nome, b, a in (("batch (100)", 100, 0.1), ("mini-batch (10)", 10, 0.05), ("stocastica (1)", 1, 0.01)):
    th, passi, costi = addestra(Xz, y, b, a, 10)
    print("%-16s passi: %4d   J dopo 1, 2, 5, 10 epoche: %s" % (nome, passi, np.round([costi[0], costi[1], costi[4], costi[9]], 3)))

costo minimo teorico: 1.8801
batch (100)      passi:   10   J dopo 1, 2, 5, 10 epoche: [93.589 76.164 41.358 15.645]
mini-batch (10)  passi:  100   J dopo 1, 2, 5, 10 epoche: [42.349 16.37   2.55   1.886]
stocastica (1)   passi: 1000   J dopo 1, 2, 5, 10 epoche: [16.991  3.891  1.885  1.881]

Cosa fa il codice. La funzione rimescola gli indici a ogni epoca (rng.permutation) e li scorre a blocchi di dim_batch; per ogni blocco calcola il gradiente solo su quegli esempi e fa un passo. La variante cambia solo dim_batch: 100 = batch, 10 = mini-batch, 1 = stocastica. A parità di epoche, il batch fa solo 10 passi in totale e non ha ancora raggiunto il minimo; mini-batch (100 passi) e stocastica (1000 passi) ci arrivano. I learning rate sono diversi perché passi più rumorosi vogliono α più piccoli. La traiettoria stocastica non si ferma mai del tutto: il costo continua a oscillare lievemente attorno al minimo, ed è per questo che nella pratica si riduce α col tempo (learning rate decay).

Verificare un gradiente: il gradient check

Quando scriverai a mano il gradiente di una funzione più complicata, un errore di segno o una dimenticanza sono facili. La difesa standard è confrontare il gradiente analitico con quello numerico (differenze finite centrali), parametro per parametro:

def costo(theta):
    return ((X @ theta - y) ** 2).sum() / (2 * m)

def gradiente_analitico(theta):
    return X.T @ (X @ theta - y) / m

def gradiente_numerico(theta, eps=1e-6):
    g = np.zeros_like(theta)
    for k in range(len(theta)):
        piu, meno = theta.copy(), theta.copy()
        piu[k] += eps; meno[k] -= eps
        g[k] = (costo(piu) - costo(meno)) / (2 * eps)
    return g

t_prova = np.array([1.0, 0.5])
ga, gn = gradiente_analitico(t_prova), gradiente_numerico(t_prova)
print("analitico:", np.round(ga, 6), " numerico:", np.round(gn, 6))
print("errore relativo: %.2e" % (np.linalg.norm(ga - gn) / np.linalg.norm(ga + gn)))

analitico: [-10.066094 -68.731636]  numerico: [-10.066094 -68.731636]
errore relativo: 5.23e-11

Cosa fa il codice. gradiente_numerico perturba un parametro alla volta in avanti e indietro di una quantità minuscola e misura la variazione del costo. L'errore relativo è dell'ordine di 10⁻¹⁰ o meno: se tra le due versioni c'è una differenza dell'ordine di 10⁻³ o più, c'è un bug nel gradiente analitico. Riutilizzerai questa tecnica per controllare la backpropagation.

La libreria: PyTorch calcola il gradiente per te

Scrivere il gradiente a mano è istruttivo, ma per modelli complessi si usa la differenziazione automatica: la libreria registra le operazioni eseguite (un «grafo» delle operazioni) e ricava il gradiente applicando la regola della catena all'indietro. La stessa regressione in PyTorch. Il listato riusa Xz, y e minimo dei blocchi precedenti; è stato eseguito due volte con PyTorch 2.12.0 su CPU (CUDA_VISIBLE_DEVICES="", 2 thread) con uscita identica.

import torch

Xt = torch.tensor(Xz, dtype=torch.float32)
yt = torch.tensor(y, dtype=torch.float32)
theta_t = torch.zeros(2, requires_grad=True)               # requires_grad: "tieni traccia delle operazioni su questo tensore"
ottimizzatore = torch.optim.SGD([theta_t], lr=0.1)

for epoca in range(100):
    ottimizzatore.zero_grad()                              # azzera i gradienti accumulati
    costo_t = ((Xt @ theta_t - yt) ** 2).sum() / (2 * len(yt))
    costo_t.backward()                                     # calcola d(costo)/d(theta) per ogni parametro
    ottimizzatore.step()                                   # theta <- theta - lr * gradiente

print("PyTorch: theta = %s   J = %.4f   (a mano: J minimo %.4f)" % (np.round(theta_t.detach().numpy(), 4), costo_t.item(), minimo))

PyTorch: theta = [13.8072  5.9826]   J = 1.8801   (a mano: J minimo 1.8801)

Il risultato coincide con la soluzione esatta calcolata con NumPy sugli stessi dati standardizzati (equazione normale, np.linalg.solve), θ = [13.8075, 5.9827] e J minimo = 1.8801, a meno della precisione di float32: la differenza sul quarto decimale (13.8072 contro 13.8075) è proprio quella. Il costo stampato è quello calcolato nell'ultima epoca, prima dell'ultimo aggiornamento.

Cosa fa il codice. requires_grad=True dice a PyTorch di costruire il grafo delle operazioni. backward() percorre il grafo all'indietro e deposita in theta_t.grad il gradiente del costo. ottimizzatore.step() esegue la regola θ ← θ − α·grad (con i valori predefiniti di torch.optim.SGD, cioè momentum 0, è esattamente la regola di questa lezione; il learning rate predefinito è 0,001, per questo il listato lo imposta a 0,1). zero_grad() è necessario perché PyTorch accumula i gradienti: senza azzerarli, a ogni epoca si sommerebbero a quelli precedenti (una delle trappole più frequenti). I valori di θ sono nelle unità della feature standardizzata, come nei blocchi precedenti. Il comportamento di zero_grad() e di backward() è descritto nelle note sull'autograd e nella pagina di torch.optim.SGD della documentazione ufficiale; cosa fa esattamente backward() lo derivi per intero in backpropagation derivata passaggio per passaggio.

Trade-off e trappole

• Learning rate: troppo grande diverge, troppo piccolo non arriva. Si esplora su scala logaritmica (0,001; 0,01; 0,1...) e si guarda la curva del costo: deve scendere in modo regolare. Se sale o diventa nan, dimezza o dividi per dieci α.
• Non standardizzare le feature: vedi sopra. È la causa più frequente di discese che «non funzionano».
• Funzioni non convesse: la regressione lineare ha un solo minimo; le reti neurali hanno una superficie di costo piena di avvallamenti, selle e zone piatte. La discesa del gradiente trova un minimo locale (o un punto stazionario), non garantisce il migliore. In pratica funziona sorprendentemente bene, ma il risultato dipende dall'inizializzazione.
• SGD senza rimescolamento: se i dati sono ordinati (per esempio per classe), i mini-batch consecutivi sono tutti della stessa classe e il gradiente è fuorviante. Rimescola a ogni epoca.
• Criterio di arresto: un numero fisso di iterazioni è semplice; meglio fermarsi quando il costo su un set di validazione smette di calare (early stopping).
• Gradienti non azzerati (PyTorch) e gradienti calcolati a mano senza verifica: due bug silenziosi; il secondo si previene con il gradient check.

Le sei voci non sono indipendenti. Learning rate, scala delle feature e criterio di arresto rispondono alla stessa domanda, cioè se il passo è adatto alla forma della superficie e quando fermarsi; la non convessità decide dove si arriva; il rimescolamento decide se il gradiente stimato sul mini-batch è affidabile. Le ultime due voci sono bug silenziosi: nessun errore, solo numeri sbagliati, ed è per questo che il gradient check e il controllo della curva del costo vanno fatti prima di fidarsi del risultato.

Il numero di condizionamento: perché la discesa del gradiente diventa lenta

Nella sezione sulla scala delle feature abbiamo visto che feature con scale molto diverse rallentano la discesa. Ora diamo a quel fenomeno un numero. Per una funzione di costo quadratica, come l'MSE della regressione lineare, la curvatura in ogni direzione è un autovalore della matrice Hessiana H (la matrice delle derivate seconde). Il learning rate più grande che non diverge è 2/λ_max, dove λ_max è l'autovalore più grande; ma la direzione con l'autovalore più piccolo λ_min progredisce a passo λ_min · η, quindi lentissimo. Il rapporto κ = λ_max / λ_min è il numero di condizionamento, e il numero di iterazioni necessarie cresce proporzionalmente a κ.

Intuizione: una valle lunga e stretta, ripidissima sui fianchi e quasi piatta lungo il fondo; il learning rate deve essere piccolo per non rimbalzare da un fianco all'altro, e per questo si avanza lungo il fondo a passi minuscoli.

Il momentum (o «metodo della palla pesante») rimedia accumulando una velocità: v ← β·v + g, w ← w − η·v (è la forma di torch.optim.SGD con momentum=β e dampening=0, come descritto nella documentazione). Le componenti del gradiente che oscillano da un passo all'altro (le pareti della valle) si annullano nella somma, mentre quelle che puntano sempre nella stessa direzione (il fondo) si sommano. Con la scelta ottimale dei parametri per una quadratica, η = 4/(√κ + 1)² e β = ((√κ − 1)/(√κ + 1))², il numero di iterazioni cresce come √κ invece che come κ. Misuriamolo su una funzione con due autovalori, 1 e κ, contando i passi per portare la norma del punto sotto 10⁻⁶.

def steps_to_converge(kappa, lr, momentum=0.0, tol=1e-6, max_steps=200000):
    curvature = np.array([1.0, kappa])                                    # la funzione è 0,5·(x₁² + κ·x₂²)
    point, velocity = np.array([1.0, 1.0]), np.zeros(2)
    for step in range(1, max_steps + 1):
        velocity = momentum * velocity + curvature * point                # il gradiente è curvatura · punto
        point = point - lr * velocity
        if np.linalg.norm(point) < tol:
            return step
    return None

print("kappa    discesa semplice (lr=1/kappa)   con momentum ottimale   rapporto")
for kappa in (10, 100, 1000):
    plain = steps_to_converge(kappa, 1.0 / kappa)
    root = np.sqrt(kappa)
    heavy = steps_to_converge(kappa, 4 / (root + 1) ** 2, ((root - 1) / (root + 1)) ** 2)
    print("%-8d  %-30d  %-21d  %.1fx" % (kappa, plain, heavy, plain / heavy))

kappa    discesa semplice (lr=1/kappa)   con momentum ottimale   rapporto
10        132                             27                     4.9x
100       1375                            95                     14.5x
1000      13809                           321                    43.0x

Cosa fa il codice. La funzione conta i passi necessari a portare il punto, partito da (1, 1), sotto una norma di 10⁻⁶, su una quadratica con curvature 1 e κ. La discesa semplice, con learning rate 1/κ (metà del massimo stabile 2/κ), richiede 132 passi con κ = 10, 1 375 con κ = 100 e 13 809 con κ = 1 000: cresce quasi esattamente di dieci volte ogni volta che κ si moltiplica per dieci, come prevede la teoria (il numero di passi è circa κ · ln(10⁶) ≈ 13,8 · κ). Con il momentum ottimale i passi sono 27, 95 e 321: moltiplicare κ per dieci moltiplica i passi per circa tre (√10 ≈ 3,2).

Il risparmio cresce con la difficoltà del problema: 4,9 volte meno passi a κ = 10, 14,5 a κ = 100, 43,0 a κ = 1 000.

Due avvertenze oneste. Primo, i parametri «ottimali» presuppongono di conoscere gli autovalori, cosa che in una rete neurale non è vera: nella pratica si fissa β = 0,9 e si regola il learning rate, e il guadagno è più modesto. Secondo, questa è una quadratica con curvature allineate agli assi, il caso più semplice: le reti neurali hanno funzioni di costo non convesse, e per esse il momentum è un'euristica ben collaudata, non un teorema.

Il metodo di Newton: un solo passo per la regressione lineare

La discesa del gradiente usa solo la pendenza. Il metodo di Newton usa anche la curvatura: w ← w − H⁻¹ g, cioè corregge il gradiente con l'inverso dell'Hessiana, e per questo «vede» la forma della valle e va dritto al fondo. Per una funzione quadratica, come l'MSE della regressione lineare, un solo passo di Newton raggiunge il minimo esatto (è l'equazione normale in un'altra veste). Lo verifichiamo su dati con due feature di scala molto diversa (deviazione 1 e 30), confrontando con la discesa del gradiente al massimo learning rate stabile.

rng_g = np.random.default_rng(0)
n_g = 200
X_g = np.column_stack([np.ones(n_g), rng_g.normal(0, 1, n_g), rng_g.normal(0, 30, n_g)])
y_g = X_g @ np.array([1.0, 2.0, -0.05]) + rng_g.normal(0, 0.5, n_g)

hessian = X_g.T @ X_g / n_g
w_star = np.linalg.solve(X_g.T @ X_g, X_g.T @ y_g)                        # soluzione esatta
eigenvalues = np.linalg.eigvalsh(hessian)
print("autovalori dell'Hessiana: %s   numero di condizionamento %.0f" % (np.round(eigenvalues, 3), eigenvalues.max() / eigenvalues.min()))

w = np.zeros(3)
for iteration in range(1, 100001):
    w = w - (1.0 / eigenvalues.max()) * (X_g.T @ (X_g @ w - y_g) / n_g)
    if np.abs(w - w_star).max() < 1e-6:
        break
print("discesa del gradiente: %d iterazioni per arrivare entro 1e-6 dalla soluzione" % iteration)

w_newton = np.zeros(3) - np.linalg.solve(hessian, X_g.T @ (X_g @ np.zeros(3) - y_g) / n_g)
print("Newton: 1 iterazione, errore massimo %.1e" % np.abs(w_newton - w_star).max())

autovalori dell'Hessiana: [9.1900e-01 9.9400e-01 9.5315e+02]   numero di condizionamento 1037
discesa del gradiente: 15002 iterazioni per arrivare entro 1e-6 dalla soluzione
Newton: 1 iterazione, errore massimo 5.8e-15

Cosa fa il codice. Costruisce 200 punti con tre colonne (un termine costante, una feature di deviazione 1 e una di deviazione 30) e un bersaglio lineare con rumore 0,5. L'Hessiana dell'errore (con la convenzione del mezzo MSE, gradiente XᵀXw/n) ha autovalori 0,919, 0,994 e 953,15: il numero di condizionamento è 1 037, quasi interamente dovuto alla feature di scala 30, il cui quadrato (900) compare nella curvatura. La discesa del gradiente, anche con il massimo learning rate sicuro 1/λ_max, impiega 15 002 iterazioni per arrivare a meno di 10⁻⁶ dalla soluzione esatta (la stessa conta degli esperimenti di sopra: circa 14 · κ). Newton, in una sola iterazione, trova la soluzione con un errore massimo di 5,8·10⁻¹⁵, cioè all'arrotondamento.

Perché allora non si usa Newton ovunque? Perché richiede di formare e invertire una matrice p × p: il costo è O(p³) e la memoria O(p²). Per una regressione con tre coefficienti è nulla; per una rete con cento milioni di parametri è impossibile. Questo è il trade-off di fondo dell'ottimizzazione: i metodi del secondo ordine sono potenti ma costosi, quelli del primo ordine (con momentum, o con scale adattive come Adam) sono scalabili. E la standardizzazione delle feature, che costa quasi nulla, è un modo di avvicinare il problema al caso ben condizionato dove anche il primo ordine funziona.

Rumore della SGD e decadimento del learning rate

x_n = rng_g.normal(0, 1, 1000)
y_n = 3 * x_n + rng_g.normal(0, 1, 1000)
w_best = (x_n @ y_n) / (x_n @ x_n)                                        # minimo esatto per il modello y = w·x

def sgd_final_error(schedule, seed, epochs=30):
    order_rng = np.random.default_rng(seed)
    w, last = 0.0, []
    for epoch in range(epochs):
        for i in order_rng.permutation(1000):
            w -= schedule(epoch) * (w * x_n[i] - y_n[i]) * x_n[i]         # un passo per esempio
        last.append(w)
    return np.abs(np.array(last[-5:]) - w_best).mean()                    # errore medio sulle ultime 5 epoche

print("soluzione esatta w = %.4f   (il valore vero è 3)" % w_best)
for name, schedule in (("learning rate fisso 0,05", lambda e: 0.05), ("decadimento 0,05/(1+epoca)", lambda e: 0.05 / (1 + e))):
    errors = [sgd_final_error(schedule, seed) for seed in range(10)]
    print("%-34s  distanza media dal minimo, 10 semi: %.4f" % (name, np.mean(errors)))

soluzione esatta w = 3.0497   (il valore vero è 3)
learning rate fisso 0,05            distanza media dal minimo, 10 semi: 0.1259
decadimento 0,05/(1+epoca)          distanza media dal minimo, 10 semi: 0.0123

Cosa fa il codice. Addestra un modello a un solo parametro, y = w·x, con la discesa stocastica (un esempio per passo, 30 epoche su 1 000 punti), e misura la distanza media dalla soluzione esatta (3,0497, vicina al valore vero 3) sulle ultime cinque epoche, su dieci ordinamenti casuali diversi.

Con learning rate fisso 0,05 la distanza media è 0,1259; con il learning rate che decade come 0,05/(1 + epoca) è 0,0123, circa dieci volte più piccola. Il motivo è la natura della SGD: il gradiente di un singolo esempio è una stima rumorosa del gradiente vero, e con un passo fisso il parametro continua a oscillare attorno al minimo con un'ampiezza proporzionale al learning rate, senza mai fermarsi. Se il passo diminuisce, l'oscillazione si restringe. È lo stesso compromesso visto per la discesa a lotti: all'inizio un passo grande per avvicinarsi in fretta, alla fine un passo piccolo per rifinire. Nelle reti reali queste idee prendono il nome di learning rate schedule (decadimento a gradini, coseno, riscaldamento iniziale). Nota la limitazione dell'esperimento: un solo problema a un parametro, e 10 semi; il rapporto di dieci volte non va preso come una costante, ma come l'ordine di grandezza dell'effetto in questo caso.

Quando la funzione non è convessa: il punto di arrivo dipende dalla partenza

f = lambda x: x ** 4 - 3 * x ** 2 + x                                     # due minimi locali
df = lambda x: 4 * x ** 3 - 6 * x + 1
print("partenza    arrivo    valore della funzione")
for start in (-2.0, -0.5, 0.5, 2.0):
    x = start
    for _ in range(2000):
        x -= 0.01 * df(x)
    print("%+.1f       %+.4f   %+.4f" % (start, x, f(x)))

partenza    arrivo    valore della funzione
-2.0       -1.3008   -3.5139
-0.5       -1.3008   -3.5139
+0.5       +1.1309   -1.0702
+2.0       +1.1309   -1.0702

Cosa fa il codice. Applica la discesa del gradiente (2 000 passi di 0,01) alla funzione f(x) = x⁴ − 3x² + x, che ha due valli. Le quattro partenze danno due soli arrivi: partendo da −2 o da −0,5 si finisce nel minimo di sinistra, x = −1,3008, dove la funzione vale −3,5139; partendo da +0,5 o da +2 si finisce in quello di destra, x = +1,1309, dove vale −1,0702. Entrambi sono punti in cui il gradiente è zero e la discesa si ferma soddisfatta, ma solo il primo è il minimo globale: l'altro è un minimo locale, con un costo di 2,44 unità più alto. L'algoritmo non ha modo di accorgersene, perché vede solo la pendenza qui e ora.

Questo è esattamente ciò che cambia quando si passa dalla regressione lineare (funzione convessa, un solo minimo, la partenza non conta) alle reti neurali (funzione non convessa, risultato dipendente dall'inizializzazione e dal rumore della SGD). Le reti molto grandi sono meno afflitte dal problema di quanto suggerisca l'intuizione a una dimensione, perché in molte dimensioni i minimi locali davvero cattivi sembrano rari rispetto ai punti di sella; ma è un'osservazione empirica, non una garanzia, e per questo si addestra spesso con più semi e si riporta la variabilità, come facciamo in queste lezioni.

Per continuare da solo

Tre prove da fare senza soluzione, ognuna sugli esperimenti già eseguiti qui sopra. Prova A. Standardizza le due feature dell'esperimento di Newton (media 0, deviazione 1) e riconta le iterazioni della discesa del gradiente: il numero di condizionamento cala? di quanto? Che cosa dice sul rapporto tra preprocessing e ottimizzazione? Prova B. Nell'esperimento di momentum prova valori di β non ottimali (0,5; 0,99): per quali la convergenza peggiora rispetto alla discesa semplice? Prova C. Nell'esperimento non convesso cerca, a mano o con una griglia di partenze, il punto di separazione tra i due bacini di attrazione: dove si trova rispetto ai due minimi?

Nella stessa serie rimandano a questa lezione anche: Analisi matematica 1: limiti e derivate, Analisi matematica 2: gradiente, integrali, discesa del gradiente, Support Vector Machine: il margine massimo e il trucco del kernel, Boosting: AdaBoost, Gradient Boosting e cenni su XGBoost, Regressione logistica: dalla retta alla probabilità, Il percettrone e le funzioni di attivazione: a cosa serve la non linearità, Reti neurali ricorrenti (RNN): dare una memoria alla rete, Backpropagation through time: gradiente che svanisce e che esplode.

Fonti

• PyTorch, Autograd mechanics - il fatto che i gradienti si accumulano nel campo .grad dei tensori foglia con requires_grad=True, consultata il 2026-11-12
• PyTorch, torch.optim.SGD - parametri predefiniti (lr 0,001, momentum 0) e la regola di aggiornamento con momentum, consultata il 2026-11-12
• NumPy, numpy.linalg.eigvalsh - gli autovalori di una matrice simmetrica, usata per λ_max e λ_min, consultata il 2026-11-12

Riferimenti storici citati a memoria e non riaperti nel corso della verifica (la pagina dell'editore di Robbins e Monro ha rifiutato un client automatico): A.-L. Cauchy, Méthode générale pour la résolution des systèmes d'équations simultanées, 1847, per la discesa del gradiente; H. Robbins e S. Monro, A Stochastic Approximation Method, Annals of Mathematical Statistics, 1951, per la versione stocastica.

Nella prossima lezione estendi la regressione a più feature e a relazioni curve in regressione multivariata, polinomiale e normalizzazione delle feature, e vedi come lo stesso gradiente continua a funzionare.