57 - Backpropagation through time: gradiente che svanisce e che esplode

Una RNN si addestra come qualunque rete: si calcola l'errore, lo si porta all'indietro con la regola della catena e si aggiornano i pesi. Il trucco è che la rete, srotolata nel tempo, è profonda quanto la sequenza è lunga, quindi il gradiente deve attraversare T passi uguali. Il gradiente è un prodotto di T matrici: se queste matrici sono in qualche senso «più piccole di 1» il prodotto tende a zero (gradiente che svanisce: la rete non impara le dipendenze lontane); se sono «più grandi di 1» può crescere senza controllo (gradiente che esplode: l'addestramento diventa instabile). Questo articolo deriva la formula, la verifica numericamente, misura entrambi i fenomeni e mostra il rimedio per l'esplosione, il gradient clipping, insieme a ciò che non cura. È il seguito di Reti neurali ricorrenti (RNN).

Le parole che servono

La perdita L è un solo numero che misura quanto la rete ha sbagliato, e addestrare significa farla scendere. La derivata dice quanto cambia un numero se ne ritocco un altro di pochissimo: se ∂L/∂w = 3, alzare w di 0.01 alza la perdita di circa 0.03. Il gradiente è l'insieme di tutte le derivate della perdita rispetto a tutti i pesi, e l'aggiornamento è w ← w − η · ∂L/∂w, dove η è il learning rate, la lunghezza del passo.

La regola della catena dice che se A influenza B e B influenza C, l'effetto di A su C è il prodotto dei due effetti: in una catena di ingranaggi, se il primo gira di un giro e il secondo ne fa due, e il terzo ne fa tre per ogni giro del secondo, il terzo fa sei giri; i fattori si moltiplicano. Quando ingresso e uscita sono vettori la «derivata» è una matrice, il Jacobiano: la riga i, colonna j dice quanto cambia la componente i dell'uscita al variare della componente j dell'ingresso; qui ne serve uno per passo, ∂h_t/∂h_{t-1}. La norma ‖v‖ è la «lunghezza» di un vettore (o la «grandezza» di una matrice) e permette di dire con un solo numero se un gradiente è grande o piccolo. La backpropagation è l'algoritmo che applica la regola della catena all'indietro, dalla perdita ai pesi, riusando i calcoli intermedi.

L'analogia: il telefono senza fili

Dieci persone in fila si passano un messaggio sussurrandolo. Se ognuna lo attenua un po', all'ultima arriva silenzio: è il gradiente che svanisce. Se ognuna lo amplifica un po', all'ultima arriva un urlo distorto: è il gradiente che esplode. Nella rete il «messaggio» è il segnale di errore che deve raggiungere un passo lontano nel tempo per dirgli «correggi il tuo contributo», e ogni passo lo moltiplica per una matrice.

La matematica: derivare la BPTT

Prendiamo il caso più semplice, many-to-one, con una sola perdita alla fine. Con a_t = W_xh x_t + W_hh h_{t-1} e h_t = tanh(a_t) la perdita è:

L = ½ (w · h_T − y)²

Il peso W_hh compare a ogni passo. Per la regola della catena il suo effetto sulla perdita è la somma dei contributi di ogni passo t:

∂L/∂W_hh = Σ_{t=1..T}  (∂L/∂h_T) · (∂h_T/∂h_t) · (∂h_t/∂W_hh)|diretto

Qui «diretto» significa la derivata di h_t rispetto a W_hh tenendo h_{t-1} fermo. Il pezzo critico è ∂h_T/∂h_t, cioè quanto una piccola variazione dello stato al tempo t arriva fino alla fine. Siccome h_j dipende da h_{j-1}, è un prodotto di Jacobiani:

∂h_T/∂h_t = Π_{j=t+1..T}  ∂h_j/∂h_{j-1}
∂h_j/∂h_{j-1} = diag(1 − h_j²) · W_hh        (derivata della tanh: 1 − tanh²)

Prendiamo le norme. La derivata della tanh vale al massimo 1, e vale quasi 0 quando la tanh è satura, cioè quando |h| è vicino a 1. Detto γ il massimo della derivata (qui 1) e σ il massimo valore singolare di W_hh, ogni fattore ha norma al più γσ e quindi:

‖∂h_T/∂h_t‖ ≤ (γ · σ)^(T − t)

Da qui due regimi, analizzati da Bengio, Simard e Frasconi (1994) e in forma pulita da Pascanu, Mikolov e Bengio (2013). Se γσ < 1 il prodotto decade esponenzialmente con la distanza T − t, e il gradiente dei passi lontani è praticamente zero: è una garanzia, succede sempre. Se σ è grande il prodotto può crescere esponenzialmente, ma è una condizione necessaria e non sufficiente, perché dipende dalle direzioni e dalla saturazione della tanh.

L'algoritmo si ricava leggendo la formula da destra a sinistra. Partendo da dh = ∂L/∂h_T, a ogni passo all'indietro si passa dalla tanh (da = dh · (1 − h²)), si accumula il contributo a W_hh (da · h_{t-1}ᵀ) e si manda il gradiente al passo precedente (dh = W_hhᵀ · da). Questa è la backpropagation through time (BPTT).

Quando c'è una perdita a ogni passo

Finora la perdita era una sola, alla fine. Nei compiti many-to-many (etichettare ogni parola, prevedere ogni valore successivo) ogni passo ha la sua perdita L_t e la totale è L = L_1 + … + L_T. La derivata rispetto a W_hh diventa una doppia somma:

∂L/∂W_hh = Σ_{t=1..T}  Σ_{k=1..t}  (∂L_t/∂h_t) · (∂h_t/∂h_k) · (∂h_k/∂W_hh)|diretto

L'errore commesso al tempo t manda un segnale indietro verso tutti i passi k ≤ t, e ogni passo k contribuisce al gradiente del peso condiviso; il termine ∂h_t/∂h_k è di nuovo il prodotto di Jacobiani, con t − k fattori. Ne seguono due conseguenze pratiche. I contributi «vicini» (t − k piccolo) arrivano forti e quelli lontani deboli, quindi la rete impara per prime le dipendenze corte e per ultime, o mai, quelle lunghe. Il costo della BPTT completa cresce con T, e anche la memoria, perché ogni stato intermedio va conservato fino al passaggio all'indietro.

Perché il limite sulla norma è vero

La disuguaglianza ‖∂h_T/∂h_t‖ ≤ (γσ)^(T−t) si dimostra in due righe. Per ogni matrice A e vettore v vale ‖A v‖ ≤ ‖A‖ · ‖v‖, dove ‖A‖ è la norma spettrale, il massimo allungamento che A può produrre; ne segue che la norma di un prodotto è al massimo il prodotto delle norme, ‖A B‖ ≤ ‖A‖ · ‖B‖. Ogni Jacobiano è diag(1 − h²) · W_hh: la diagonale ha norma al massimo γ (il valore più grande di 1 − h², cioè 1) e W_hh ha norma σ.

Quindi ogni fattore ha norma al più γσ, e un prodotto di T − t fattori ha norma al più (γσ)^(T−t). Se γσ < 1 questo tende a zero esponenzialmente, ed è una garanzia matematica e non una tendenza. Se invece γσ > 1 la disuguaglianza non dice nulla, perché è solo un limite superiore e il prodotto vero può restare piccolo: per questo l'esplosione è una condizione necessaria e non sufficiente.

Nel disegno le frecce rosse vanno da destra a sinistra: il segnale di errore nasce nella perdita e deve arrivare a h1 passando per quattro moltiplicazioni. Se ogni fattore vale 0.5, all'arrivo il segnale è moltiplicato per 0.5⁴ = 0.0625; con cento passi diventerebbe 0.5¹⁰⁰ ≈ 8·10⁻³¹, un numero che il computer rappresenta senza problemi ma che per l'ottimizzatore equivale a zero.

Codice a mano: BPTT con verifica numerica

Il modo più sicuro per sapere se una derivazione è giusta è confrontarla con le differenze finite, (L(w+ε) − L(w−ε)) / 2ε. Il listato implementa la BPTT per la rete many-to-one con tre neuroni e sei passi, poi ritocca ogni elemento di W_hh di ±10⁻⁶ e confronta.

import numpy as np

rng = np.random.default_rng(1)
n, d, T = 3, 2, 6
Wxh, Whh = rng.normal(0, .6, (n, d)), rng.normal(0, .6, (n, n))
w_out = rng.normal(0, .6, n)
xs, y = rng.normal(size=(T, d)), 0.7

def forward(Whh):
    hs = [np.zeros(n)]
    for x in xs:
        hs.append(np.tanh(Wxh @ x + Whh @ hs[-1]))
    loss = 0.5 * (w_out @ hs[-1] - y) ** 2
    return loss, hs

def bptt(Whh):
    loss, hs = forward(Whh)
    dh = (w_out @ hs[-1] - y) * w_out           # dL/dh_T
    dWhh = np.zeros_like(Whh)
    for t in range(T, 0, -1):                   # walk BACKWARDS in time
        da = dh * (1 - hs[t] ** 2)              # through the tanh
        dWhh += np.outer(da, hs[t - 1])         # this step's contribution to dL/dWhh
        dh = Whh.T @ da                         # pass the gradient to h_{t-1}
    return dWhh

# Check against finite differences
analytic = bptt(Whh)
numeric = np.zeros_like(Whh)
eps = 1e-6
for i in range(n):
    for j in range(n):
        P, M = Whh.copy(), Whh.copy()
        P[i, j] += eps; M[i, j] -= eps
        numeric[i, j] = (forward(P)[0] - forward(M)[0]) / (2 * eps)
print(np.abs(analytic - numeric).max())    # ~1e-10: the manual BPTT is correct

# Output:
# 2.942041055220557e-11

forward è il passaggio in avanti già visto, con in più la perdita finale ½(w·h_T − y)², e restituisce tutti gli stati perché il passaggio indietro ne ha bisogno. bptt parte dalla derivata della perdita rispetto all'ultimo stato (dh) e cammina da t = T fino a t = 1. Dentro il ciclo ci sono le tre righe della derivazione: da = dh * (1 - hs[t]**2) attraversa la tanh, dWhh += np.outer(da, hs[t-1]) accumula il contributo di questo passo al peso condiviso (si somma perché lo stesso peso è usato a ogni passo; con = al posto di += il gradiente sarebbe quello del solo primo passo all'indietro) e dh = Whh.T @ da manda il gradiente al passo precedente.

L'errore massimo è circa 3·10⁻¹¹: la BPTT scritta a mano coincide con la derivata numerica, quindi la derivazione è giusta. Il ciclo che scorre t all'indietro fa capire anche il costo: la memoria cresce con T perché bisogna conservare tutti gli stati hs.

Il caso più piccolo: un solo neurone

Prima delle matrici, un neurone solo senza ingresso, h_t = tanh(w · h_{t−1}). La derivata di h_T rispetto a h_0 è il prodotto dei fattori w · (1 − h_t²), uno per passo: la regola della catena in versione numerica. La funzione scalar_gradient moltiplica un fattore per iterazione e il ciclo finale la prova con quattro pesi a 1, 5, 10 e 20 passi.

import numpy as np

# One-neuron RNN without input:  h_t = tanh(w * h_{t-1}).
# d h_T / d h_0 is the product of the per-step factors  w * (1 - h_t^2).
def scalar_gradient(w, steps, h0=0.5):
    h, grad = h0, 1.0
    for _ in range(steps):
        h = np.tanh(w * h)
        grad *= w * (1 - h ** 2)         # chain rule: multiply one factor per step
    return grad

for w in (0.5, 0.9, 1.1, 2.0):
    print(f"w={w}:", [f"{scalar_gradient(w, s):.2e}" for s in (1, 5, 10, 20)])

# Output:
# w=0.5: ['4.70e-01', '2.88e-02', '9.00e-04', '8.79e-07']
# w=0.9: ['7.40e-01', '3.30e-01', '1.66e-01', '5.42e-02']
# w=1.1: ['8.24e-01', '3.79e-01', '1.43e-01', '2.01e-02']
# w=2.0: ['8.40e-01', '1.68e-03', '2.14e-07', '3.48e-15']

Un risultato che sorprende: anche con w = 2.0, che «sembra» amplificare, il gradiente svanisce (da 0.84 a 3.5·10⁻¹⁵). La ragione è la saturazione: con w grande lo stato viene spinto verso ±1, dove 1 − h² è quasi zero, e il fattore w(1 − h²) scende sotto 1. In questo esperimento a un neurone solo con tanh l'esplosione non compare per nessuno dei quattro pesi; per esplodere servono matrici con direzioni favorevoli (W_hh grande in qualche direzione mentre lo stato non è saturo). È un motivo per cui nella pratica il vanishing si vede sempre, mentre l'exploding compare a tratti, come picchi improvvisi.

Vedere il decadimento: la norma del prodotto di Jacobiani

Misuriamo direttamente ‖∂h_T/∂h_0‖ per una rete ricorrente con matrice casuale riscalata in modo che l'autovalore più grande in modulo valga 0.9 oppure 1.1, sia lineare sia con tanh. jacobian_norm costruisce W e la divide per il suo autovalore più grande in modulo (il «raggio spettrale»), così che valga esattamente radius: 0.9 contrae leggermente, 1.1 espande leggermente. Poi simula T passi e moltiplica via via il Jacobiano del passo (diag(pendenza) @ W) nel prodotto J; con activation="linear" togliamo la tanh per isolare l'effetto della sola matrice. Alla fine restituisce np.linalg.norm(J, 2), la norma 2, cioè il massimo fattore per cui J amplifica un vettore: la misura giusta per «quanto segnale passa».

import numpy as np

def jacobian_norm(T, n=16, radius=1.0, activation="tanh", seed=0):
    """Spectral norm of d h_T / d h_0 for a recurrent net h_t = f(W h_{t-1})."""
    rng = np.random.default_rng(seed)
    W = rng.normal(size=(n, n))
    W = radius * W / max(abs(np.linalg.eigvals(W)))     # rescale: largest |eigenvalue| = radius
    h = rng.normal(size=n) * 0.5
    J = np.eye(n)                                        # running product of Jacobians
    for _ in range(T):
        h = np.tanh(W @ h) if activation == "tanh" else W @ h
        slope = 1 - h**2 if activation == "tanh" else np.ones(n)   # f'(a)
        J = np.diag(slope) @ W @ J                       # chain rule, one step
    return np.linalg.norm(J, 2)

for act in ("linear", "tanh"):
    for radius in (0.9, 1.1):
        vals = [f"{jacobian_norm(T, radius=radius, activation=act):.1e}" for T in (1, 10, 50, 100)]
        print(f"{act:6s} radius={radius}", vals)

# Output:
# linear radius=0.9 ['1.4e+00', '4.4e-01', '6.6e-03', '3.4e-05']
# linear radius=1.1 ['1.7e+00', '3.3e+00', '1.5e+02', '1.8e+04']
# tanh   radius=0.9 ['1.2e+00', '3.6e-01', '5.3e-03', '2.7e-05']
# tanh   radius=1.1 ['1.5e+00', '2.0e+00', '1.4e-02', '4.2e-07']

Nel caso lineare un raggio spettrale appena sotto 1 fa sparire il gradiente in 100 passi (da 1.4 a 3.4·10⁻⁵), e appena sopra 1 lo fa esplodere (1.8·10⁴). Con la tanh succede una cosa istruttiva: anche con raggio 1.1 il gradiente prima cresce (2.0 a T=10) e poi svanisce (4.2·10⁻⁷ a T=100), perché gli stati si saturano verso ±1 e la derivata 1 − h² si avvicina a zero. La tanh frena l'esplosione, ma negli intervalli in cui lo stato non è saturo il gradiente può comunque salire. È un solo seme, quindi conta l'ordine di grandezza e non la cifra esatta. In pratica si osservano entrambe le cose: picchi di gradiente alternati a lunghi tratti piatti.

Gli autovalori, e perché l'inizializzazione ortogonale aiuta

Nel caso lineare il conto è esatto. Se h_t = W h_{t−1}, allora h_T = W^T h_0. Una matrice «normale» si può scrivere con i suoi autovalori λ_i (i fattori per cui stira alcune direzioni speciali): dopo T passi la componente lungo la direzione i è moltiplicata per λ_i^T. Con |λ| = 0.9 si ha 0.9^100 = 2.7·10⁻⁵, con |λ| = 1.1 si ha 1.1^100 = 1.4·10⁴: due numeri solo apparentemente vicini a 1 producono nove ordini di grandezza di differenza, perché T è un esponente.

Il rimedio suggerito è scegliere W_hh in modo che tutti i suoi «fattori di stiramento» valgano esattamente 1. Una matrice ortogonale (QᵀQ = I) conserva la lunghezza di ogni vettore. Il listato lo verifica: np.linalg.qr produce una Q ortogonale da una matrice casuale (ogni matrice quadrata si scompone in una ortogonale per una triangolare), si controlla che ‖Qx‖ = ‖x‖, si moltiplica mille volte per Q e si confronta con una matrice casuale a raggio spettrale 1.

import numpy as np

rng = np.random.default_rng(0)
n = 16

# 1) the eigenvalue picture, in the LINEAR case: h_T = W^T h_0
W = rng.normal(size=(n, n))
W = 0.9 * W / max(abs(np.linalg.eigvals(W)))
print("0.9^100 =", 0.9 ** 100, "   1.1^100 =", 1.1 ** 100)

# 2) an orthogonal matrix keeps the length of every vector: |Q x| = |x|
Q, _ = np.linalg.qr(rng.normal(size=(n, n)))          # QR gives an orthogonal Q
x = rng.normal(size=n)
print("|x| =", np.linalg.norm(x).round(6), "  |Qx| =", np.linalg.norm(Q @ x).round(6))

# 3) the product of Jacobians of a LINEAR recurrent net with orthogonal W: norm stays 1 forever
J = np.eye(n)
for _ in range(1000):
    J = Q @ J
print("norm of Q^1000:", np.linalg.norm(J, 2).round(6))

# versus a random matrix rescaled to spectral radius 1.0 (not normal -> norms can still drift)
R = rng.normal(size=(n, n)); R = R / max(abs(np.linalg.eigvals(R)))
J = np.eye(n)
for T in (10, 100, 1000):
    J = np.linalg.matrix_power(R, T)
    print(f"random, radius 1: norm of R^{T} = {np.linalg.norm(J, 2):.3e}")

# Output:
# 0.9^100 = 2.6561398887587544e-05    1.1^100 = 13780.61233982238
# |x| = 4.320329   |Qx| = 4.320329
# norm of Q^1000: 1.0
# random, radius 1: norm of R^10 = 1.803e+00
# random, radius 1: norm of R^100 = 1.467e+00
# random, radius 1: norm of R^1000 = 1.651e+00

La prima riga mostra l'effetto dell'esponente. Dopo mille moltiplicazioni per Q la norma resta 1.0: nessun decadimento e nessuna esplosione, per quanti passi. La matrice casuale a raggio spettrale 1 non è altrettanto pulita, con norme intorno a 1.5-1.8: le matrici non normali possono amplificare temporaneamente anche con raggio spettrale 1. Per questo una pratica comune è inizializzare la matrice ricorrente come ortogonale (in PyTorch esiste torch.nn.init.orthogonal_; qui non l'abbiamo usato).

Un'avvertenza onesta: l'ortogonalità vale per il prodotto di matrici. Il Jacobiano vero contiene anche diag(1 − h²), che è minore di 1, quindi con la tanh il decadimento ritorna (lento, ma ritorna); e durante l'addestramento i pesi si allontanano dall'ortogonalità. È un buon punto di partenza, non una garanzia.

Il rimedio per l'esplosione: gradient clipping

Se il gradiente ha norma enorme, un singolo aggiornamento può lanciare i pesi lontanissimo e rovinare l'addestramento (la perdita diventa nan). Il clipping per norma, proposto da Pascanu et al. (2013), mantiene la direzione ma limita la lunghezza: g ← g · min(1, τ / ‖g‖), dove τ è una soglia da scegliere osservando le norme che compaiono nell'addestramento.

import numpy as np

def clip_by_global_norm(grads, max_norm):
    """Rescale all gradients together if their joint norm exceeds max_norm."""
    total = np.sqrt(sum((g ** 2).sum() for g in grads))
    factor = min(1.0, max_norm / (total + 1e-6))
    return [g * factor for g in grads], total

grads = [np.array([30.0, -40.0]), np.array([0.0, 0.0])]     # exploding: norm = 50
clipped, before = clip_by_global_norm(grads, max_norm=5.0)
print(before, clipped[0])                                    # 50.0, about [3, -4]: same direction, norm 5

# Output:
# 50.0 [ 2.99999994 -3.99999992]

La funzione calcola la norma globale (radice della somma dei quadrati di tutti i gradienti, come se fossero un unico vettore lungo) e riscala tutto con lo stesso fattore min(1, soglia/norma). Due scelte da notare. Il fattore non supera mai 1, quindi i gradienti piccoli non vengono ingranditi. Si usa una norma globale invece di tagliare ogni tensore separatamente, così le proporzioni fra i gradienti dei diversi pesi restano intatte e la direzione di discesa non cambia. Con norma 50 e soglia 5 tutto viene moltiplicato per 0.1: l'uscita è (3, −4) con un'incertezza di 10⁻⁸ dovuta al termine 1e-6 aggiunto per evitare la divisione per zero.

In PyTorch è una riga fra backward() e step(): torch.nn.utils.clip_grad_norm_(params, max_norm), che restituisce la norma prima del clipping (utile da registrare per vedere se la soglia scatta spesso, come si legge nel sorgente). Esiste anche clip_grad_value_, che taglia ogni componente separatamente: fa una cosa diversa e non conserva la direzione del gradiente. Le abbiamo eseguite (PyTorch 2.12, CPU, 2 thread, due volte con la stessa uscita) sugli stessi due gradienti dell'esempio NumPy, norma globale 50.

import torch

torch.set_num_threads(2)

def gradienti():
    """Gli stessi due gradienti dell'esempio NumPy: norma globale 50."""
    p = [torch.zeros(2, requires_grad=True), torch.zeros(2, requires_grad=True)]
    p[0].grad = torch.tensor([30.0, -40.0])
    p[1].grad = torch.tensor([0.0, 0.0])
    return p

p = gradienti()
prima = torch.nn.utils.clip_grad_norm_(p, max_norm=5.0)    # restituisce la norma PRIMA del taglio
g = p[0].grad
print("clip_grad_norm_: norma prima =", prima.item(), " dopo =", round(g.norm().item(), 6))
print("  gradiente:", [round(v, 6) for v in g.tolist()], " direzione:", [round(v, 4) for v in (g / g.norm()).tolist()])

p = gradienti()
ret = torch.nn.utils.clip_grad_value_(p, clip_value=5.0)   # taglia ogni componente a [-5, 5]
g = p[0].grad
print("clip_grad_value_: restituisce", ret, " norma dopo =", round(g.norm().item(), 6))
print("  gradiente:", g.tolist(), " direzione:", [round(v, 4) for v in (g / g.norm()).tolist()])

# Output:
# clip_grad_norm_: norma prima = 50.0  dopo = 5.0
#   gradiente: [3.0, -4.0]  direzione: [0.6, -0.8]
# clip_grad_value_: restituisce None  norma dopo = 7.071068
#   gradiente: [5.0, -5.0]  direzione: [0.7071, -0.7071]

Che cosa mostra l'uscita. clip_grad_norm_ restituisce 50.0, la norma prima del taglio, e porta il gradiente a (3, −4): norma 5 e direzione (0.6, −0.8) invariata, come la nostra funzione (a meno dell'arrotondamento a 32 bit). clip_grad_value_ con soglia 5 porta lo stesso gradiente a (5, −5): la norma diventa 7.07, quindi più della soglia, e la direzione (0.7071, −0.7071) non è più quella di partenza; inoltre restituisce None, quindi non dà la norma da registrare. Gli stessi 30 e −40 producono due vettori diversi perché i due tagli rispondono a domande diverse: uno limita la lunghezza, l'altro limita ogni componente.

Il gradiente che svanisce non si cura con il clipping

Il clipping taglia solo le norme troppo grandi: se il gradiente è già piccolo, scalarlo non lo riporta in vita. Misuriamolo. La rete deve ricordare il primo elemento di una sequenza casuale e restituirlo alla fine, dopo T − 1 passi di rumore (y = x[:, 0]). Predire sempre 0 dà un errore quadratico medio di circa 1, la varianza di una normale standard: «1» significa «non ho imparato niente». Ogni passo genera un batch fresco, quindi la rete non può memorizzare esempi.

Il listato scrive la BPTT a mano e l'ottimizzatore Adam per esteso (è la stessa struttura dell'articolo precedente, con 32 unità). Il punteggio non è la perdita dell'ultimo batch, che è rumorosa, ma l'errore su 2000 sequenze nuove; e l'esperimento è ripetuto su 5 semi, perché con un solo seme si poteva raccontare qualunque storia. Il clipping è disattivato (clip=None) nelle prime tre righe, per isolare il vanishing. L'ultima riga rifà il caso più difficile con clip=1.0.

import numpy as np

def train(T, steps=400, clip=None, seed=0, hidden=32):
    """RNN(1 -> hidden) + linear head, target = FIRST element of the sequence, BPTT by hand, Adam."""
    rng = np.random.default_rng(seed)
    k = 1 / np.sqrt(hidden)
    p = {"Wx": rng.uniform(-k, k, hidden), "Wh": rng.uniform(-k, k, (hidden, hidden)),
         "b": np.zeros(hidden), "w": rng.uniform(-k, k, hidden), "c": np.zeros(())}
    m = {n: np.zeros_like(v) for n, v in p.items()}
    s = {n: np.zeros_like(v) for n, v in p.items()}
    lr, b1, b2, eps = 1e-2, 0.9, 0.999, 1e-8
    for step in range(1, steps + 1):
        x = rng.standard_normal((64, T))
        y = x[:, 0]                                     # target = FIRST element: needs T-1 steps of memory
        hs = [np.zeros((64, hidden))]
        for t in range(T):
            hs.append(np.tanh(x[:, t:t + 1] * p["Wx"] + hs[-1] @ p["Wh"].T + p["b"]))
        err = hs[-1] @ p["w"] + p["c"] - y
        loss = (err ** 2).mean()
        g = {n: np.zeros_like(v) for n, v in p.items()}
        dout = 2 * err / len(err)
        g["w"], g["c"] = hs[-1].T @ dout, dout.sum()
        dh = dout[:, None] * p["w"]
        for t in range(T, 0, -1):                       # BPTT: T steps backwards
            da = dh * (1 - hs[t] ** 2)
            g["Wx"] += (da * x[:, t - 1:t]).sum(axis=0)
            g["Wh"] += da.T @ hs[t - 1]
            g["b"] += da.sum(axis=0)
            dh = da @ p["Wh"]
        if clip:                                        # clip by global norm, as in the block above
            total = np.sqrt(sum((v ** 2).sum() for v in g.values()))
            f = min(1.0, clip / (total + 1e-6))
            g = {n: v * f for n, v in g.items()}
        for n in p:
            m[n] = b1 * m[n] + (1 - b1) * g[n]
            s[n] = b2 * s[n] + (1 - b2) * g[n] ** 2
            p[n] = p[n] - lr * (m[n] / (1 - b1 ** step)) / (np.sqrt(s[n] / (1 - b2 ** step)) + eps)
    # honest score: fresh sequences, not the last training batch
    x = rng.standard_normal((2000, T))
    h = np.zeros((2000, hidden))
    for t in range(T):
        h = np.tanh(x[:, t:t + 1] * p["Wx"] + h @ p["Wh"].T + p["b"])
    return ((h @ p["w"] + p["c"] - x[:, 0]) ** 2).mean()

# Predicting 0 always gives MSE ~ 1.0 (the variance of a standard normal): the "I learned nothing" baseline.
for T in (5, 20, 60):
    scores = [train(T, seed=s) for s in range(5)]
    print(f"T={T:3d}  test MSE over 5 seeds: min {min(scores):.3f}  median {sorted(scores)[2]:.3f}  max {max(scores):.3f}")

# Clipping on the hardest case: does it rescue the long dependency?
scores = [train(60, clip=1.0, seed=s) for s in range(5)]
print(f"T= 60  with clip=1.0:  min {min(scores):.3f}  median {sorted(scores)[2]:.3f}  max {max(scores):.3f}")

# Output:
# T=  5  test MSE over 5 seeds: min 0.002  median 0.003  max 0.005
# T= 20  test MSE over 5 seeds: min 0.026  median 0.280  max 0.973
# T= 60  test MSE over 5 seeds: min 0.072  median 0.230  max 1.014
# T= 60  with clip=1.0:  min 0.186  median 0.805  max 1.027

A distanza 5 la rete impara quasi perfettamente (errore mediano 0.003). A distanza 20 il risultato dipende molto dal seme: da 0.026 a 0.973, cioè da «ha imparato» a «non ha imparato nulla». A distanza 60 con 400 passi l'errore mediano è 0.23, con un seme che arriva a 0.072 e uno che resta a 1.014. Quindi la dipendenza lontana non è impossibile da imparare, ma è molto più difficile e meno affidabile: il segnale che arriva al primo passo è debole e l'esito dipende dall'inizializzazione. L'ultima riga è la lezione sul clipping: con soglia 1.0 l'errore mediano a T=60 sale a 0.805, peggio che senza. Non ha curato il vanishing. Ha fatto di peggio, perché con una soglia bassa scatta spesso e rallenta ogni aggiornamento, come se il learning rate fosse più piccolo. Sono 5 semi, 400 passi e un solo compito: indicano un ordine di grandezza, non una legge.

Gli altri rimedi, e i loro limiti

• BPTT troncata. Si srotola solo per k passi (ad esempio 35) e si interrompe il gradiente prima; lo stato si porta avanti ma senza gradiente (h.detach() in PyTorch). Costa meno memoria, ma per costruzione la rete non può imparare dipendenze più lunghe di k passi.
• Inizializzazione ortogonale o identità di W_hh (Le, Jaitly e Hinton, 2015, per l'identità con ReLU): si parte con un prodotto di Jacobiani che non decade. Aiuta all'inizio, non garantisce nulla durante l'addestramento.
• Celle con porte (LSTM, GRU). Cambiano la struttura in modo che il gradiente abbia un percorso quasi lineare nel tempo. È il rimedio strutturale e il tema degli articoli successivi.

I tre rimedi agiscono a livelli diversi: il primo limita fin dove il gradiente può viaggiare, il secondo migliora il punto di partenza, il terzo cambia l'architettura. Solo l'ultimo risolve il problema alla radice, e infatti è quello che si usa in pratica. Il clipping, a differenza di questi tre, non tocca il vanishing: limita la lunghezza del gradiente che viaggia, non fin dove arriva.

Come scegliere la soglia del clipping

La procedura pratica è: addestra una volta senza clipping per poche centinaia di passi, registra la norma del gradiente prima del taglio a ogni passo, guardane la distribuzione e scegli τ poco sopra il 95-99° percentile. Un esempio ipotetico, con numeri inventati a scopo didattico: se in 10.000 passi la norma sta nel 95% dei casi fra 0.2 e 1.5, con rarissimi picchi a 80, 300 e 2000, una soglia sensata è τ = 2 o τ = 5. Con τ = 5 un picco da 300 viene ridotto di un fattore 60 e i passi normali non vengono toccati; con τ = 0.1 verrebbero ridotti tutti, e staresti addestrando con un learning rate più piccolo di fatto, che è esattamente ciò che l'esperimento con clip=1.0 ha mostrato in piccolo.

Le trappole in cui si cade davvero

La più frequente è una perdita che diventa nan dopo molte epoche tranquille: quasi sempre è un picco di gradiente, e registrare la norma pre-clipping a ogni passo lo rivela (un salto da 1 a 10⁴ è la firma). La seconda è una soglia troppo bassa, come sopra. La terza è pensare che «la perdita scende» significhi «la rete ha memoria lunga»: può ridurre l'errore sfruttando solo gli ultimi elementi, e serve un test che richieda davvero la dipendenza lontana, come quello del primo elemento. L'ultima riguarda la BPTT troncata: se riusi lo stato fra batch senza staccarlo dal grafo (detach()), PyTorch tiene in memoria tutta la storia e dopo qualche batch esaurisci la RAM o ottieni un errore sul secondo backward().

Verificati con quattro conti

Una RNN lineare con autovalore massimo 0.95 dopo 200 passi dà un fattore 0.95^200 ≈ 3.5·10⁻⁵: il gradiente è svanito. Il clipping non risolve il vanishing perché agisce solo se la norma supera la soglia, e un gradiente vicino a zero ha norma già sotto la soglia. La BPTT troncata limita fin dove il gradiente viaggia all'indietro, il clipping limita quanto è grande il gradiente che viaggia: rimedi diversi per problemi diversi. Nel codice della BPTT si usa += per dWhh perché lo stesso peso compare a ogni passo e il gradiente totale è la somma dei contributi di tutti i passi.

Che cosa sai fare adesso

Sai derivare la BPTT, scriverla in NumPy e collaudarla con le differenze finite. Sai spiegare perché il gradiente di una RNN è un prodotto di Jacobiani diag(1 − h²) · W_hh, misurare il decadimento o la crescita della sua norma e distinguere il vanishing, che il clipping non cura, dall'esplosione, che il clipping contiene. Il passo successivo è applicare tutto a dati veri: RNN su una serie temporale; poi la struttura che cambia il problema alla radice, la cella LSTM. Per i prerequisiti, la backpropagation passo per passo e la discesa del gradiente.

Nella stessa serie rimandano a questa lezione anche: LSTM alla prova: il problema dell'addizione, strati multipli e padding, L'attenzione, derivata per intero: query, chiave, valore e il blocco Transformer scritto da zero.