47 - Backpropagation derivata per intero, passaggio per passaggio

Una catena di montaggio produce un oggetto difettoso. Chi ne è responsabile? Non basta dire «l'ultima postazione»: forse ha lavorato bene su un pezzo già rovinato dalla precedente. Per correggere la catena devi risalire, postazione dopo postazione, calcolando per ciascuna quanto ha contribuito al difetto finale. Se sai che la postazione 3 ha contribuito per un certo grado, e che la postazione 3 dipende dalla 2, puoi ricavare il contributo della 2 combinando il contributo della 3 con quanto la 3 dipende dalla 2. È una contabilità delle responsabilità che procede all'indietro. La backpropagation (propagazione all'indietro dell'errore) è esattamente questa contabilità, e il suo linguaggio matematico è uno solo: la regola della catena del calcolo differenziale.

Questa è la lezione cardine del corso. Una rete neurale si addestra con la discesa del gradiente, e la discesa del gradiente ha bisogno del gradiente della funzione di costo rispetto a ogni peso della rete, che può essere un numero di milioni. La backpropagation li calcola tutti con un costo dello stesso ordine di un solo passaggio in avanti. L'algoritmo è stato reso celebre da D. Rumelhart, G. Hinton e R. Williams nel 1986 (la forma di calcolo all'indietro, nota in generale come modalità inversa della differenziazione automatica, era già stata descritta in precedenza, per esempio da S. Linnainmaa nel 1970). Non ha nulla di misterioso: è la regola della catena applicata con ordine. In questa lezione la derivi per intero, la verifichi con un esempio numerico svolto a mano, la implementi, e controlli il risultato in due modi indipendenti.

Prima di cominciare: la regola della catena

Hai già usato la regola della catena nelle lezioni precedenti. La ripassiamo con precisione, perché tutto si regge su di lei. Se y dipende da u, e u dipende da x, cioè y = f(u) e u = g(x), la derivata di y rispetto a x è il prodotto delle derivate lungo la catena:

dy/dx = (dy/du) · (du/dx)

L'intuizione: se u cambia 3 volte più in fretta di x, e y cambia 2 volte più in fretta di u, allora y cambia 6 volte più in fretta di x. Le «velocità di cambiamento» si moltiplicano. Verifica numerica minima: y = (3x)² con u = 3x, y = u². Allora dy/du = 2u = 6x e du/dx = 3, quindi dy/dx = 18x. Derivando direttamente y = 9x² si ottiene 18x: coincide.

Con più variabili entra una novità che userai in continuazione. Se y dipende da due variabili intermedie u₁ e u₂, che dipendono entrambe da x, la derivata totale è la somma dei contributi di tutti i cammini:

dy/dx = (∂y/∂u₁)(du₁/dx) + (∂y/∂u₂)(du₂/dx)

Ogni neurone di uno strato influisce sulla perdita attraverso tutti i neuroni dello strato successivo a cui è collegato: per questo nelle formule che seguono compariranno somme (che diventeranno prodotti di matrici).

Il grafo computazionale: un solo neurone

Un calcolo si può disegnare come un grafo computazionale: ogni nodo è un'operazione elementare, ogni freccia porta un valore all'operazione successiva. Prendiamo il caso più piccolo che contiene tutti gli ingredienti, un neurone con attivazione sigmoide e perdita cross-entropia, con un solo ingresso:

z = w · x + b              (somma pesata)
a = σ(z)                   (attivazione sigmoide)
L = −log(a)  se y = 1      (cross-entropia per l'etichetta y = 1)

Scegliamo numeri concreti: x = 1,5, w = 0,8, b = −0,5, y = 1. Il passaggio in avanti dà z = 0,8·1,5 − 0,5 = 0,7; a = σ(0,7) ≈ 0,6682; L = −log(0,6682) ≈ 0,4032. Questo è l'unico calcolo che serve fare «in avanti»; ogni valore intermedio va ricordato, perché servirà al ritorno.

Il disegno riporta in blu sopra ogni nodo il valore del passaggio in avanti e in rosso sotto ogni nodo la derivata della perdita rispetto a quel valore, calcolata all'indietro. Ripercorriamo i calcoli rossi, da destra a sinistra, applicando la regola della catena un passo alla volta. Si parte dalla fine: la derivata della perdita rispetto a sé stessa vale 1.

Il primo passo, dalla perdita al valore a, usa L = −log a, quindi ∂L/∂a = −1/a = −1/0,6682 ≈ −1,497: se a aumentasse di una piccola quantità, la perdita scenderebbe di 1,497 volte quella quantità. Il secondo passo attraversa la sigmoide: ∂a/∂z = a(1 − a) = 0,6682·0,3318 ≈ 0,2217 (la derivata della sigmoide, dimostrata nella lezione sulla regressione logistica), e per la regola della catena ∂L/∂z = (∂L/∂a)·(∂a/∂z) = −1,497 · 0,2217 ≈ −0,3318. Questo coincide con a − y = 0,6682 − 1 = −0,3318, la semplificazione «uscita meno bersaglio» che avevi derivato per la cross-entropia, ed è il motivo per cui conviene fermarsi su z invece di andare avanti.

Da z il gradiente si distribuisce agli ingredienti della somma pesata, e qui l'ordine non è più rilevante perché i tre calcoli sono indipendenti. Per il bias, z = wx + b dà ∂z/∂b = 1, quindi ∂L/∂b = (∂L/∂z)·1 = −0,3318: la somma «lascia passare» il gradiente invariato. Per il peso, ∂z/∂w = x = 1,5, quindi ∂L/∂w = −0,3318·1,5 ≈ −0,4977: il prodotto manda il gradiente al peso scalato dall'altro fattore, l'ingresso. Per l'ingresso, ∂z/∂x = w = 0,8, quindi ∂L/∂x ≈ −0,2654; nell'unico neurone l'ingresso è un dato e non serve, ma in una rete x è l'uscita dello strato precedente, e questo è il gradiente che si propaga all'indietro.

Due regole pratiche emergono, e sono tutto ciò che serve per qualsiasi grafo: una somma distribuisce il gradiente identico a tutti gli addendi; un prodotto lo distribuisce moltiplicandolo per l'altro fattore. Verifichiamo i numeri con il codice, e controlliamoli in modo indipendente con le differenze finite viste nella lezione sulla discesa del gradiente:

import numpy as np

sigmoide = lambda z: 1 / (1 + np.exp(-z))
x, w, b, y = 1.5, 0.8, -0.5, 1.0

def perdita(w, b, x):                              # cross-entropia con etichetta y = 1
    return -np.log(sigmoide(w * x + b))

# passaggio in avanti, ricordando i valori intermedi
z = w * x + b
a = sigmoide(z)
L = perdita(w, b, x)

# passaggio all'indietro, nodo per nodo
dL_da = -1 / a
da_dz = a * (1 - a)
dL_dz = dL_da * da_dz                              # regola della catena
dL_dw, dL_db, dL_dx = dL_dz * x, dL_dz * 1.0, dL_dz * w

print("avanti:    z = %.4f   a = %.4f   L = %.4f" % (z, a, L))
print("indietro:  dL/da = %.4f   da/dz = %.4f   dL/dz = %.4f  (a - y = %.4f)" % (dL_da, da_dz, dL_dz, a - y))
print("gradienti: dL/dw = %.5f   dL/db = %.5f   dL/dx = %.5f" % (dL_dw, dL_db, dL_dx))

eps = 1e-6
num = ((perdita(w + eps, b, x) - perdita(w - eps, b, x)) / (2 * eps),
       (perdita(w, b + eps, x) - perdita(w, b - eps, x)) / (2 * eps),
       (perdita(w, b, x + eps) - perdita(w, b, x - eps)) / (2 * eps))
print("numerici:  dL/dw = %.5f   dL/db = %.5f   dL/dx = %.5f" % num)

avanti:    z = 0.7000   a = 0.6682   L = 0.4032
indietro:  dL/da = -1.4966   da/dz = 0.2217   dL/dz = -0.3318  (a - y = -0.3318)
gradienti: dL/dw = -0.49772   dL/db = -0.33181   dL/dx = -0.26545
numerici:  dL/dw = -0.49772   dL/db = -0.33181   dL/dx = -0.26545

Cosa fa il codice. Prima calcola il passaggio in avanti salvando z e a (le due quantità intermedie che servono al ritorno). Poi percorre il grafo all'indietro: ogni riga è una moltiplicazione di una derivata «che arriva» da destra per una derivata «locale» del nodo. I tre gradienti finali coincidono con quelli ottenuti perturbando i parametri di una quantità minuscola (differenze finite centrali), fino all'ultima cifra stampata (cinque decimali). È la prima conferma: la regola della catena lavora come promesso.

Una rete a strati: la derivazione generale

Ora la stessa idea su una rete con L strati e un batch di m esempi. Ripetiamo la notazione della lezione precedente, con la prima dimensione sempre il batch:

A⁰ = X                                  (m, n₀)
Z^l = A^{l−1} · W^l + b^l               (m, n_l)       preattivazione dello strato l
A^l = f_l(Z^l)                          (m, n_l)       attivazione (applicata elemento per elemento)
J   = (1/m) Σ_i ℓ(A^L_i, y_i)                          perdita media sul batch

Vogliamo ∂J/∂W^l e ∂J/∂b^l per ogni strato l. La chiave è introdurre una quantità intermedia, il termine d'errore (o «delta») dello strato l: la derivata della perdita rispetto alla preattivazione.

δ^l = ∂J/∂Z^l          forma (m, n_l)   (una riga per esempio, una colonna per neurone)

Scegliere la preattivazione Z (e non l'attivazione A) come variabile di riferimento è comodo perché Z è il punto in cui i pesi «entrano»: una volta che conosciamo δ^l, i gradienti dei pesi dello strato l escono subito. Procediamo in tre passi, esaminando un singolo esempio e un singolo neurone, e poi riscrivendo in forma di matrice.

Passo 1: i gradienti dei pesi, dato δ

La preattivazione del neurone k dello strato l è Z^l_k = Σ_j A^{l−1}_j W^l_{jk} + b^l_k. Dipende dal peso W^l_{jk} solo attraverso il termine A^{l−1}_j · W^l_{jk}, quindi ∂Z^l_k/∂W^l_{jk} = A^{l−1}_j e ∂Z^l_k/∂b^l_k = 1. Per la regola della catena (e poiché W^l_{jk} influisce sulla perdita solo tramite Z^l_k):

∂J/∂W^l_{jk} = δ^l_k · A^{l−1}_j          e          ∂J/∂b^l_k = δ^l_k

In parole: il gradiente di un peso è il prodotto dell'errore del neurone a cui il peso arriva per l'attivazione del neurone da cui il peso parte. Per tutto il batch si sommano (e si mediano) i contributi di tutti gli esempi. Nella forma a matrice:

∂J/∂W^l = (1/m) · (A^{l−1})ᵀ · δ^l              forma:  (n_{l−1}, m) · (m, n_l) = (n_{l−1}, n_l)   ✓ come W^l
∂J/∂b^l = (1/m) · Σ_righe δ^l                    forma:  (n_l,)                                      ✓ come b^l

La trasposta e la somma sulle righe nascono dalla somma su tutti gli esempi del batch. Il controllo delle forme è già una verifica: il gradiente di un parametro deve avere sempre la stessa forma del parametro.

Passo 2: l'errore dello strato finale

Nell'ultimo strato, δ^L = ∂J/∂Z^L si calcola direttamente. Per le tre combinazioni attivazione-perdita standard hai già derivato il risultato nelle lezioni precedenti, ed è sempre lo stesso:

uscita lineare + errore quadratico (½(ŷ − y)²):       δ^L = ŷ − y
uscita sigmoide + cross-entropia binaria:              δ^L = a − y
uscita softmax + cross-entropia categorica:            δ^L = P − Y     (Y = etichette in codifica one-hot)

(Il fattore 1/m della media sul batch si applica una volta sola, nella formula del gradiente dei pesi.) Questo è il punto di partenza del cammino all'indietro: uscita meno bersaglio.

Passo 3: la ricorrenza all'indietro

Abbiamo δ^L. Come si ottiene δ^l partendo da δ^{l+1}? Il neurone j dello strato l influisce sulla perdita attraverso tutti i neuroni k dello strato successivo (somma sui cammini), perché Z^{l+1}_k = Σ_j A^l_j W^{l+1}_{jk} + b^{l+1}_k dipende da A^l_j con coefficiente W^{l+1}_{jk}. Applichiamo la regola della catena in due tratti, prima fino ad A^l_j e poi fino a Z^l_j (dove A^l_j = f_l(Z^l_j)):

∂J/∂A^l_j = Σ_k (∂J/∂Z^{l+1}_k) · (∂Z^{l+1}_k/∂A^l_j) = Σ_k δ^{l+1}_k · W^{l+1}_{jk}

δ^l_j = ∂J/∂Z^l_j = (∂J/∂A^l_j) · f_l'(Z^l_j) = ( Σ_k δ^{l+1}_k · W^{l+1}_{jk} ) · f_l'(Z^l_j)

Nella forma di matrice, la somma su k è un prodotto di δ^{l+1} per la trasposta di W^{l+1}, e il prodotto per f' è elemento per elemento (indicato con ⊙):

δ^l = ( δ^{l+1} · (W^{l+1})ᵀ ) ⊙ f_l'(Z^l)          forma: (m, n_{l+1}) · (n_{l+1}, n_l) = (m, n_l)  ✓

Questa è l'equazione della backpropagation. Leggila come una frase: «l'errore di uno strato è l'errore dello strato dopo, rimandato indietro attraverso i pesi (la trasposta), e filtrato dalla derivata della non linearità». Un neurone con derivata nulla (per esempio una ReLU spenta) blocca del tutto il passaggio dell'errore: da qui il fenomeno delle ReLU «morte». E i pesi fanno da canali di ritorno: lo stesso peso che trasportava il segnale in avanti trasporta la responsabilità all'indietro.

Le quattro equazioni, in sintesi

(1) δ^L     = A^L − Y                                       (uscita meno bersaglio, per le tre coppie standard)
(2) δ^l     = ( δ^{l+1} · (W^{l+1})ᵀ ) ⊙ f_l'(Z^l)         (per l = L−1, ..., 1)
(3) ∂J/∂W^l = (1/m) · (A^{l−1})ᵀ · δ^l
(4) ∂J/∂b^l = (1/m) · somma sulle righe di δ^l

Tutto qui. L'algoritmo è: (a) passaggio in avanti salvando tutte le Z^l e le A^l; (b) calcolare δ^L con la (1); (c) per ogni strato, dall'ultimo al primo, calcolare i gradienti con la (3) e (4) e poi, se l > 1, propagare δ all'indietro con la (2). Il costo del passaggio all'indietro è dello stesso ordine di quello in avanti, perché ogni strato richiede prodotti di matrici delle stesse dimensioni.

Un esempio svolto a mano: la rete 2 → 2 → 1

Riprendiamo la rete della lezione precedente: 2 ingressi, 2 neuroni nascosti con ReLU, 1 uscita lineare, con x = (1, 2), W¹ = [[0,5; −1,0], [0,25; 0,5]], b¹ = (0, 1), W² = [[2,0], [−1,0]], b² = 0,5. Avevamo calcolato Z¹ = (1, 1), A¹ = (1, 1) e uscita ŷ = 1,5. Aggiungiamo il bersaglio y = 1 e la perdita quadratica J = ½(ŷ − y)² = ½·0,25 = 0,125.

Applichiamo le quattro equazioni, una alla volta (i valori sono quelli scritti in rosso nel disegno):

L'ordine dei passi conta, perché ognuno usa il risultato del precedente. L'errore d'uscita (equazione 1) è δ² = ŷ − y = 1,5 − 1 = 0,5, l'unica quantità che parte dalla perdita. Con le equazioni (3-4) i gradienti dello strato 2 sono ∂J/∂W² = (A¹)ᵀ δ² = (1, 1)ᵀ · 0,5 = (0,5; 0,5)ᵀ, di forma (2, 1) come W², e ∂J/∂b² = 0,5.

Per propagare l'errore (equazione 2) si calcola δ² · (W²)ᵀ = 0,5 · (2,0; −1,0) = (1,0; −0,5): il neurone nascosto 1 riceve il 100% del segnale d'errore moltiplicato per il peso +2, il neurone 2 riceve il −50% a causa del peso −1. Moltiplicando per la derivata della ReLU in Z¹ = (1, 1), che vale 1 in entrambi i casi (Z positivo), si ha δ¹ = (1,0; −0,5). Con (3-4) sullo strato 1: ∂J/∂W¹ = xᵀ δ¹ = (1; 2)ᵀ · (1,0; −0,5) = [[1,0; −0,5], [2,0; −1,0]] e ∂J/∂b¹ = (1,0; −0,5).

Un'interpretazione leggibile: il gradiente del peso che collega l'ingresso 2 al neurone nascosto 1 vale 2,0, il doppio di quello dell'ingresso 1 verso lo stesso neurone (1,0), perché l'ingresso 2 vale il doppio. Il segno negativo dei gradienti verso il neurone nascosto 2 dice che per ridurre la perdita quel neurone dovrebbe aumentare la sua attivazione (la rete sta sottostimando in quel ramo). Verifichiamo tutto con il codice, ancora una volta con le differenze finite come arbitro indipendente:

x = np.array([[1.0, 2.0]]); y = np.array([[1.0]])
W1 = np.array([[0.5, -1.0], [0.25, 0.5]]); b1 = np.array([0.0, 1.0])
W2 = np.array([[2.0], [-1.0]]);            b2 = np.array([0.5])

def passaggio_avanti(W1, b1, W2, b2):
    Z1 = x @ W1 + b1; A1 = np.maximum(0, Z1); Z2 = A1 @ W2 + b2
    return Z1, A1, Z2, 0.5 * ((Z2 - y) ** 2).mean()

Z1, A1, Z2, J = passaggio_avanti(W1, b1, W2, b2)

# passaggio all'indietro, con le equazioni (1)-(4)
d2 = Z2 - y                                        # (1) uscita meno bersaglio:               (1, 1)
dW2 = A1.T @ d2                                    # (3) (A1)^T d2:                          (2, 1)
db2 = d2.sum(axis=0)                               # (4)
d1 = (d2 @ W2.T) * (Z1 > 0)                        # (2) d2 W2^T, poi * derivata ReLU:       (1, 2)
dW1 = x.T @ d1                                     # (3) x^T d1:                             (2, 2)
db1 = d1.sum(axis=0)                               # (4)

print("J = %.3f   delta2 = %s   delta1 = %s" % (J, d2.ravel(), d1.ravel()))
print("dW2 =", dW2.ravel(), "  db2 =", db2, "\ndW1 =", dW1.tolist(), "  db1 =", db1)

def gradiente_numerico(indice, eps=1e-6):          # perturba un parametro alla volta
    params = [W1.copy(), b1.copy(), W2.copy(), b2.copy()]
    out = np.zeros_like(params[indice])
    for i in np.ndindex(out.shape):
        piu = [p.copy() for p in params]; meno = [p.copy() for p in params]
        piu[indice][i] += eps; meno[indice][i] -= eps
        out[i] = (passaggio_avanti(*piu)[3] - passaggio_avanti(*meno)[3]) / (2 * eps)
    return out

for nome, an, k in (("W1", dW1, 0), ("b1", db1, 1), ("W2", dW2, 2), ("b2", db2, 3)):
    print("%s: massima differenza tra analitico e numerico = %.2e" % (nome, np.abs(an - gradiente_numerico(k)).max()))

J = 0.125   delta2 = [0.5]   delta1 = [ 1.  -0.5]
dW2 = [0.5 0.5]   db2 = [0.5] 
dW1 = [[1.0, -0.5], [2.0, -1.0]]   db1 = [ 1.  -0.5]
W1: massima differenza tra analitico e numerico = 6.04e-11
b1: massima differenza tra analitico e numerico = 1.98e-11
W2: massima differenza tra analitico e numerico = 7.44e-12
b2: massima differenza tra analitico e numerico = 7.44e-12

Cosa fa il codice. Le righe con i commenti (1)-(4) sono le quattro equazioni, trascritte quasi alla lettera. Nota (Z1 > 0): è la derivata della ReLU (vero dove Z è positivo, falso altrove), che NumPy converte in 1 e 0 quando la moltiplica. La funzione gradiente_numerico perturba un parametro alla volta in avanti e indietro, richiamando il passaggio in avanti: è lento (due passaggi in avanti per ogni parametro) ma indipendente dalla teoria, e per questo è l'arbitro ideale. Le differenze massime sono dell'ordine di 10⁻¹¹ o meno: i gradienti calcolati a mano coincidono con quelli numerici.

Implementazione generale e gradient check su una rete 4 → 8 → 3

Mettiamo le quattro equazioni in una classe con un numero qualsiasi di strati. Per l'ultimo strato usiamo softmax e cross-entropia (equazione 1: δ = P − Y), per gli strati nascosti ReLU o sigmoide. Dopo il passaggio in avanti, la classe conserva tutte le Z e le A: sono indispensabili al ritorno, ed è il motivo per cui l'addestramento di una rete consuma molta più memoria dell'uso in inferenza.

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

Xi, yi = load_iris(return_X_y=True)
Xi = StandardScaler().fit_transform(Xi)

def softmax(Z):
    Z = Z - Z.max(axis=1, keepdims=True)
    E = np.exp(Z)
    return E / E.sum(axis=1, keepdims=True)

class Rete:
    def __init__(self, dimensioni, attivazione="relu", seed=0):
        rng = np.random.default_rng(seed)
        self.attivazione = attivazione
        self.W = [rng.normal(size=(i, o)) * np.sqrt(2.0 / i) for i, o in zip(dimensioni[:-1], dimensioni[1:])]
        self.b = [np.zeros(o) for o in dimensioni[1:]]

    def f(self, Z):
        return np.maximum(0, Z) if self.attivazione == "relu" else 1 / (1 + np.exp(-Z))

    def df(self, Z, A):                                    # derivata di f, in termini di Z (ReLU) o di A (sigmoide: A(1-A))
        return (Z > 0).astype(float) if self.attivazione == "relu" else A * (1 - A)

    def avanti(self, X):
        self.Z, self.A = [], [X]                           # si conservano tutte le preattivazioni e attivazioni
        A = X
        for l, (W, b) in enumerate(zip(self.W, self.b)):
            Z = A @ W + b
            A = softmax(Z) if l == len(self.W) - 1 else self.f(Z)
            self.Z.append(Z); self.A.append(A)
        return A

    def perdita(self, X, y):
        P = self.avanti(X)
        return -np.log(P[np.arange(len(y)), y]).mean()

    def indietro(self, y):
        m = len(y)
        delta = self.A[-1] - np.eye(self.A[-1].shape[1])[y]            # (1): uscita meno bersaglio (softmax + cross-entropia)
        gW, gb = [None] * len(self.W), [None] * len(self.W)
        for l in range(len(self.W) - 1, -1, -1):
            gW[l] = self.A[l].T @ delta / m                            # (3)
            gb[l] = delta.sum(axis=0) / m                              # (4)
            if l > 0:                                                  # (2): propagare all'indietro, tranne che oltre il primo strato
                delta = (delta @ self.W[l].T) * self.df(self.Z[l - 1], self.A[l])
        return gW, gb

rete = Rete([4, 8, 3])
rete.perdita(Xi, yi)
gW, gb = rete.indietro(yi)

def gradiente_numerico_rete(rete, X, y, eps=1e-5):
    risultato = []
    for P in rete.W + rete.b:
        g = np.zeros_like(P)
        for idx in np.ndindex(P.shape):
            vecchio = P[idx]
            P[idx] = vecchio + eps; piu = rete.perdita(X, y)
            P[idx] = vecchio - eps; meno = rete.perdita(X, y)
            P[idx] = vecchio; g[idx] = (piu - meno) / (2 * eps)
        risultato.append(g)
    return risultato

num = gradiente_numerico_rete(rete, Xi, yi)
print("errore relativo analitico/numerico (norma della differenza / norma della somma):")
for nome, a, n in zip(["W1 (4x8)", "W2 (8x3)", "b1 (8,)", "b2 (3,)"], gW + gb, num):
    print("  %-9s %.2e" % (nome, np.linalg.norm(a - n) / np.linalg.norm(a + n)))
print("parametri: %d  ->  il gradiente numerico ha richiesto %d passaggi in avanti, quello analitico un passaggio in avanti più uno indietro" % (
    sum(W.size + b.size for W, b in zip(rete.W, rete.b)), 2 * sum(W.size + b.size for W, b in zip(rete.W, rete.b))))

errore relativo analitico/numerico (norma della differenza / norma della somma):
  W1 (4x8)  5.40e-11
  W2 (8x3)  2.86e-11
  b1 (8,)   4.79e-11
  b2 (3,)   7.45e-12
parametri: 67  ->  il gradiente numerico ha richiesto 134 passaggi in avanti, quello analitico un passaggio in avanti più uno indietro

Cosa fa il codice. avanti conserva in self.Z e self.A tutti i valori intermedi. indietro è un ciclo che va dall'ultimo strato al primo e per ciascuno applica (3) e (4); poi, se non siamo al primo strato, propaga δ all'indietro con (2): delta @ self.W[l].T è il prodotto con la trasposta, e self.df(self.Z[l - 1], self.A[l]) è la derivata dell'attivazione dello strato precedente. L'indice l - 1 per le Z (e l per le A, che hanno un elemento in più, l'ingresso) è la parte che di solito si sbaglia: l'allineamento degli indici merita un controllo, ed è proprio quello che fa il gradient check.

L'errore relativo è dell'ordine di 10⁻¹⁰ per tutti e quattro i gruppi di parametri: la backpropagation implementata è corretta. L'ultima riga contrappone i costi: per i 67 parametri della rete il gradiente numerico ha richiesto 134 passaggi in avanti; quello analitico, uno in avanti e uno all'indietro. Con un milione di parametri sarebbero due milioni contro due: è questa la ragione per cui la backpropagation ha reso praticabile l'addestramento delle reti.

Un secondo arbitro: la differenziazione automatica di PyTorch

Il gradient check con le differenze finite è un arbitro indipendente dalla teoria. Un secondo arbitro è la differenziazione automatica di una libreria, che costruisce il grafo mentre calcola in avanti e applica la regola della catena all'indietro. Il listato copia gli stessi pesi in tensori float64 con requires_grad=True (float64 perché la nostra rete lavora in doppia precisione; la libreria userebbe float32 di suo), rifà il passaggio in avanti con la perdita di cross-entropia della libreria (che contiene la softmax) e chiama backward(). Riusa rete, Xi, yi, gW e gb del blocco precedente. È stato eseguito con PyTorch 2.12.0 su CPU (CUDA_VISIBLE_DEVICES="", 2 thread), due volte, con uscita identica; le differenze, essendo arrotondamenti, dipendono dalla macchina e dalla versione, e non hanno un valore da ricordare.

import torch

tW = [torch.tensor(W, requires_grad=True) for W in rete.W]       # stessi pesi, in float64
tb = [torch.tensor(b, requires_grad=True) for b in rete.b]
Xt, yt = torch.tensor(Xi), torch.tensor(yi)

Z1 = torch.relu(Xt @ tW[0] + tb[0])
Z2 = Z1 @ tW[1] + tb[1]
J = torch.nn.functional.cross_entropy(Z2, yt)                    # applica softmax + cross-entropia internamente
J.backward()                                                     # backpropagation automatica

print("perdita PyTorch %.6f   perdita a mano %.6f" % (J.item(), rete.perdita(Xi, yi)))
for nome, mio, p in zip(["W1", "W2", "b1", "b2"], gW + gb, tW + tb):
    print("  %s  massima differenza tra i miei gradienti e quelli di PyTorch: %.2e" % (nome, np.abs(mio - p.grad.numpy()).max()))

perdita PyTorch 1.572833   perdita a mano 1.572833
  W1  massima differenza tra i miei gradienti e quelli di PyTorch: 9.71e-17
  W2  massima differenza tra i miei gradienti e quelli di PyTorch: 1.67e-16
  b1  massima differenza tra i miei gradienti e quelli di PyTorch: 5.55e-17
  b2  massima differenza tra i miei gradienti e quelli di PyTorch: 1.04e-16

Cosa fa il codice. Copia in tensori PyTorch gli stessi pesi, rifà il passaggio in avanti con operazioni di PyTorch e chiama J.backward(), che percorre il grafo computazionale costruito durante il calcolo in avanti applicando la regola della catena nodo per nodo. Le perdite coincidono (1,572833) e i gradienti coincidono con i nostri fino all'arrotondamento, con differenze massime nell'ordine di 10⁻¹⁶. Sono due misure diverse (differenza massima in valore assoluto qui, errore relativo sulle norme nel gradient check), ma dicono la stessa cosa: la backpropagation scritta a mano è corretta.

Quello che fa backward() non è diverso da ciò che hai scritto: è lo stesso algoritmo applicato automaticamente a qualsiasi grafo di operazioni differenziabili. Il vantaggio è che non devi derivare a mano ogni nuova architettura; il rischio è non capire che cosa succede quando qualcosa va storto.

Perché il gradiente svanisce nelle reti profonde: lo vedi nella ricorrenza

L'equazione (2) mostra il meccanismo del gradiente che svanisce (e di quello che esplode). Per arrivare dall'ultimo strato al primo, δ viene moltiplicato per una matrice (W)ᵀ e per una derivata f' a ogni strato: in sostanza un prodotto di L−1 fattori. Se i fattori tipici sono minori di 1, il prodotto cala in modo esponenziale; se sono maggiori di 1, cresce. La derivata della sigmoide non supera mai 0,25 (lezione precedente), quindi un tipico fattore per strato è minore di 1 e i gradienti dei primi strati sono minuscoli. Misuriamolo in una rete con 10 strati, confrontando sigmoide e ReLU (con He):

rng = np.random.default_rng(2)
Xr = rng.normal(size=(200, 20))
yr = rng.integers(0, 3, size=200)

for attivazione in ("sigmoid", "relu"):
    profonda = Rete([20] + [30] * 9 + [3], attivazione=attivazione, seed=1)
    if attivazione == "sigmoid":                                # per la sigmoide uso la scala di Xavier, come in letteratura
        profonda.W = [np.random.default_rng(i).normal(size=W.shape) * np.sqrt(1.0 / W.shape[0]) for i, W in enumerate(profonda.W)]
    profonda.perdita(Xr, yr)
    gW_prof, _ = profonda.indietro(yr)
    print("%-8s norma del gradiente dei pesi, dal primo all'ultimo strato:" % attivazione)
    print("         ", ["%.1e" % np.linalg.norm(g) for g in gW_prof])

sigmoid  norma del gradiente dei pesi, dal primo all'ultimo strato:
          ['4.0e-07', '1.4e-06', '6.4e-06', '2.5e-05', '1.1e-04', '5.3e-04', '1.9e-03', '1.0e-02', '4.3e-02', '2.5e-01']
relu     norma del gradiente dei pesi, dal primo all'ultimo strato:
          ['2.5e-01', '4.5e-01', '5.0e-01', '6.3e-01', '8.4e-01', '1.0e+00', '1.8e+00', '2.1e+00', '1.4e+00', '1.6e+00']

Cosa fa il codice. Crea una rete con 10 strati di pesi (20 → 30 → ... → 30 → 3) e calcola il gradiente su 200 esempi casuali, stampando la norma del gradiente dei pesi di ogni strato, dal primo (vicino all'ingresso) all'ultimo. Con la sigmoide la norma cresce di un fattore tra 3,5 e 6 a ogni strato verso l'uscita: il primo strato ha un gradiente di 4·10⁻⁷ contro 0,25 dell'ultimo, un rapporto di circa seicentomila a uno.

Gli strati iniziali praticamente non imparano. Con la ReLU e l'inizializzazione di He tutte le norme stanno tra 0,25 e 2,1: il segnale d'errore arriva ovunque. Quello che qui hai misurato su un esempio è, in piccolo, la ragione per cui la ReLU e la corretta inizializzazione hanno reso addestrabili le reti profonde.

Trade-off e trappole reali

• Non salvare le attivazioni: il passaggio all'indietro ne ha bisogno. L'addestramento consuma memoria proporzionale a profondità × batch × larghezza; è il motivo per cui un modello che «entra» in memoria in inferenza può non entrare in addestramento. Tecniche come il gradient checkpointing scambiano calcolo con memoria.
• Sbagliare l'allineamento degli indici (Z^l contro A^l contro A^{l−1}): è l'errore tipico di un'implementazione a mano. Il gradient check lo cattura sempre; non scrivere una backpropagation senza averlo eseguito almeno una volta.
• Usare la derivata della sigmoide in funzione di Z invece che di A: la formula A(1 − A) usa l'attivazione, non la preattivazione. Applicare σ(1 − σ) a un valore che è già un'attivazione dà un risultato sbagliato senza segnalarlo.
• Dimenticare la divisione per m (la media sul batch): i gradienti risultano m volte più grandi, il learning rate effettivo cambia con la dimensione del batch, e il gradient check fallisce.
• Ridurre male le dimensioni del bias: il gradiente di b si ottiene sommando δ sulle righe (asse 0); sommare sull'asse sbagliato dà forme errate o, peggio, forme «giuste» per caso quando n_l = m.
• Gradient check con la ReLU: vicino a Z = 0 la derivata non esiste, e le differenze finite possono discostarsi. Se un controllo fallisce di poco su un singolo elemento, prova con un altro seme prima di cercare un bug.
• Il gradiente non è l'aggiornamento: la backpropagation calcola ∂J/∂θ. È l'ottimizzatore (SGD, Adam) a decidere come usarlo per aggiornare i pesi.
• Gradienti che esplodono: il caso opposto del gradiente che svanisce. Si limita con il gradient clipping (riscalare il gradiente se la sua norma supera una soglia) e con un'inizializzazione corretta.
• Le librerie accumulano i gradienti: in PyTorch, senza optimizer.zero_grad() a ogni passo, backward() li somma ai precedenti (è voluto, per sommare su più batch, ma è una fonte di bug quando si dimentica; comportamento documentato della libreria, non provocato qui).

Le nove voci si raggruppano in tre famiglie. Alcune sono errori di implementazione che il gradient check intercetta sempre (indici, derivata della sigmoide, divisione per m, asse del bias). Altre sono costi e limiti dell'algoritmo (memoria per le attivazioni, gradienti che svaniscono o esplodono). Le ultime sono confusioni concettuali: il gradiente non è l'aggiornamento, e il gradient check con la ReLU può fallire di poco vicino a zero senza che ci sia un bug.

Domande frequenti

La backpropagation è un algoritmo di apprendimento? No, è un algoritmo per calcolare gradienti. L'apprendimento è la discesa del gradiente (o una sua variante) che usa quei gradienti. Si confondono spesso perché si presentano insieme.

Perché si va all'indietro e non in avanti? Ci sarebbe anche un modo «in avanti» di applicare la regola della catena (differenziazione automatica in modalità diretta), ma richiede un passaggio per ogni parametro. In una rete c'è una sola perdita scalare e moltissimi parametri: la modalità inversa calcola la derivata di un output rispetto a tutti gli input in un solo passaggio. È un caso di scelta dell'ordine di moltiplicazione che rende il calcolo economico.

Il cervello usa la backpropagation? È una questione di ricerca aperta, e non è necessario risolverla per usare l'algoritmo. Non c'è nessuna evidenza semplice che il cervello propaghi gradienti esatti all'indietro attraverso le sinapsi nel modo descritto qui; la backpropagation va considerata un metodo di ingegneria e non un modello del cervello.

Devo implementarla a mano ogni volta? No. Nel lavoro quotidiano si usa la differenziazione automatica di PyTorch o TensorFlow. Saperla scrivere ti serve per capire che cosa fa la libreria, per diagnosticare i problemi (gradienti che svaniscono, ReLU morte, forme sbagliate) e per leggere con spirito critico il codice altrui.

Che cosa sai fare adesso

Sai derivare la backpropagation con la regola della catena, applicarla a mano a un neurone e a una rete 2 → 2 → 1, e implementarla per una rete con un numero qualsiasi di strati controllando la forma di ogni gradiente. Sai verificare il risultato con le differenze finite e riconoscere dalla ricorrenza perché il gradiente svanisce con la sigmoide: nell'esperimento a dieci strati la norma del primo strato era circa seicentomila volte più piccola di quella dell'ultimo. Le lezioni di riferimento sono le reti a strati e il forward pass e la discesa del gradiente; nella prossima, una rete neurale da zero in NumPy, metti insieme il tutto su dati veri.

Nella stessa serie rimandano a questa lezione anche: Analisi matematica 1: limiti e derivate, Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero, Il percettrone e le funzioni di attivazione: a cosa serve la non linearità, La convoluzione: filtri, kernel e perché guardare localmente, Reti neurali ricorrenti (RNN): dare una memoria alla rete, Backpropagation through time: gradiente che svanisce e che esplode, L'attenzione, derivata per intero: query, chiave, valore e il blocco Transformer scritto da zero.