25 - Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero

Corso Machine Learning Engineer · Modulo 02 · Matematica e statistica · Lezione 10

Obiettivo. Derivare una funzione rispetto a un vettore e rispetto a una matrice senza sbagliare le forme, e capire come fa una libreria a calcolare da sola i gradienti di milioni di parametri. Alla fine sai dichiarare una convenzione di layout e rispettarla, sai ricavare da solo le cinque identità che ricorrono in ogni rete, sai distinguere il modo diretto dal modo inverso in base al costo, e sai scrivere un autograd di centoventi righe i cui gradienti coincidono con le differenze finite centrate entro un errore relativo di 10⁻⁷.

Falso mito da sfatare. «La derivata rispetto a una matrice è una cosa esotica, roba da matematici.» No: è la stessa derivata parziale di sempre, ripetuta per ogni casella e rimessa a posto in una tabella della stessa forma della matrice. L'unica difficoltà vera è decidere come impilare i risultati, e il 90% degli errori di segno e di forma nasce dal non aver dichiarato quella decisione.

Introduzione

Nella lezione sul gradiente e la discesa del gradiente hai derivato la catena su un neurone solo: un numero dentro, un numero fuori, una derivata. Nella lezione su vettori, matrici e sistemi hai imparato a moltiplicare matrici e a leggere le forme. Queste due cose, messe insieme, non bastano ancora: nella lezione sulla backpropagation ti servirà ∂L/∂W con W di forma (128, 64), cioè la derivata di un numero solo rispetto a 8192 numeri contemporaneamente. Questa lezione è il ponte.

È anche la lezione che chiude un debito: finora abbiamo usato ŷ − y senza derivarlo e abbiamo chiamato «autograd» una scatola nera. Alla fine non c'è più niente di nero dentro quella scatola.

Risultato finale di questa lezione. Un file di centoventi righe che costruisce un grafo di calcolo, lo percorre all'indietro una sola volta e produce i gradienti di una rete 2→16→16→2, verificati contro due oracoli indipendenti: le formule chiuse derivate a mano e le differenze finite centrate.

Roadmap della lezione.

1. Derivare rispetto a un vettore, e perché la convenzione di layout va dichiarata.
2. Le cinque identità che servono davvero, derivate per intero, componente per componente.
3. Il Jacobiano, il prodotto vettore-Jacobiano e il conto della memoria.
4. L'Hessiano, la curvatura, il passo di Newton e perché nel deep learning non si usa.
5. Modo diretto e modo inverso: lo stesso risultato, due costi diversi.
6. Un autograd a modo inverso scritto da zero.
7. Il gradient check fatto bene: differenze centrate, scelta di h, errore relativo.
8. L'esperimento: una spirale a due bracci, addestrata con il nostro autograd.

Il codice di questa lezione è stato eseguito con Python 3.14.4 e NumPy 2.3.5, e usa soltanto NumPy: nessuna libreria di deep learning, perché il punto è proprio costruirne il pezzo centrale. Dove compaiono dati casuali il generatore è inizializzato con un seme fisso, quindi gli output sono riproducibili riga per riga; con versioni diverse di NumPy possono cambiare le ultime cifre, non le conclusioni.

Derivare rispetto a un vettore

L'analogia: il cruscotto con molte manopole

Immagina una centralina con 8192 manopole e un unico quadrante che mostra la perdita. La derivata parziale rispetto a una manopola è di quanto si muove l'ago per millimetro di quella manopola, tenendo ferme tutte le altre.

Il gradiente non è una di queste letture: è la lettura simultanea di tutte e 8192, annotata accanto a ogni manopola. E il punto è che quel foglio ha la forma del cruscotto: se le manopole stanno in una griglia 128 per 64, anche le annotazioni ci stanno. Non è estetica: è ciò che permette di scrivere l'aggiornamento come W ← W − η · ∂L/∂W, una sottrazione fra due tabelle della stessa forma.

Definizione. Data una funzione a valori reali f(x) con x ∈ Rⁿ, il gradiente è il vettore delle derivate parziali

∇f(x) = ∂f/∂x = (∂f/∂x₁, ∂f/∂x₂, …, ∂f/∂xₙ).

Leggila come una frase: la componente i-esima dice quanto cresce f se muovi di poco la sola coordinata i, lasciando ferme le altre.

La convenzione di layout, dichiarata una volta per tutte

Qui comincia la parte che i libri danno per scontata e che costa ore di debug. Quando derivi uno scalare rispetto a un vettore ottieni n numeri, e devi decidere se metterli in colonna o in riga. Le due scelte hanno un nome.

In questa lezione, e in tutto il corso, adottiamo la convenzione del denominatore per i gradienti: ∂f/∂x ha sempre la forma di x. La ragione è pratica: l'aggiornamento è una sottrazione fra il parametro e il suo gradiente, e se le due forme non coincidono devi trasporre a ogni passo — e prima o poi dimentichi.

Per i Jacobiani usiamo invece la convenzione del numeratore: J ha forma (m, n), una riga per uscita e una colonna per ingresso, la forma in cui la regola della catena si scrive come prodotto di matrici nell'ordine naturale. Le due convenzioni convivono perché riguardano oggetti diversi: il gradiente di uno scalare e il Jacobiano di una funzione vettoriale. Dichiararlo in testa al file vale più di qualunque commento successivo.

Il gradiente trasposto che quasi funziona

Ecco l'errore più costoso della lezione, e vale la pena vederlo prima delle formule. Calcoli il gradiente rispetto a una matrice di pesi e lo scrivi trasposto. Se la matrice è rettangolare, NumPy solleva un errore di forma alla prima sottrazione e in trenta secondi hai trovato il bug. Se è quadrata — e negli strati nascosti lo è quasi sempre, 512→512 — nessuno protesta: il codice gira, la perdita magari scende anche un po', e tu passi due giorni a chiederti perché la rete impara male.

La regola da portare via. Un gradiente trasposto su una matrice quadrata non produce nessun errore: produce un numero sbagliato. È per questo che il gradient check non è un lusso da perfezionisti ma la prima cosa da scrivere, prima ancora del ciclo di addestramento.

Le cinque identità che servono davvero

Esistono manuali con centinaia di identità di calcolo matriciale — The Matrix Cookbook ne raccoglie moltissime. Nella pratica del machine learning ne servono cinque, e le deriviamo tutte e cinque componente per componente, senza mai scrivere «si dimostra che». Il metodo è sempre lo stesso: scrivi la funzione con gli indici e le somme esplicite, deriva rispetto alla componente k-esima come faresti con una funzione di una variabile, poi riconosci nel risultato un prodotto fra matrici.

Le prime due: il prodotto scalare e la forma quadratica

Identità 1. f(x) = aᵀx = Σᵢ aᵢ xᵢ. Deriviamo rispetto a x_k: nella somma c'è un solo addendo che contiene x_k, ed è a_k x_k; tutti gli altri sono costanti e spariscono. Quindi ∂f/∂x_k = a_k. Impilando le n componenti:

∂(aᵀx)/∂x = a.

In parole: la derivata di una funzione lineare è il suo vettore di coefficienti, esattamente come la derivata di ax è a nel caso scalare. ∎

Identità 2. f(x) = xᵀAx = Σᵢ Σⱼ xᵢ Aᵢⱼ xⱼ. Qui x_k compare in due famiglie di addendi: quelli con i = k (cioè Σⱼ x_k A_kⱼ xⱼ) e quelli con j = k (cioè Σᵢ xᵢ Aᵢ_k x_k), più il termine A_kk x_k² che appartiene a entrambe le famiglie. Deriviamo con la regola del prodotto, trattando separatamente i due ruoli:

∂f/∂x_k = Σⱼ A_kⱼ xⱼ + Σᵢ xᵢ Aᵢ_k = (Ax)_k + (Aᵀx)_k

Il doppio conteggio del termine diagonale si sistema da sé: ∂(A_kk x_k²)/∂x_k = 2A_kk x_k, che è proprio il contributo che le due somme danno insieme. Impilando:

∂(xᵀAx)/∂x = (A + Aᵀ)x, che diventa 2Ax solo se A è simmetrica. ∎

Quel «solo se» è la seconda trappola della lezione. La formula 2Ax è così comune da sembrare universale, ma vale per le matrici simmetriche — covarianze, AᵀA, matrici di Gram. Su una matrice qualunque sbaglia, e sbaglia in modo silenzioso, perché la forma del risultato è giusta.

import numpy as np

def riga(v, cifre=4):
    """Stampa un array a cifre fisse: due esecuzioni danno lo stesso testo."""
    return " ".join(f"{val:+.{cifre}f}" for val in np.ravel(v))

def gradiente_numerico(f, x, h=1e-5):
    """Differenze centrate, una componente alla volta. Funziona anche sulle matrici."""
    g = np.zeros_like(x, dtype=float)
    for i in range(x.size):
        passo = np.zeros_like(x, dtype=float)
        passo.flat[i] = h
        g.flat[i] = (f(x + passo) - f(x - passo)) / (2 * h)
    return g

a = np.array([1.0, 2.0, 3.0])
x = np.array([4.0, 5.0, 6.0])
print("f(x) = aᵀx")
print("  analitico a        :", riga(a))
print("  numerico           :", riga(gradiente_numerico(lambda v: a @ v, x)))

A = np.array([[2.0, 1.0],
              [0.0, 3.0]])                       # NON simmetrica, di proposito
x2 = np.array([1.0, 2.0])
print("f(x) = xᵀAx, con A non simmetrica")
print("  (A + Aᵀ)x          :", riga((A + A.T) @ x2))
print("  2Ax (qui sbagliato):", riga(2 * (A @ x2)))
print("  numerico           :", riga(gradiente_numerico(lambda v: v @ A @ v, x2)))

Output ottenuto:

f(x) = aᵀx
  analitico a        : +1.0000 +2.0000 +3.0000
  numerico           : +1.0000 +2.0000 +3.0000
f(x) = xᵀAx, con A non simmetrica
  (A + Aᵀ)x          : +6.0000 +13.0000
  2Ax (qui sbagliato): +8.0000 +12.0000
  numerico           : +6.0000 +13.0000

Cosa fa il codice. gradiente_numerico perturba una componente alla volta di ±h e calcola il rapporto incrementale centrato; usa passo.flat[i] così la stessa funzione vale per vettori e per matrici, mentre riga formatta a cifre fisse perché la stampa di NumPy cambia allineamento a seconda dei valori. Per aᵀx il numerico coincide con a. Per la forma quadratica con A non simmetrica dà (6, 13), che è (A + Aᵀ)x, non 2Ax = (8, 12): la scorciatoia sbaglia su entrambe le componenti, con la forma giusta. Le differenze centrate fanno qui da giudice imparziale, e sono esatte sui polinomi di secondo grado — proprietà che dimostreremo parlando di gradient check.

Minimi quadrati e prodotto esterno

Identità 3. f(x) = ‖Ax − b‖². Scriviamo il residuo con gli indici: rᵢ = Σⱼ Aᵢⱼ xⱼ − bᵢ, e f = Σᵢ rᵢ². Per la regola della catena, ∂f/∂x_k = Σᵢ 2 rᵢ · ∂rᵢ/∂x_k, e ∂rᵢ/∂x_k = Aᵢ_k perché bᵢ è costante. Dunque ∂f/∂x_k = 2 Σᵢ Aᵢ_k rᵢ, che è la k-esima componente di Aᵀr:

∂‖Ax − b‖²/∂x = 2Aᵀ(Ax − b). ∎

Questa è la derivata dei minimi quadrati visti nella lezione di statistica. Annullarla dà AᵀAx = Aᵀb, cioè l'equazione normale della lezione su vettori e matrici: il punto di minimo è la soluzione di un sistema lineare, e tutta la regressione lineare sta in queste due righe.

Identità 4. Ora deriviamo rispetto a una matrice. Sia z = Wx con W di forma (m, n), e sia L una perdita scalare che dipende da W solo attraverso z. Chiamiamo δ = ∂L/∂z (forma (m,), il gradiente che arriva dall'alto). Per la catena, ∂L/∂W_pq = Σᵢ δᵢ · ∂zᵢ/∂W_pq. Ma zᵢ = Σⱼ Wᵢⱼ xⱼ, quindi ∂zᵢ/∂W_pq è zero a meno che i = p, e in quel caso vale x_q. La somma collassa su un solo termine:

∂L/∂W_pq = δ_p x_q, cioè ∂L/∂W = δ xᵀ, il prodotto esterno fra il gradiente che scende e l'ingresso che era salito. ∎

Vale la pena fermarsi, perché è il cuore della backpropagation: il gradiente di un peso è il prodotto di due numeri soli, quello che è entrato in quel peso e quello che è uscito dall'errore. Niente Jacobiani giganti: un prodotto esterno di costo m·n, quanto il prodotto in avanti.

A = np.array([[1.0, 0.0],
              [1.0, 1.0],
              [0.0, 2.0]])
b = np.array([1.0, 1.0, 1.0])
x = np.array([1.0, 2.0])

r = A @ x - b
print("residuo r = Ax − b   :", riga(r))
print("‖Ax − b‖²            :", f"{float(r @ r):.4f}")
print("analitico 2Aᵀ(Ax − b):", riga(2 * A.T @ r))
print("numerico             :", riga(gradiente_numerico(lambda v: float((A @ v - b) @ (A @ v - b)), x)))

x_stella = np.linalg.solve(A.T @ A, A.T @ b)      # equazione normale AᵀA x = Aᵀb
print("soluzione di AᵀA x = Aᵀb:", riga(x_stella, 6))
print("il gradiente nella soluzione è nullo:", np.allclose(2 * A.T @ (A @ x_stella - b), 0.0))

W = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0]])                   # forma (2, 3)
x3 = np.array([1.0, 0.0, -1.0])                   # forma (3,)
delta = np.array([10.0, 20.0])                    # ∂L/∂(Wx), forma (2,)
print("forma di W:", W.shape, "| forma di δxᵀ:", np.outer(delta, x3).shape)
print("analitico δxᵀ:", riga(np.outer(delta, x3)))
print("numerico     :", riga(gradiente_numerico(lambda M: delta @ (M @ x3), W)))
print("xδᵀ avrebbe forma", np.outer(x3, delta).shape, "→ con W quadrata l'errore non si vedrebbe")

Output ottenuto:

residuo r = Ax − b   : +0.0000 +2.0000 +3.0000
‖Ax − b‖²            : 13.0000
analitico 2Aᵀ(Ax − b): +4.0000 +16.0000
numerico             : +4.0000 +16.0000
soluzione di AᵀA x = Aᵀb: +0.777778 +0.444444
il gradiente nella soluzione è nullo: True
forma di W: (2, 3) | forma di δxᵀ: (2, 3)
analitico δxᵀ: +10.0000 +0.0000 -10.0000 +20.0000 +0.0000 -20.0000
numerico     : +10.0000 +0.0000 -10.0000 +20.0000 +0.0000 -20.0000
xδᵀ avrebbe forma (3, 2) → con W quadrata l'errore non si vedrebbe

Cosa fa il codice. La prima parte verifica l'identità 3: residuo (0, 2, 3), perdita 13, gradiente analitico (4, 16) uguale al numerico; poi risolve l'equazione normale, e il punto (7/9, 4/9) annulla il gradiente — la definizione stessa di soluzione ai minimi quadrati. La seconda verifica l'identità 4 derivando rispetto a una matrice: gradiente_numerico perturba una casella alla volta di W e ricostruisce una tabella (2, 3) identica a δxᵀ. L'ultima riga stampa la forma della versione trasposta, (3, 2): incompatibile qui, perfettamente compatibile se W fosse quadrata.

Softmax ed entropia incrociata: da dove esce ŷ − y

L'ultima identità è quella che in ogni corso compare come dono del cielo. Il modello produce K punteggi grezzi z (i logit); la softmax li trasforma in probabilità ŷᵢ = e^{zᵢ} / Σ_j e^{z_j}; l'entropia incrociata con l'etichetta vera y (un vettore con un 1 e tutti 0) vale L = −Σᵢ yᵢ log ŷᵢ. Deriviamo L rispetto a z_k.

Primo passo: il logaritmo della softmax si semplifica. log ŷᵢ = zᵢ − log(Σ_j e^{z_j}), quindi L = −Σᵢ yᵢ zᵢ + (Σᵢ yᵢ) · log(Σ_j e^{z_j}). Poiché Σᵢ yᵢ = 1, resta L = −Σᵢ yᵢ zᵢ + log(Σ_j e^{z_j}).

Secondo passo: deriviamo i due addendi. Il primo è lineare, e dà −y_k. Il secondo è il logaritmo di una somma: ∂/∂z_k log(Σ_j e^{z_j}) = e^{z_k} / Σ_j e^{z_j} = ŷ_k. Sommando:

∂L/∂z = ŷ − y. ∎

In parole: il gradiente sui logit è la differenza fra quello che il modello ha detto e quello che doveva dire. Due conseguenze: non serve mai costruire il Jacobiano K × K della softmax, perché la composizione con l'entropia incrociata collassa in una sottrazione; e la somma delle componenti è zero, perché sia ŷ sia y sommano a 1 — se un logit sale, un altro scende.

def softmax(z):
    z = z - z.max()                              # stabilità numerica: vedi la lezione su IEEE 754
    e = np.exp(z)
    return e / e.sum()

z = np.array([1.0, 2.0, 0.0])
y = np.array([0.0, 1.0, 0.0])                    # la classe vera è la seconda
p = softmax(z)
perdita = lambda v: float(-np.log(softmax(v) @ y))

print("ŷ = softmax(z)  :", riga(p))
print("perdita −log ŷ_y:", f"{perdita(z):.4f}")
print("analitico ŷ − y :", riga(p - y))
print("numerico        :", riga(gradiente_numerico(perdita, z)))
print("la somma del gradiente è zero:", np.allclose((p - y).sum(), 0.0))

Output ottenuto:

ŷ = softmax(z)  : +0.2447 +0.6652 +0.0900
perdita −log ŷ_y: 0.4076
analitico ŷ − y : +0.2447 -0.3348 +0.0900
numerico        : +0.2447 -0.3348 +0.0900
la somma del gradiente è zero: True

Cosa fa il codice. La sottrazione del massimo dentro softmax non cambia il risultato matematico ma evita l'overflow di np.exp sui logit grandi: è l'attenzione alla virgola mobile della lezione sui numeri in virgola mobile. Il gradiente analitico ŷ − y coincide col numerico fino alla quarta cifra e somma a zero; la classe giusta ha gradiente negativo (il suo logit deve salire), le altre due positivo.

Trade-off. La versione fusa «softmax più entropia incrociata» è più stabile e più veloce dei due pezzi separati, ed è il motivo per cui le librerie espongono una funzione sola (in PyTorch è cross_entropy, che si aspetta i logit e non le probabilità). Il prezzo è non poter infilare niente in mezzo: per una calibrazione o una perdita pesata per classe devi tornare ai due pezzi e accettare la stabilità peggiore.

Il Jacobiano e il prodotto vettore-Jacobiano

Finora la funzione in cima era sempre scalare: una perdita. Ma uno strato di rete non produce un numero, produce un vettore. Serve l'oggetto che descrive la derivata di una funzione vettoriale.

Definizione. Data F: Rⁿ → Rᵐ, il Jacobiano in x è la matrice J di forma (m, n) con

Jᵢⱼ = ∂Fᵢ/∂xⱼ.

Leggila per righe e per colonne. La riga i è il gradiente dell'uscita i-esima: dice da quali ingressi dipende. La colonna j dice che effetto ha l'ingresso j-esimo su tutte le uscite insieme. Per m = 1 il Jacobiano ha una riga sola ed è il gradiente trasposto: il gradiente è un caso particolare del Jacobiano, non un oggetto diverso.

La catena come prodotto di Jacobiani

La regola della catena, in forma vettoriale, è una moltiplicazione di matrici. Se z = G(y) e y = F(x), allora

J_{G∘F}(x) = J_G(y) · J_F(x), con forme (p, m) · (m, n) = (p, n).

È la formula che giustifica l'intera backpropagation: una rete di L strati è una composizione di L funzioni, e la derivata della perdita è il prodotto di L Jacobiani. Le forme combaciano automaticamente, ed è un controllo gratuito: se non si incastrano, hai sbagliato l'ordine dei fattori. Il prodotto non è commutativo: scambiare due Jacobiani non dà «lo stesso risultato ma trasposto», dà un errore di forma o un numero sbagliato.

Il VJP, e perché nessuno materializza un Jacobiano

Qui la matematica si separa dall'ingegneria: scrivere la catena come prodotto di Jacobiani è corretto, ma calcolarli uno per uno e poi moltiplicarli sarebbe un suicidio di memoria — e non serve mai. Nella backpropagation il Jacobiano ti serve sempre moltiplicato a sinistra per un vettore, quello che arriva dallo strato successivo. L'operazione u ↦ uᵀJ si chiama prodotto vettore-Jacobiano (vector-Jacobian product, VJP) ed è l'unica cosa che una libreria implementa davvero. Per uno strato lineare z = Wx il Jacobiano rispetto a x è W stessa, e il VJP è uᵀW: un prodotto matrice-vettore di costo m·n, lo stesso del passaggio in avanti, non il suo quadrato. Il duale è il JVP (Jacobian-vector product), w ↦ Jw, e serve al modo diretto.

def G(v):
    """G: R² → R³. Tre uscite che dipendono dai due ingressi in modi diversi."""
    x, y = v
    return np.array([x**2 * y, x + y**3, x * y])

def jacobiano_analitico(v):
    x, y = v
    return np.array([[2*x*y, x**2],
                     [1.0,   3*y**2],
                     [y,     x]])

def jacobiano_numerico(f, v, h=1e-5):
    colonne = []
    for i in range(v.size):
        passo = np.zeros_like(v)
        passo[i] = h
        colonne.append((f(v + passo) - f(v - passo)) / (2 * h))
    return np.stack(colonne, axis=1)              # una colonna per ingresso

v = np.array([2.0, 3.0])
J = jacobiano_analitico(v)
print("forma di J:", J.shape, "(righe = uscite, colonne = ingressi)")
print("J analitico:", riga(J))
print("J numerico :", riga(jacobiano_numerico(G, v)))

u = np.array([1.0, 1.0, 1.0])
print("VJP uᵀJ, una passata all'indietro :", riga(u @ J))
print("JVP Jw con w = (1, 0), una in avanti:", riga(J @ np.array([1.0, 0.0])))

for n in (256, 1024, 4096, 16384):
    print(f"strato {n} → {n}: Jacobiano di {n*n} numeri = {n*n*4/2**20:.2f} MiB in float32")

Output ottenuto:

forma di J: (3, 2) (righe = uscite, colonne = ingressi)
J analitico: +12.0000 +4.0000 +1.0000 +27.0000 +3.0000 +2.0000
J numerico : +12.0000 +4.0000 +1.0000 +27.0000 +3.0000 +2.0000
VJP uᵀJ, una passata all'indietro : +16.0000 +33.0000
JVP Jw con w = (1, 0), una in avanti: +12.0000 +1.0000 +3.0000
strato 256 → 256: Jacobiano di 65536 numeri = 0.25 MiB in float32
strato 1024 → 1024: Jacobiano di 1048576 numeri = 4.00 MiB in float32
strato 4096 → 4096: Jacobiano di 16777216 numeri = 64.00 MiB in float32
strato 16384 → 16384: Jacobiano di 268435456 numeri = 1024.00 MiB in float32

Cosa fa il codice, e che cosa ci dice. jacobiano_numerico costruisce la matrice per colonne: perturba un ingresso alla volta e guarda come si muovono tutte le uscite insieme — è il modo diretto implementato in modo ingenuo, e costa una passata per ingresso. Il Jacobiano analitico coincide; il VJP contrae le tre uscite e restituisce due numeri, il JVP contrae i due ingressi e ne restituisce tre. L'ultimo ciclo conta i byte, ed è la parte che pesa: uno strato 4096 → 4096 ha un Jacobiano da 64 MiB in precisione singola, uno 16384 → 16384 arriva a un gigabyte, per un esempio e uno strato; con un minibatch da 32 e dieci strati si moltiplica per 320. Nessuna libreria lo fa, e non per pigrizia: il Jacobiano di uno strato lineare è la matrice dei pesi, che già hai in memoria, e il VJP la usa senza copiarla. La regola da portare via: non costruire mai un Jacobiano che ti serve solo moltiplicato.

Hessiano, curvatura e il passo di Newton

Il gradiente dice in che direzione la funzione sale, ma non dice niente su quanto rapidamente cambia quella direzione. Quell'informazione sta nelle derivate seconde.

Definizione. Per f: Rⁿ → R, l'Hessiano è la matrice n × n delle derivate seconde,

Hᵢⱼ = ∂²f/∂xᵢ∂xⱼ.

È il Jacobiano del gradiente: H = J_{∇f}. Con derivate seconde continue il teorema di Schwarz garantisce ∂²f/∂xᵢ∂xⱼ = ∂²f/∂xⱼ∂xᵢ, quindi H è simmetrica. Non è un dettaglio: all'Hessiano si applica il teorema spettrale della lezione su autovalori, SVD e PCA — autovalori reali, autovettori perpendicolari, H = QΛQᵀ.

Il test degli autovalori: minimo, massimo o sella

In un punto dove il gradiente è nullo, la natura del punto si legge negli autovalori dell'Hessiano. Il motivo è lo sviluppo di Taylor al secondo ordine: muovendosi di d dal punto critico, f(x + d) ≈ f(x) + ½ dᵀHd, e il segno di dᵀHd decide se ci si alza o si scende. Scrivendo d nella base degli autovettori, dᵀHd = Σᵢ λᵢ cᵢ²: ogni autovalore pesa il quadrato della componente lungo il suo autovettore.

• Tutti gli autovalori positivi: dᵀHd > 0 in ogni direzione, quindi si sale sempre. È un minimo locale.
• Tutti negativi: si scende sempre. È un massimo.
• Segni misti: si sale lungo certe direzioni e si scende lungo altre. È una sella, e nelle reti profonde sono la grande maggioranza dei punti critici.
• Qualche autovalore nullo: il secondo ordine non decide, servono le derivate successive. In pratica è una direzione piatta, e la discesa del gradiente ci si trascina.

Il rapporto λ_max/λ_min è il numero di condizionamento dell'Hessiano, e misura quanto la conca è schiacciata. È la stessa «valle stretta» della lezione sulla discesa del gradiente: più il rapporto è grande, più la discesa rimbalza fra le pareti invece di scendere lungo il fondo.

Il passo di Newton, derivato da Taylor

Se conosci la curvatura, puoi fare di meglio che muoverti un pochino nella direzione opposta al gradiente: puoi saltare al minimo dell'approssimazione quadratica. Sviluppiamo f al secondo ordine intorno a x:

f(x + d) ≈ f(x) + ∇f(x)ᵀd + ½ dᵀHd

Questa è una funzione di d di cui sappiamo derivare il gradiente con le identità 1 e 2: ∂/∂d = ∇f(x) + Hd, perché H è simmetrica e quindi ∂(½dᵀHd)/∂d = Hd. Annullandolo si ottiene Hd = −∇f(x), cioè

x ← x − H⁻¹∇f(x), il passo di Newton. ∎

In parole: invece di scegliere un passo a mano, lo si deduce dalla curvatura. Su una funzione esattamente quadratica un solo passo porta nel minimo, da qualunque punto si parta. E H⁻¹ non va invertita: si risolve Hd = −∇f, che costa meno ed è più stabile — la lezione dell'equazione normale, di nuovo.

def f2(v):
    """f(x₀, x₁) = 2x₀² + 3x₀x₁ + 4x₁² − 5x₀ − 6x₁"""
    return 2*v[0]**2 + 3*v[0]*v[1] + 4*v[1]**2 - 5*v[0] - 6*v[1]

def gradiente_f2(v):
    return np.array([4*v[0] + 3*v[1] - 5, 3*v[0] + 8*v[1] - 6])

def hessiano_numerico(f, x, h=1e-4):
    n = x.size
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            ei = np.zeros(n); ei[i] = h
            ej = np.zeros(n); ej[j] = h
            H[i, j] = (f(x+ei+ej) - f(x+ei-ej) - f(x-ei+ej) + f(x-ei-ej)) / (4*h*h)
    return H

H = np.array([[4.0, 3.0],
              [3.0, 8.0]])                        # ∂²f/∂xᵢ∂xⱼ, calcolate a mano
x0 = np.array([1.0, 1.0])
print("H analitico:", riga(H))
print("H numerico :", riga(hessiano_numerico(f2, x0)))
print("simmetrico, come vuole Schwarz:", np.allclose(H, H.T))

autov = np.linalg.eigvalsh(H)
print("autovalori di H:", riga(autov), "→ tutti positivi:", bool(np.all(autov > 0)), "→ minimo")
print("condizionamento λ_max/λ_min:", f"{autov[-1] / autov[0]:.4f}")
sella = lambda v: v[0]**2 - v[1]**2
print("autovalori di x₀² − x₁²:", riga(np.linalg.eigvalsh(hessiano_numerico(sella, x0))), "→ sella")

x = np.array([0.0, 0.0])
print("partenza:", riga(x, 6), "| f =", f"{f2(x):.6f}")
x_newton = x - np.linalg.solve(H, gradiente_f2(x))
print("dopo UN passo di Newton:", riga(x_newton, 6), "| f =", f"{f2(x_newton):.6f}")
print("il gradiente nel nuovo punto è nullo:", np.allclose(gradiente_f2(x_newton), 0.0))

eta = 0.1
contrazione = max(abs(1 - eta*autov[0]), abs(1 - eta*autov[-1]))
print(f"discesa del gradiente con η = {eta}: l'errore si contrae di {contrazione:.4f} per passo")
print("passi per ridurre l'errore di 1000 volte:", int(np.ceil(np.log(1e-3) / np.log(contrazione))))
for n in (8192, 10**6):
    print(f"n = {n}: Hessiano di {n*n} numeri ({n*n*4/2**30:.2f} GiB in float32), n³/3 = {n**3/3:.2e} operazioni")

Output ottenuto:

H analitico: +4.0000 +3.0000 +3.0000 +8.0000
H numerico : +4.0000 +3.0000 +3.0000 +8.0000
simmetrico, come vuole Schwarz: True
autovalori di H: +2.3944 +9.6056 → tutti positivi: True → minimo
condizionamento λ_max/λ_min: 4.0116
autovalori di x₀² − x₁²: -2.0000 +2.0000 → sella
partenza: +0.000000 +0.000000 | f = 0.000000
dopo UN passo di Newton: +0.956522 +0.391304 | f = -3.565217
il gradiente nel nuovo punto è nullo: True
discesa del gradiente con η = 0.1: l'errore si contrae di 0.7606 per passo
passi per ridurre l'errore di 1000 volte: 26
n = 8192: Hessiano di 67108864 numeri (0.25 GiB in float32), n³/3 = 1.83e+11 operazioni
n = 1000000: Hessiano di 1000000000000 numeri (3725.29 GiB in float32), n³/3 = 3.33e+17 operazioni

Cosa fa il codice. hessiano_numerico usa la formula a quattro punti — due coordinate perturbate alla volta, tutte le combinazioni di segno, diviso 4h² — e conferma l'Hessiano calcolato a mano. eigvalsh dà gli autovalori, entrambi positivi: siamo in una conca, e il condizionamento dice quanto è schiacciata; x₀² − x₁² fa da controprova con un autovalore per segno, cioè una sella. Infine il passo di Newton parte dall'origine e atterra nel minimo esatto (22/23, 9/23) in un passo, con gradiente nullo, mentre la discesa del gradiente con η = 0,1 contrae l'errore di max|1 − ηλᵢ| per passo e ne serve 26.

Perché il passo di Newton non si usa nel deep learning

Dopo una dimostrazione così, la domanda è obbligatoria: se Newton converge in un passo e la discesa del gradiente in decine, perché tutti usano la seconda? Le ragioni sono tre, e sono tutte quantitative.

La memoria. L'Hessiano ha n² caselle. Per la sola matrice W di forma (128, 64) della prossima lezione, n = 8192 e l'Hessiano occupa 0,25 GiB in precisione singola; per un modello da un milione di parametri — un modello piccolo — sono mille miliardi di numeri, quasi quattro terabyte. Non è questione di hardware migliore: è crescita quadratica contro crescita lineare.

Il tempo. Risolvere Hd = −∇f costa dell'ordine di n³/3 operazioni: con n attorno al milione sono 3·10¹⁷ operazioni per passo, contro le poche decine di milioni di un passo di discesa. Un passo migliore non vale centomila passi.

Il condizionamento. Anche potendo, l'Hessiano di una rete profonda è malcondizionato e spesso indefinito: in una sella ha autovalori negativi, e lì −H⁻¹∇f non punta in discesa ma verso l'alto. Newton, preso alla lettera, corre verso il punto critico più vicino, che in una rete profonda è quasi sempre una sella.

Trade-off. Quello che sopravvive del secondo ordine è l'idea, non la matrice: i metodi quasi-Newton (L-BFGS) ne tengono un'approssimazione a rango basso con poche decine di vettori, e gli ottimizzatori adattivi tipo Adam usano una stima diagonale della curvatura, cioè n numeri invece di n². Costano quasi come la discesa del gradiente e catturano la parte di curvatura che conta. Il prezzo è che l'approssimazione può essere sbagliata e introduce altri iperparametri da scegliere.

Modo diretto e modo inverso

L'analogia: la distinta base contro il preventivo

Un'azienda assembla un prodotto finito partendo da molte materie prime. Arriva la notizia che il rame è aumentato di un centesimo al chilo, e il direttore vuole sapere di quanto cresce il costo del prodotto finito. Si prende la distinta base e si propaga l'aumento in avanti: il rame entra nei cavi, i cavi nel modulo, il modulo nel prodotto. Una passata dalla materia prima al prodotto, e hai la risposta — per il rame. Se domani aumenta l'alluminio, rifai tutta la passata.

Il giorno dopo il direttore cambia domanda: «il prodotto finito è aumentato di un euro; di chi è la colpa?» Adesso la distinta base si percorre all'indietro: l'euro si ripartisce fra i sottoassiemi, ogni sottoassieme ripartisce la sua quota fra i componenti, e in una sola passata hai la responsabilità di tutte le materie prime insieme. È il preventivo rovesciato, e nel calcolo si chiama modo inverso.

L'analogia regge anche su un punto meno confortante: per ripartire all'indietro devi ricordarti com'era fatto l'assemblaggio. La distinta base in avanti si può percorrere dimenticando quello che c'è dietro; quella all'indietro richiede di aver conservato la struttura. È esattamente il costo di memoria del modo inverso, il famoso «nastro».

Il conto esatto, e da dove viene la scelta

Formalizziamo. Sia f: Rⁿ → Rᵐ valutabile con C operazioni elementari.

• Modo diretto: una passata propaga in avanti, accanto a ogni valore, la sua derivata rispetto a un ingresso scelto. Una passata dà una colonna del Jacobiano, cioè Jw per w versore. Per tutto il Jacobiano servono n passate: costo ≈ n · C.
• Modo inverso: una passata in avanti registra i valori, una passata all'indietro propaga la derivata di una uscita scelta rispetto a tutto. Una passata dà una riga del Jacobiano, cioè uᵀJ. Per tutto il Jacobiano servono m passate: costo ≈ m · C, con un fattore costante fra 2 e 4 e una memoria proporzionale alla lunghezza del calcolo.

Qui finisce la matematica e comincia la scelta ingegneristica. La funzione da differenziare è la perdita: un numero solo. Quindi m = 1 e n è il numero di parametri, fra i milioni e i miliardi. Il modo inverso costa come due o tre valutazioni della perdita e restituisce tutti i gradienti; il modo diretto ne costerebbe un miliardo. Non è una preferenza di stile: è un rapporto di un miliardo a uno.

Il simmetrico vale altrettanto: con un parametro solo e mille uscite (n = 1, m = 1000) il modo diretto sarebbe mille volte più conveniente e non avrebbe bisogno di nastro. Il modo inverso non è «migliore»: è migliore quando le uscite sono poche e gli ingressi molti.

Ora lo misuriamo. Prendiamo f(x) = Σᵢ xᵢ², con n ingressi e una uscita, la implementiamo nei due modi e contiamo le operazioni elementari. Contare è meglio che cronometrare: un cronometro misura anche il computer, un contatore misura solo l'algoritmo ed è riproducibile su qualunque macchina.

class Duale:
    """Numero duale a + b·ε con ε² = 0: porta il valore e UNA derivata direzionale."""

    def __init__(self, valore, derivata, conta):
        self.valore, self.derivata, self.conta = valore, derivata, conta

    def __add__(self, altro):
        self.conta[0] += 2                        # una per il valore, una per la derivata
        return Duale(self.valore + altro.valore, self.derivata + altro.derivata, self.conta)

    def __mul__(self, altro):
        self.conta[0] += 2
        return Duale(self.valore * altro.valore,
                     self.valore * altro.derivata + self.derivata * altro.valore, self.conta)

def modo_diretto(valori, i_seme, conta):
    """f(x) = Σ xᵢ²: una passata dà UNA derivata parziale, quella del seme."""
    s = Duale(0.0, 0.0, conta)
    for i, v in enumerate(valori):
        xi = Duale(v, 1.0 if i == i_seme else 0.0, conta)
        s = s + xi * xi
    return s.derivata

def modo_inverso(valori, conta):
    """La stessa f: una passata avanti e una all'indietro danno TUTTE le derivate."""
    s = 0.0
    for v in valori:                              # passata in avanti
        q = v * v;  conta[0] += 1
        s = s + q;  conta[0] += 1
    g = np.zeros(len(valori))
    aggiunto_s = 1.0                              # il seme: ∂f/∂f = 1
    for i, v in enumerate(valori):                # passata all'indietro
        aggiunto_q = aggiunto_s;      conta[0] += 1          # ∂s/∂qᵢ = 1
        g[i] = aggiunto_q * (2 * v);  conta[0] += 1          # ∂qᵢ/∂xᵢ = 2xᵢ
    return g

rng = np.random.default_rng(7)
print(f"{'n':>4} | {'diretto':>8} | {'inverso':>8} | {'rapporto':>8} | stesso gradiente")
for n in (2, 10, 100):
    valori = rng.normal(size=n)
    c_d = [0]; g_d = np.array([modo_diretto(valori, i, c_d) for i in range(n)])
    c_i = [0]; g_i = modo_inverso(valori, c_i)
    print(f"{n:>4} | {c_d[0]:>8} | {c_i[0]:>8} | {c_d[0]/c_i[0]:>8.1f} | {np.allclose(g_d, g_i)}")

Output ottenuto:

   n |  diretto |  inverso | rapporto | stesso gradiente
   2 |       16 |        8 |      2.0 | True
  10 |      400 |       40 |     10.0 | True
 100 |    40000 |      400 |    100.0 | True

Cosa fa il codice, e che cosa ci dice. Il numero duale è il modo diretto in forma di oggetto: ogni valore porta appiccicata la propria derivata rispetto all'ingresso seminato, e __add__ e __mul__ sono la regola della somma e quella del prodotto scritte una volta per tutte. Per il gradiente completo la passata va rifatta n volte; il modo inverso ne fa una in avanti e una all'indietro e produce tutte le derivate in un colpo. I due gradienti coincidono, ed è lo stesso oggetto calcolato in due ordini diversi. Il rapporto fra i due conteggi è esattamente n, per costruzione e non per caso. In cambio il modo inverso ha dovuto conservare qualcosa — qui s e il valore di ogni ingresso, su una rete profonda tutte le attivazioni intermedie. La regola da portare via: il modo inverso compra tempo con la memoria.

Un autograd da zero

Adesso abbiamo tutto per scrivere la cosa vera. Un sistema di differenziazione automatica a modo inverso ha bisogno di quattro ingredienti, e nessuno dei quattro è difficile: un oggetto che ricorda da dove viene, una regola locale per ogni operazione, un ordine in cui percorrere il grafo all'indietro, e un seme da cui partire.

La classe Nodo, e perché la chiusura

Ogni nodo conserva quattro cose: il dato (il valore calcolato in avanti), il grad (la derivata della perdita rispetto a quel valore, inizialmente zero), i genitori, e _indietro, una funzione senza argomenti che distribuisce il proprio gradiente ai genitori. L'ultima è il trucco che rende il codice breve: la regola di derivazione viene scritta mentre l'operazione viene eseguita, quando gli operandi sono ancora a portata di mano, e conservata come chiusura da chiamare più tardi. Nessuna tabella di derivate, nessun dispatch su stringhe: la derivata è già impacchettata nel nodo.

def somma_su_broadcast(grad, forma):
    """Rimette il gradiente nella forma dell'operando: somma gli assi aggiunti dal broadcast."""
    while grad.ndim > len(forma):
        grad = grad.sum(axis=0)
    for asse, dim in enumerate(forma):
        if dim == 1 and grad.shape[asse] != 1:
            grad = grad.sum(axis=asse, keepdims=True)
    return grad

class Nodo:
    """Un array che ricorda da dove viene: dato, gradiente, genitori, regola all'indietro."""

    def __init__(self, dato, genitori=(), operazione="foglia"):
        self.dato = np.asarray(dato, dtype=float)
        self.grad = np.zeros_like(self.dato)
        self.genitori = genitori
        self.operazione = operazione
        self._indietro = lambda: None             # le foglie non hanno niente da distribuire

    def __add__(self, altro):
        altro = altro if isinstance(altro, Nodo) else Nodo(altro)
        fuori = Nodo(self.dato + altro.dato, (self, altro), "+")
        def propaga():                            # la somma distribuisce il gradiente intatto
            self.grad += somma_su_broadcast(fuori.grad, self.dato.shape)
            altro.grad += somma_su_broadcast(fuori.grad, altro.dato.shape)
        fuori._indietro = propaga
        return fuori

    def __mul__(self, altro):
        altro = altro if isinstance(altro, Nodo) else Nodo(altro)
        fuori = Nodo(self.dato * altro.dato, (self, altro), "*")
        def propaga():                            # regola del prodotto, elemento per elemento
            self.grad += somma_su_broadcast(altro.dato * fuori.grad, self.dato.shape)
            altro.grad += somma_su_broadcast(self.dato * fuori.grad, altro.dato.shape)
        fuori._indietro = propaga
        return fuori

    def __matmul__(self, altro):
        fuori = Nodo(self.dato @ altro.dato, (self, altro), "@")
        def propaga():
            self.grad += fuori.grad @ altro.dato.T        # δWᵀ: il VJP dello strato lineare
            altro.grad += self.dato.T @ fuori.grad        # Xᵀδ: il prodotto esterno, identità 4
        fuori._indietro = propaga
        return fuori

Come leggere il codice. Le tre operazioni sono le tre righe di matematica già derivate. La somma ha derivata 1 rispetto a ciascun addendo e passa il gradiente intatto; l'unica complicazione è il broadcast, perché sommando un bias (16,) a una matrice (600, 16) il gradiente torna con 600 righe e va ricompattato sommando l'asse aggiunto. Il prodotto matriciale è l'identità 4 e il suo gemello: gradiente sull'ingresso δWᵀ, sui pesi Xᵀδ, il prodotto esterno generalizzato al minibatch. E i gradienti si sommano (+=), non si assegnano: se un nodo è usato due volte i due contributi devono accumularsi.

Gli operatori, l'ordinamento topologico e il seme

Mancano le non linearità, la riduzione finale e il motore che percorre il grafo. L'ordine non è libero: un nodo può distribuire il proprio gradiente solo dopo aver ricevuto i contributi di tutti i figli. L'ordine giusto è l'ordinamento topologico rovesciato, cioè quello che produce una visita in profondità che accoda un nodo dopo i suoi genitori.

    def relu(self):
        fuori = Nodo(np.maximum(self.dato, 0.0), (self,), "relu")
        def propaga():
            self.grad += (self.dato > 0) * fuori.grad     # in 0 la derivata si SCEGLIE: qui 0
        fuori._indietro = propaga
        return fuori

    def tanh(self):
        t = np.tanh(self.dato)
        fuori = Nodo(t, (self,), "tanh")
        def propaga():
            self.grad += (1.0 - t**2) * fuori.grad        # tanh' = 1 − tanh²
        fuori._indietro = propaga
        return fuori

    def log(self):
        fuori = Nodo(np.log(self.dato), (self,), "log")
        def propaga():
            self.grad += fuori.grad / self.dato
        fuori._indietro = propaga
        return fuori

    def sum(self):
        fuori = Nodo(self.dato.sum(), (self,), "sum")
        def propaga():
            self.grad += np.ones_like(self.dato) * fuori.grad
        fuori._indietro = propaga
        return fuori

    def indietro(self):
        """Ordinamento topologico, poi il grafo percorso all'inverso partendo dal seme."""
        ordine, visti = [], set()
        def visita(nodo):
            if id(nodo) in visti:
                return
            visti.add(id(nodo))
            for genitore in nodo.genitori:
                visita(genitore)
            ordine.append(nodo)                   # dopo i genitori: ordine topologico
        visita(self)
        self.grad = np.ones_like(self.dato)       # il seme: ∂L/∂L = 1
        for nodo in reversed(ordine):
            nodo._indietro()

def entropia_incrociata(logit, y_onehot):
    """Softmax ed entropia incrociata in un solo nodo: la derivata è (ŷ − y)/m, identità 5."""
    z = logit.dato - logit.dato.max(axis=1, keepdims=True)
    e = np.exp(z)
    p = e / e.sum(axis=1, keepdims=True)
    m = p.shape[0]
    fuori = Nodo(-np.log((p * y_onehot).sum(axis=1)).mean(), (logit,), "entropia_incrociata")
    def propaga():
        logit.grad += (p - y_onehot) / m * fuori.grad
    fuori._indietro = propaga
    return fuori

Come leggere il codice. Le quattro non linearità seguono lo stesso schema: calcola il valore in avanti, costruisci il nodo, impacchetta la derivata locale in una chiusura. relu merita una nota: in 0 non è derivabile, e il codice sceglie 0 scrivendo self.dato > 0 invece di >= — una convenzione, non un teorema, ed è proprio lì che il controllo numerico fallirà. Il metodo indietro fa tre cose in dieci righe: visita il grafo accodando ogni nodo dopo i suoi genitori, mette il seme 1 sulla radice, chiama le chiusure in ordine rovesciato. entropia_incrociata è un nodo fuso: usa direttamente l'identità 5 invece di derivare softmax e logaritmo separatamente.

Proviamolo su un esempio verificabile a mano: X = (1, −2, 3) come riga, W di forma (3, 2), perdita L = somma(relu(XW)). Il prodotto dà Z = (4, 0): primo canale attivo, secondo esattamente sul ginocchio della ReLU. Il gradiente a monte è quindi δ = (1, 0), e le due identità danno ∂L/∂W = Xᵀδ e ∂L/∂X = δWᵀ: numeri che si contano sulle dita.

X = Nodo([[1.0, -2.0, 3.0]])                      # (1, 3)
W = Nodo([[1.0, 2.0],
          [0.0, 1.0],
          [1.0, 0.0]])                            # (3, 2)
L = (X @ W).relu().sum()
L.indietro()
print("Z = X @ W :", riga(X.dato @ W.dato))
print("L         :", f"{float(L.dato):.4f}")
print("∂L/∂W     :", riga(W.grad), "  (atteso Xᵀδ con δ = (1, 0))")
print("∂L/∂X     :", riga(X.grad), "  (atteso δWᵀ)")

a = Nodo(3.0)
b = a * a + a                                     # a compare due volte: i contributi si sommano
b.indietro()
print("d(a² + a)/da con a = 3:", f"{float(a.grad):.4f}", "(atteso 2a + 1 = 7)")

Output ottenuto:

Z = X @ W : +4.0000 +0.0000
L         : 4.0000
∂L/∂W     : +1.0000 +0.0000 -2.0000 +0.0000 +3.0000 +0.0000   (atteso Xᵀδ con δ = (1, 0))
∂L/∂X     : +1.0000 +0.0000 +1.0000   (atteso δWᵀ)
d(a² + a)/da con a = 3: 7.0000 (atteso 2a + 1 = 7)

Cosa fa il codice. Costruisce il grafo della figura e lo percorre all'indietro una volta. Il gradiente su W è il prodotto esterno fra ingresso e δ: prima colonna (1, −2, 3), seconda tutta nulla perché la ReLU ha chiuso il secondo canale. Quello su X è δWᵀ. Il secondo esempio verifica l'accumulo: a entra due volte in a² + a, e il += somma i due contributi ottenendo 2a + 1 = 7. Con = al posto di += uno dei due andrebbe perso e il risultato sarebbe 6 o 1 secondo l'ordine di visita: un bug che dipende dall'ordine è il peggiore.

Il gradient check fatto bene

Un autograd scritto a mano è codice come un altro, e il suo output è un numero che nessuno sa leggere a occhio. Serve un oracolo indipendente, e l'unico disponibile senza installare niente sono le differenze finite. Ma vanno usate con tre accortezze, e ciascuna delle tre si può derivare.

Prima accortezza: differenze centrate, non in avanti. Sviluppiamo f in serie di Taylor nei due punti x ± h:

f(x + h) = f(x) + h f'(x) + (h²/2) f''(x) + (h³/6) f'''(x) + …

f(x − h) = f(x) − h f'(x) + (h²/2) f''(x) − (h³/6) f'''(x) + …

Sottraendo, i termini pari si cancellano: f(x + h) − f(x − h) = 2h f'(x) + (h³/3) f'''(x) + …, da cui

(f(x + h) − f(x − h)) / 2h = f'(x) + (h²/6) f'''(x) + …

L'errore è O(h²). La differenza in avanti (f(x + h) − f(x))/h, invece, non cancella il termine in f'' e sbaglia di O(h). A parità di h la versione centrata costa una valutazione in più e guadagna cinque ordini di grandezza con h = 1e-5. È il miglior affare della lezione. ∎

Quanto deve valere h

Qui si scontrano due errori che vanno in direzioni opposte. Il troncamento è quello appena calcolato, ≈ (h²/6)·|f'''|, e diminuisce con h. L'arrotondamento è quello della lezione sui numeri in virgola mobile: f(x + h) e f(x − h) sono due numeri vicinissimi, e sottrarli è una cancellazione catastrofica; l'incertezza relativa su ciascuno è l'epsilon di macchina ε ≈ 2,2·10⁻¹⁶, quindi l'errore assoluto sulla differenza è ≈ ε·|f| e, diviso per 2h, cresce come ε·|f|/h.

Come leggere la tabella, riga per riga. Gli ordini di grandezza vengono dalle due formule, con |f| ≈ 1 e |f'''| ≈ 6. Il minimo della somma si trova derivando: d/dh (h²/3 + ε/h) = 0 dà h ≈ ∛ε ≈ 5·10⁻⁶. Da qui la ricetta h = 1e-5, che non è superstizione ma la radice cubica dell'epsilon di macchina — e da qui il controsenso apparente: rimpicciolire h oltre quel punto peggiora la stima.

Seconda accortezza: il criterio è l'errore relativo. Confrontare |numerico − analitico| con una soglia assoluta non ha senso, perché i gradienti di strati diversi differiscono di ordini di grandezza. Si usa |a − b| / max(|a|, |b|, η) con un piccolo η a proteggere la divisione per zero, e si accetta sotto 10⁻⁶.

Terza accortezza: non controllare dove la funzione non è derivabile. È l'errore che vedremo nell'esperimento, e vale la pena anticiparlo: su una ReLU valutata esattamente in 0, le differenze finite restituiscono un numero sensato e diverso da qualunque convenzione analitica.

def gradient_check(perdita_di, parametri, h=1e-5, eta=1e-12):
    """Confronta i gradienti già in .grad con le differenze centrate. Criterio: errore RELATIVO."""
    peggiore = 0.0
    for p in parametri:
        for i in range(p.dato.size):
            originale = p.dato.flat[i]
            p.dato.flat[i] = originale + h; piu = perdita_di()
            p.dato.flat[i] = originale - h; meno = perdita_di()
            p.dato.flat[i] = originale                       # SEMPRE ripristinare
            numerico = (piu - meno) / (2 * h)
            analitico = float(p.grad.flat[i])
            scala = max(abs(numerico), abs(analitico), eta)
            peggiore = max(peggiore, abs(numerico - analitico) / scala)
    return peggiore

# (a) perdita liscia: L = Σ (XW)², un polinomio di secondo grado nei pesi
X = Nodo([[1.0, -2.0, 3.0]])
W = Nodo([[1.0, 2.0], [0.0, 1.0], [1.0, 0.0]])
Z = X @ W
L = (Z * Z).sum()
L.indietro()
perdita_liscia = lambda: float(((X.dato @ W.dato) ** 2).sum())
print("(a) perdita liscia L = Σ (XW)²")
print("    analitico:", riga(W.grad))
print("    numerico :", riga(gradiente_numerico(lambda M: float(((X.dato @ M) ** 2).sum()), W.dato)))
print("    errore relativo massimo:", f"{gradient_check(perdita_liscia, [W]):.2e}")

# (b) la stessa rete con la ReLU, valutata ESATTAMENTE sul ginocchio
X = Nodo([[1.0, -2.0, 3.0]])
W = Nodo([[1.0, 2.0], [0.0, 1.0], [1.0, 0.0]])
L = (X @ W).relu().sum()
L.indietro()
perdita_relu = lambda: float(np.maximum(X.dato @ W.dato, 0.0).sum())
print("(b) perdita con ReLU, nel punto di non derivabilità")
print("    analitico:", riga(W.grad))
print("    numerico :", riga(gradiente_numerico(lambda M: float(np.maximum(X.dato @ M, 0.0).sum()), W.dato)))
print("    errore relativo massimo:", f"{gradient_check(perdita_relu, [W]):.4f}")

# (c) la U del passo: troncamento a sinistra, cancellazione a destra, su f(t) = t³ in t = 1
print("(c) scelta di h su f(t) = t³ in t = 1, derivata esatta 3")
for h in (1e-1, 1e-2, 1e-3, 1e-5, 1e-7, 1e-9, 1e-11, 1e-13):
    stima = ((1.0 + h) ** 3 - (1.0 - h) ** 3) / (2 * h)
    print(f"    h = {h:.0e}: stima {stima:.12f} | errore relativo {abs(stima - 3.0) / 3.0:.3e}")

Output ottenuto:

(a) perdita liscia L = Σ (XW)²
    analitico: +8.0000 +0.0000 -16.0000 +0.0000 +24.0000 +0.0000
    numerico : +8.0000 +0.0000 -16.0000 +0.0000 +24.0000 +0.0000
    errore relativo massimo: 2.68e-11
(b) perdita con ReLU, nel punto di non derivabilità
    analitico: +1.0000 +0.0000 -2.0000 +0.0000 +3.0000 +0.0000
    numerico : +1.0000 +0.5000 -2.0000 -1.0000 +3.0000 +1.5000
    errore relativo massimo: 1.0000
(c) scelta di h su f(t) = t³ in t = 1, derivata esatta 3
    h = 1e-01: stima 3.010000000000 | errore relativo 3.333e-03
    h = 1e-02: stima 3.000100000000 | errore relativo 3.333e-05
    h = 1e-03: stima 3.000001000000 | errore relativo 3.333e-07
    h = 1e-05: stima 3.000000000097 | errore relativo 3.246e-11
    h = 1e-07: stima 3.000000000086 | errore relativo 2.876e-11
    h = 1e-09: stima 3.000000081688 | errore relativo 2.723e-08
    h = 1e-11: stima 3.000000248221 | errore relativo 8.274e-08
    h = 1e-13: stima 2.999267501025 | errore relativo 2.442e-04

Cosa fa il codice. Il blocco (a) controlla l'autograd su una perdita liscia: analitico 2·Xᵀ(XW) e numerico coincidono, con errore dell'ordine dell'arrotondamento. Il blocco (b) ripete il controllo con la ReLU su un ingresso che manda il secondo canale esattamente in 0: l'autograd risponde 0 per convenzione, le differenze centrate con la media fra la pendenza a sinistra (zero) e quella a destra, e l'errore relativo arriva a 1. Non è un bug dell'autograd, è una domanda mal posta. Il blocco (c) mostra la U del passo.

Che cosa ci dice la scansione di h. Da 1e-1 a 1e-3 l'errore relativo cala esattamente come h²: 3,3·10⁻³, poi 3,3·10⁻⁵, poi 3,3·10⁻⁷, cento volte più piccolo a ogni passo — è la previsione di Taylor, confermata cifra per cifra. Poi la discesa si ferma: fra 1e-5 e 1e-7 l'errore è piatto intorno a 3·10⁻¹¹, perché il troncamento è ormai sotto l'arrotondamento e non comanda più. Da lì in poi si risale: 2,7·10⁻⁸ a 1e-9, 8,3·10⁻⁸ a 1e-11, e 2,4·10⁻⁴ a 1e-13, cioè quasi mille volte peggio del minimo. Si noti che a h = 1e-13 la stima è 2,999267…: ha perso quattro cifre significative. La regola da portare via: il fondo della U è largo tre ordini di grandezza, e qualunque h fra 10⁻⁶ e 10⁻⁴ va bene; sotto 10⁻⁸ stai misurando il rumore.

L'esperimento: una spirale a due bracci

Serve un problema che una retta non possa risolvere, altrimenti il gradiente non viene messo alla prova. Costruiamo due spirali intrecciate: due bracci che partono dall'origine sfasati di mezzo giro. I dati sono simulati da noi con numpy.random.default_rng, e lo dichiariamo: quello che misuriamo è il nostro codice, non la natura.

La rete è 2→16→16→2 con tanh nei due strati nascosti ed entropia incrociata in uscita: 354 parametri, zero librerie di deep learning. L'ultimo strato parte a zero, e l'effetto si prevede esattamente: con i logit nulli la softmax dà (0,5; 0,5) per ogni esempio e la perdita iniziale vale log 2 = 0,693147, qualunque cosa facciano i primi due strati. È il primo collaudo: se quel numero non esce, c'è un errore nel passaggio in avanti e non serve guardare i gradienti.

def spirale_a_due_bracci(per_braccio=300, rumore=0.18, seme=42):
    """Due bracci intrecciati. Dati SIMULATI da noi, non un dataset reale."""
    rng_dati = np.random.default_rng(seme)
    t = np.linspace(0.3, 3.2, per_braccio)
    blocchi, etichette = [], []
    for classe in (0, 1):
        angolo = 2.5 * t + classe * np.pi         # il secondo braccio è sfasato di mezzo giro
        px = t * np.cos(angolo) + rng_dati.normal(0.0, rumore, per_braccio)
        py = t * np.sin(angolo) + rng_dati.normal(0.0, rumore, per_braccio)
        blocchi.append(np.stack([px, py], axis=1))
        etichette.append(np.full(per_braccio, classe))
    return np.concatenate(blocchi), np.concatenate(etichette)

X_dati, y_dati = spirale_a_due_bracci()
Y = np.eye(2)[y_dati]                             # etichette in forma one-hot, (600, 2)
print("dati simulati: X", X_dati.shape, "| y", y_dati.shape,
      "| classe 0:", int((y_dati == 0).sum()), "| classe 1:", int((y_dati == 1).sum()))

from sklearn.linear_model import LogisticRegression          # solo per misurare il muro lineare
lineare = LogisticRegression().fit(X_dati, y_dati)
print("accuratezza di un modello LINEARE sugli stessi dati:", f"{lineare.score(X_dati, y_dati):.4f}")

rng = np.random.default_rng(0)
def strato(n_in, n_out, nullo=False):
    pesi = np.zeros((n_in, n_out)) if nullo else rng.normal(0.0, np.sqrt(2.0 / n_in), (n_in, n_out))
    return Nodo(pesi), Nodo(np.zeros(n_out))

W1, b1 = strato(2, 16)
W2, b2 = strato(16, 16)
W3, b3 = strato(16, 2, nullo=True)                # ultimo strato a zero: logit nulli all'inizio
parametri = [W1, b1, W2, b2, W3, b3]
print("parametri della rete 2→16→16→2:", sum(p.dato.size for p in parametri))

def avanti(X):
    H1 = (Nodo(X) @ W1 + b1).tanh()
    H2 = (H1 @ W2 + b2).tanh()
    return H2 @ W3 + b3

perdita = entropia_incrociata(avanti(X_dati), Y)
print("perdita iniziale:", f"{float(perdita.dato):.6f}", "| log 2 =", f"{float(np.log(2.0)):.6f}")
print("accuratezza iniziale:", f"{float((avanti(X_dati).dato.argmax(axis=1) == y_dati).mean()):.4f}")
perdita.indietro()
print("forme dei gradienti:", [p.grad.shape for p in parametri])
print("con l'ultimo strato nullo, ∂L/∂W1 è esattamente zero:", bool(np.abs(W1.grad).max() == 0.0))

Output ottenuto:

dati simulati: X (600, 2) | y (600,) | classe 0: 300 | classe 1: 300
accuratezza di un modello LINEARE sugli stessi dati: 0.5567
parametri della rete 2→16→16→2: 354
perdita iniziale: 0.693147 | log 2 = 0.693147
accuratezza iniziale: 0.5000
forme dei gradienti: [(2, 16), (16,), (16, 16), (16,), (16, 2), (2,)]
con l'ultimo strato nullo, ∂L/∂W1 è esattamente zero: True

Cosa fa il codice, e che cosa ci dice. La funzione genera 300 punti per braccio lungo una spirale di raggio crescente, con rumore di deviazione standard 0,18 — abbastanza perché i bracci si sfiorino senza sovrapporsi; la regressione logistica misura il muro, cioè quanto poco si ottiene senza non linearità. Poi costruisce la rete e verifica la perdita iniziale contro log 2. Qui c'è una lezione dentro la lezione: con l'ultimo strato nullo il gradiente su W1 è esattamente zero, non approssimativamente, perché il gradiente verso il basso passa per W3 e se W3 = 0 non passa niente. Ecco perché si inizializza a zero solo l'ultimo strato, e solo per questo collaudo: tenerla come inizializzazione definitiva congelerebbe i due strati nascosti al primo passo. Un'inizializzazione tutta nulla non è «una scelta neutra»: è una rete che non impara mai.

La trasposizione sbagliata, fatta apposta

Adesso diamo all'ultimo strato pesi veri, verifichiamo tutti e 354 i gradienti con il gradient check, e introduciamo deliberatamente l'errore della prima sezione: scambiare Xᵀδ con la sua trasposta. Su W2, che è 16×16, il programma non protesterà.

W3.dato = rng.normal(0.0, np.sqrt(2.0 / 16), (16, 2))        # l'ultimo strato non è più nullo
for p in parametri:
    p.grad = np.zeros_like(p.dato)                           # azzerare PRIMA di ogni passata
perdita = entropia_incrociata(avanti(X_dati), Y)
perdita.indietro()

perdita_corrente = lambda: float(entropia_incrociata(avanti(X_dati), Y).dato)
errore = gradient_check(perdita_corrente, parametri)
print("gradient check su tutti i 354 parametri: errore relativo massimo", f"{errore:.3e}")
print("sotto la soglia di 1e-6:", errore < 1e-6)

giusto = W2.grad.copy()
print("la versione trasposta ha la stessa forma:", giusto.T.shape == giusto.shape)
print("ma è lo stesso gradiente?", bool(np.allclose(giusto, giusto.T)))
print("scarto massimo fra il gradiente e la sua trasposta:", f"{np.abs(giusto - giusto.T).max():.3e}")
print("su W1, che è (2, 16), la trasposizione dà subito:", end=" ")
try:
    W1.dato - W1.grad.T
except ValueError as errore_di_forma:
    print(type(errore_di_forma).__name__, "— le forme non combaciano")

for epoca in range(1, 3001):                                 # addestramento col SOLO autograd a mano
    perdita = entropia_incrociata(avanti(X_dati), Y)
    for p in parametri:
        p.grad = np.zeros_like(p.dato)
    perdita.indietro()
    for p in parametri:
        p.dato -= 0.5 * p.grad
    if epoca % 500 == 0:
        previste = avanti(X_dati).dato.argmax(axis=1)
        print(f"epoca {epoca:>4}: perdita {float(perdita.dato):.4f} | accuratezza {(previste == y_dati).mean():.4f}")

Output ottenuto:

gradient check su tutti i 354 parametri: errore relativo massimo 1.078e-07
sotto la soglia di 1e-6: True
la versione trasposta ha la stessa forma: True
ma è lo stesso gradiente? False
scarto massimo fra il gradiente e la sua trasposta: 1.074e-01
su W1, che è (2, 16), la trasposizione dà subito: ValueError — le forme non combaciano
epoca  500: perdita 0.0582 | accuratezza 0.9917
epoca 1000: perdita 0.0117 | accuratezza 0.9967
epoca 1500: perdita 0.0069 | accuratezza 0.9983
epoca 2000: perdita 0.0047 | accuratezza 0.9983
epoca 2500: perdita 0.0035 | accuratezza 1.0000
epoca 3000: perdita 0.0027 | accuratezza 1.0000

Cosa fa il codice, e che cosa ci dice. Il blocco rimette pesi casuali nell'ultimo strato, azzera i gradienti, rifà una passata completa e poi controlla tutti e 354 i parametri contro le differenze centrate. L'errore relativo massimo è 1,078·10⁻⁷: autograd e oracolo numerico concordano su sette cifre significative, ed è la prova che il codice è giusto — non «sembra funzionare», è verificato. La parte centrale documenta la trasposizione: su W2 la forma coincide e nessuno protesta, ma il gradiente differisce dal suo trasposto fino a 0,107, enorme rispetto ai gradienti stessi, che hanno modulo massimo nell'ordine di 10⁻⁴; su W1, rettangolare, NumPy solleva subito un ValueError. L'ultima parte addestra per 3000 epoche con i soli gradienti del nostro autograd: dalla perdita iniziale 0,693147 e dall'accuratezza 0,5000 si arriva a 1,0000 sui 600 punti, contro lo 0,5567 del miglior separatore lineare. Quella differenza fra 0,56 e 1,00 è il valore aggiunto di due strati non lineari, ottenuto con centoventi righe di NumPy. L'accuratezza è misurata sugli stessi dati dell'addestramento e non dice nulla sulla generalizzazione: serve solo a mostrare che l'ottimizzazione funziona.

Trade-off. Il gradient check su tutti i parametri costa 2 × 354 passaggi in avanti: qui pochi secondi, su una rete vera impraticabile. Nella pratica si controlla un campione casuale (venti o trenta parametri bastano), su un minibatch minuscolo, e una volta sola, appena scritto il codice: non è un test di regressione da lasciare acceso, è un collaudo di messa in opera.

Errori tipici (sintomo, causa, rimedio)

1. Il gradiente trasposto che converge lo stesso. Sintomo: la rete impara peggio del previsto, e solo sugli strati quadrati. Causa: xδᵀ al posto di δxᵀ. Rimedio: gradient check su un minibatch piccolo, non controllo delle forme.
2. h troppo piccolo nel gradient check. Sintomo: l'errore relativo cresce man mano che rimpicciolisci il passo. Causa: cancellazione catastrofica nella sottrazione. Rimedio: h = 1e-5, e se il test non passa cerca il bug, non un h più piccolo.
3. Mutazione in posto che spezza il grafo. Sintomo: gradienti nulli o sbagliati, senza nessun errore. Causa: la chiusura all'indietro usa il valore attuale, non quello di allora. Rimedio: aggiornare i parametri solo dopo indietro(), e ricostruire il grafo a ogni passata.
4. Gradienti non azzerati fra due passi. Sintomo: la perdita diverge dopo poche epoche. Causa: il += dentro le chiusure accumula sulla passata precedente. Rimedio: p.grad = np.zeros_like(p.dato) prima di ogni indietro() — è lo zero_grad() delle librerie, e non è una formalità.
5. Derivata della ReLU in zero. Sintomo: il gradient check fallisce su qualche peso e passa su tutti gli altri. Causa: in 0 le differenze finite restituiscono la media delle due pendenze. Rimedio: perturbare i dati di un epsilon, o controllare con tanh; non «aggiustare» l'autograd.
6. Materializzare un Jacobiano. Sintomo: memoria esaurita su una rete che in avanti girava benissimo. Causa: aver costruito J invece di implementare uᵀJ. Rimedio: pensare sempre in termini di VJP.

Le sei righe hanno in comune una cosa: nessuna produce un errore esplicito. Il gradiente trasposto, il passo troppo piccolo, la mutazione in posto e i gradienti non azzerati danno un numero sbagliato ma plausibile, e per questo il gradient check vale più di qualunque revisione del codice a occhio. Le ultime due, la ReLU in zero e il Jacobiano materializzato, sono invece falsi allarmi o costi nascosti: la prima fa fallire un controllo corretto, la seconda esaurisce la memoria su una rete che in avanti girava.

Best practice

• Dichiara la convenzione di layout in testa al file, in un commento, e non cambiarla mai a metà progetto.
• Scrivi il gradient check prima del ciclo di addestramento, non dopo il primo risultato strano.
• Controlla le forme con un assert all'ingresso delle funzioni di gradiente: assert grad_W.shape == W.shape costa nulla e intercetta la trasposizione.
• Verifica con due oracoli indipendenti: la formula chiusa derivata a mano e le differenze centrate. Se coincidono entrambe col tuo codice, il codice è giusto.
• Azzera i gradienti prima di ogni passata all'indietro, e rendilo parte del ciclo, non un gesto occasionale.
• Dichiara seme e versioni (default_rng(42), NumPy 2.3.5): un gradient check non riproducibile non è una verifica.

Le sei regole si leggono in due gruppi. Le prime tre (convenzione dichiarata, controllo scritto prima dell'addestramento, assert sulle forme) prevengono l'errore silenzioso più costoso della lezione; le ultime tre (due oracoli, gradienti azzerati, seme e versioni dichiarati) rendono la verifica ripetibile da chiunque. Se ne puoi adottare una sola oggi, scegli il controllo sulle forme: costa una riga.

Glossario minimo

• Gradiente (gradient): il vettore delle derivate parziali di uno scalare; ha la forma dell'oggetto rispetto a cui si deriva.
• Jacobiano (Jacobian): la matrice (m, n) delle derivate di una funzione vettoriale; una riga per uscita.
• VJP (vector-Jacobian product): l'operazione u ↦ uᵀJ, l'unica che serve davvero al modo inverso; il duale w ↦ Jw è il JVP.
• Hessiano (Hessian): la matrice simmetrica delle derivate seconde; i suoi autovalori dicono minimo, massimo o sella.
• Modo diretto (forward mode): propaga le derivate insieme ai valori; costa n passate, non serve memoria.
• Modo inverso (reverse mode): una passata in avanti e una all'indietro; costa m passate e la memoria del nastro.
• Nastro (tape): la registrazione delle operazioni eseguite, necessaria per percorrere il grafo all'indietro.
• Seme (seed della derivazione): il valore 1 messo sulla radice prima di scendere; è ∂L/∂L.

Gli otto termini si tengono insieme con una regola sola: ogni oggetto ha una forma, e ogni operazione dichiara quale forma restituisce. Gradiente e Jacobiano differiscono per quella forma, VJP e JVP per il lato da cui si moltiplica, modo diretto e inverso per il verso in cui si percorre il grafo. Se ricordi le forme, il resto del glossario discende da lì.

Mini-progetto

Estendi l'autograd della lezione e verificalo. Tutto con NumPy, niente librerie di deep learning.

1. Aggiungi al Nodo gli operatori __sub__, __truediv__ e __pow__ con esponente costante, scrivendo per ciascuno la derivata locale nella chiusura; aggiungi sigmoide, e deriva a mano che σ' = σ(1 − σ);
2. scrivi controlla(nodo_perdita, parametri, quanti=20), che campiona quanti parametri a caso (con seme dichiarato) invece di controllarli tutti e restituisce l'errore relativo peggiore con l'indice del parametro colpevole;
3. rompi deliberatamente una derivata locale — per esempio 1 + t**2 invece di 1 − t**2 nel tanh — e verifica che il controllo trovi l'errore e indichi lo strato giusto. Documenta il caso limite: che cosa succede se il parametro campionato ha gradiente esattamente nullo?

Consegna: uno script che stampa l'errore relativo peggiore col codice corretto, quello col tanh rotto, e il parametro su cui il controllo si è acceso.

Challenge

• Easy. Carta e penna: deriva ∂(xᵀAx)/∂x componente per componente per n = 2 con A = [[2, 1], [0, 3]], e verifica che il risultato sia (A + Aᵀ)x e non 2Ax. Soluzione: f = 2x₀² + x₀x₁ + 3x₁², quindi ∂f/∂x₀ = 4x₀ + x₁ e ∂f/∂x₁ = x₀ + 6x₁, cioè [[4, 1], [1, 6]]x = (A + Aᵀ)x. Con 2Ax si otterrebbe [[4, 2], [0, 6]]x: sbagliato su entrambe le componenti, con la forma giusta.
• Medium. Carta e penna: dimostra ∂‖Ax − b‖²/∂x = 2Aᵀ(Ax − b) espandendo (Ax − b)ᵀ(Ax − b), poi mostra che annullare il gradiente dà AᵀAx = Aᵀb. Soluzione: f = xᵀAᵀAx − 2bᵀAx + bᵀb; il primo termine è una forma quadratica con AᵀA simmetrica, quindi deriva in 2AᵀAx; il secondo è lineare e dà −2Aᵀb; il terzo è costante. Sommando, 2Aᵀ(Ax − b). Annullando si ottiene l'equazione normale, con soluzione unica quando le colonne di A sono indipendenti.
• Hard. Implementa il modo diretto con i numeri duali su una funzione f: Rⁿ → R a tua scelta e misura che il costo del gradiente completo cresce linearmente in n, contando le operazioni elementari per n = 2, 5, 10, 20, 50, 100. Confronta col modo inverso e verifica che il rapporto sia esattamente n.
• Expert. Calcola una riga dell'Hessiano della perdita della spirale rispetto ai pesi dell'ultimo strato differenziando due volte: ottieni ∂L/∂W3 con l'autograd, poi applica le differenze centrate al gradiente invece che alla perdita, (∇L(W + h e_k) − ∇L(W − h e_k)) / 2h. Confronta con la stessa riga calcolata dalla formula a quattro punti sulla perdita, e spiega perché la prima è più accurata a parità di h.
• Enterprise. Prendi una funzione di costo del tuo lavoro — un prezzo che dipende da una tabella di parametri, un punteggio di rischio — e derivala a mano rispetto a tutti i suoi parametri, dichiarando la convenzione di layout. Verificala con il gradient check di questa lezione e riporta l'errore relativo massimo. Se fallisce su qualche parametro, cerca il punto in cui la funzione non è derivabile: un max, un arrotondamento, una soglia.

I cinque livelli salgono di difficoltà per il tipo di competenza richiesta: i primi due si risolvono con carta e penna e consolidano le identità, il terzo e il quarto richiedono di implementare e di confrontare due modi di derivare, l'ultimo porta il metodo su un caso del tuo lavoro. Sono da fare in ordine, perché ciascuno usa il risultato del precedente come oracolo.

Gli esercizi di codice di questa lezione sono in un solo file, esercizi.py (circa 20 KB): le cinque identità verificate con le differenze centrate, il conteggio delle operazioni nei due modi, l'autograd da completare e il gradient check su una rete intera. Richiede Python 3.12 o successivo e NumPy 2.3. Eseguilo con python3 esercizi.py: gira anche con le parti da completare lasciate vuote e ti dice quali mancano.

Recap

Mappa mentale. Dichiara il layout: ∂f/∂x ha la forma di x. Cinque identità: ∂(aᵀx)/∂x = a; ∂(xᵀAx)/∂x = (A + Aᵀ)x; ∂‖Ax − b‖²/∂x = 2Aᵀ(Ax − b); ∂L/∂W = δxᵀ; softmax più entropia incrociata dà ŷ − y. Jacobiano (m, n), catena come prodotto di Jacobiani, ma in pratica solo VJP. Hessiano simmetrico, autovalori per minimo/massimo/sella, Newton x ← x − H⁻¹∇f fuori scala a O(n³). Diretto costa n passate, inverso m: in ML m = 1. Autograd: dato, grad, genitori, chiusura; ordinamento topologico rovesciato; seme 1. Gradient check: differenze centrate, h = ∛ε ≈ 1e-5, errore relativo.

Cheat sheet.

3 concetti: il gradiente ha la forma dell'oggetto rispetto a cui deriva; il VJP è l'unica cosa che si calcola davvero; il modo inverso vince perché l'uscita è una sola. 3 errori: il gradiente trasposto su una matrice quadrata; h troppo piccolo nel gradient check; gradienti non azzerati fra due passi. 3 best practice: dichiarare la convenzione di layout; scrivere il gradient check prima del ciclo di addestramento; verificare con due oracoli indipendenti. 1 esercizio: il mini-progetto di questa lezione.

Prossima lezione: l'ottimizzazione con i vincoli, cioè da dove escono i moltiplicatori che la SVM usa senza spiegarli.

Domande per andare avanti

• Come si massimizza una funzione quando le variabili non sono libere ma devono soddisfare un vincolo, come ‖u‖ = 1 nella PCA? (Lezione 11 del modulo 02: moltiplicatori di Lagrange, dualità e condizioni KKT.)
• Perché la perdita di classificazione è proprio l'entropia incrociata, e non una qualunque misura di distanza fra distribuzioni? (Lezione 12 del modulo 02: entropia, entropia incrociata e divergenza KL.)
• Come si riduce la memoria del nastro quando la rete è troppo profonda per starci? (Ricalcolo dei valori intermedi, gradient checkpointing: si baratta tempo con memoria, nella direzione opposta a quella di questa lezione.)

Le tre domande sono i ganci verso le lezioni successive: le prime due aprono la lezione 11 e la lezione 12 del modulo 02, la terza è l'altra faccia del compromesso tempo-memoria incontrato nel confronto fra modo diretto e modo inverso. Non servono per ripassare: servono per decidere che cosa leggere dopo.

Nella stessa serie trovi anche: 26 - Ottimizzazione vincolata: moltiplicatori di Lagrange, dualità e condizioni KKT (la lezione successiva).