22 - Analisi matematica 2: gradiente, integrali, discesa del gradiente

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

Obiettivo. Estendere la derivata alle funzioni di molte variabili (derivate parziali e gradiente), dimostrare che il gradiente indica la direzione di salita più ripida, trasformare questa idea nell'algoritmo che addestra quasi ogni modello (la discesa del gradiente), vedere la retropropagazione come regola della catena applicata a un neurone, e introdurre gli integrali come aree. Tutto con codice che gira e che si può verificare.

Falso mito da sfatare. «Un passo più grande di discesa del gradiente accelera sempre l'addestramento.» C'è una soglia oltre la quale i passi divergono: la deriveremo e la misureremo.

Introduzione

Nella lezione precedente la funzione dipendeva da una sola variabile. Una rete neurale ha milioni di parametri: la sua «funzione di perdita» (l'errore da minimizzare) è una funzione di milioni di variabili. Servono tre strumenti: come derivare rispetto a una variabile alla volta, come combinare queste derivate in un vettore che indica dove andare, e come usare quel vettore per scendere passo dopo passo.

Termini che userai (definiti qui al primo uso): funzione di perdita (loss): numero che misura quanto il modello sbaglia, da rendere piccolo; parametri: i numeri che il modello regola durante l'addestramento; minimo: un punto in cui la funzione non scende ulteriormente.

Il percorso parte dalle derivate parziali e dal gradiente, dimostra che il gradiente indica la salita più ripida e ne ricava la discesa del gradiente, di cui misura il passo massimo e l'effetto delle valli strette. Applica poi l'algoritmo alla regressione lineare e alla retropropagazione su un neurone, tratta convessità e minimi locali, e chiude con gli integrali come aree. Dalla lezione Analisi 1 servono la derivata e la regola della catena; l'algoritmo completo è poi ripreso in La discesa del gradiente.

Funzioni di più variabili e derivate parziali

L'analogia: il paesaggio

Una funzione f(x, y) assegna a ogni punto di un piano un'altitudine: è un paesaggio. Se sei in un punto e ti muovi verso est, la quota sale o scende con una certa pendenza; se ti muovi verso nord ne trovi un'altra. La derivata parziale rispetto a x, ∂f/∂x, è la pendenza che trovi muovendoti solo lungo x, tenendo ferma y. Si calcola come la derivata di una variabile, trattando le altre come costanti.

Esempio: f(x, y) = x² + 3y² + xy. Allora ∂f/∂x = 2x + y (il termine 3y² è costante rispetto a x; xy ha derivata y) e ∂f/∂y = 6y + x.

Il gradiente è il vettore delle derivate parziali:

∇f(x, y) = (∂f/∂x, ∂f/∂y)

(∇ si legge «nabla»). Per una funzione di d variabili è un vettore di d componenti.

import numpy as np

f = lambda p: p[0] ** 2 + 3 * p[1] ** 2 + p[0] * p[1]
grad_f = lambda p: np.array([2 * p[0] + p[1], 6 * p[1] + p[0]])        # derivate parziali calcolate a mano

def grad_numerico(funzione, p, h=1e-5):
    """Gradiente per differenze centrali: una derivata parziale per coordinata."""
    p = np.asarray(p, dtype=float)
    g = np.zeros_like(p)
    for i in range(p.size):
        passo = np.zeros_like(p)
        passo[i] = h                                  # ci muoviamo solo lungo la coordinata i
        g[i] = (funzione(p + passo) - funzione(p - passo)) / (2 * h)
    return g

p0 = np.array([1.5, -0.7])
print("gradiente analitico:", grad_f(p0))
print("gradiente numerico :", np.round(grad_numerico(f, p0), 6))

Output ottenuto:

gradiente analitico: [ 2.3 -2.7]
gradiente numerico : [ 2.3 -2.7]

Cosa fa il codice. grad_f implementa le due formule ricavate sopra. grad_numerico calcola ogni componente per differenza centrale (lezione precedente), muovendosi di ±h solo lungo la coordinata i: è la definizione di derivata parziale, tradotta in codice. I due vettori coincidono: questo confronto è il gradient check, che ritroveremo.

Il gradiente indica la salita più ripida

Perché proprio il gradiente? Muoviamoci da p nella direzione di un vettore di lunghezza 1 (versore) u per un tratto piccolo t. Per lo sviluppo di Taylor al primo ordine (la linearizzazione della lezione precedente, ora in più variabili):

f(p + t·u) ≈ f(p) + t · (∇f(p) · u)

dove · è il prodotto scalare (la somma dei prodotti delle componenti). La quantità D_u f = ∇f · u è la derivata direzionale: di quanto cresce f per unità di cammino nella direzione u.

Quale u la rende massima? Il prodotto scalare di due vettori vale |∇f| · |u| · cos θ, con θ l'angolo tra loro. Poiché |u| = 1:

D_u f = |∇f| · cos θ ≤ |∇f|,

con uguaglianza solo per cos θ = 1, cioè θ = 0: u parallelo al gradiente. Quindi:

• il gradiente punta nella direzione di massima crescita, e la sua lunghezza è la pendenza massima;
• −∇f punta nella direzione di massima discesa;
• nelle direzioni perpendicolari al gradiente (cos θ = 0) la funzione, al primo ordine, non cambia: sono le curve di livello (le linee di quota costante di una mappa).

Le tre conseguenze sono la stessa disuguaglianza letta in tre direzioni: salita massima lungo il gradiente, discesa massima lungo il suo opposto, nessuna variazione al primo ordine lungo le curve di livello. Per l'ottimizzazione conta la seconda, perché dice come scegliere il passo successivo conoscendo solo la pendenza locale; il listato seguente la verifica numericamente provando 3601 direzioni.

angoli = np.linspace(0, 2 * np.pi, 3601)
versori = np.stack([np.cos(angoli), np.sin(angoli)], axis=1)         # 3601 direzioni di lunghezza 1
g = grad_f(p0)
derivate_direzionali = versori @ g                                    # D_u f = ∇f · u per ogni u

migliore = versori[np.argmax(derivate_direzionali)]
print("direzione di massima crescita:", np.round(migliore, 4))
print("gradiente normalizzato       :", np.round(g / np.linalg.norm(g), 4))
print("pendenza massima:", round(float(derivate_direzionali.max()), 4), "| |∇f| =", round(float(np.linalg.norm(g)), 4))

Output ottenuto:

direzione di massima crescita: [ 0.6481 -0.7615]
gradiente normalizzato       : [ 0.6485 -0.7612]
pendenza massima: 3.5468 | |∇f| = 3.5468

Cosa fa il codice. Genera 3601 direzioni di lunghezza 1 (le righe di versori), calcola per ciascuna la derivata direzionale con un solo prodotto matrice-vettore (@), e cerca la direzione migliore. Coincide con il gradiente normalizzato e la pendenza massima è |∇f|, come dimostrato.

La discesa del gradiente

L'analogia: l'escursionista nella nebbia

Sei su una montagna, nella nebbia fitta, e vuoi scendere a valle. Non vedi il paesaggio: senti solo la pendenza sotto i piedi. Strategia: ad ogni passo guarda in quale direzione il terreno scende più ripido e fai un passo in quella direzione. Ripeti. È la discesa del gradiente (Cauchy la propose nel 1847):

xₖ₊₁ = xₖ − η · ∇f(xₖ)

dove η («eta», il tasso di apprendimento, learning rate) è la lunghezza del passo in proporzione alla pendenza. Il segno meno perché si scende.

Quanto grande dev'essere il passo? Un caso risolto fino in fondo

Prendiamo la funzione più semplice con un minimo, f(x) = (x − 3)², con minimo in x = 3. La derivata è f'(x) = 2(x − 3). Un passo dà

xₖ₊₁ − 3 = xₖ − 3 − η·2(xₖ − 3) = (1 − 2η)(xₖ − 3).

La distanza dal minimo viene moltiplicata per (1 − 2η) a ogni passo. Quindi:

• se |1 − 2η| < 1, cioè 0 < η < 1, la distanza tende a zero: converge;
• se η = 0.5 il fattore vale 0: arriva al minimo in un solo passo;
• se 0.5 < η < 1 il fattore è negativo: si oscilla attorno al minimo, ma ci si avvicina;
• se η > 1 il fattore ha modulo > 1: la distanza cresce, i passi divergono.

I quattro regimi dipendono solo dal modulo e dal segno del fattore (1 − 2η): monotono se positivo, oscillante se negativo, divergente se il modulo supera 1. Il risultato generale è lo stesso per qualunque parabola, con 2 sostituito dalla curvatura, ed è la ragione per cui il passo massimo dipende dalla curvatura della funzione e non può essere scelto una volta per tutte.

def discesa_1d(eta, x0=-5.0, passi=8):
    x, traiettoria = x0, [x0]
    for _ in range(passi):
        x = x - eta * 2 * (x - 3)              # x_{k+1} = x_k - η f'(x_k)
        traiettoria.append(x)
    return traiettoria

for eta in (0.1, 0.5, 0.9, 1.05):
    t = discesa_1d(eta)
    print(f"η = {eta:<4}: fattore (1-2η) = {1 - 2 * eta:+.2f} | traiettoria: {[round(v, 2) for v in t[:7]]}")

Output ottenuto:

η = 0.1 : fattore (1-2η) = +0.80 | traiettoria: [-5.0, -3.4, -2.12, -1.1, -0.28, 0.38, 0.9]
η = 0.5 : fattore (1-2η) = +0.00 | traiettoria: [-5.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
η = 0.9 : fattore (1-2η) = -0.80 | traiettoria: [-5.0, 9.4, -2.12, 7.1, -0.28, 5.62, 0.9]
η = 1.05: fattore (1-2η) = -1.10 | traiettoria: [-5.0, 11.8, -6.68, 13.65, -8.71, 15.88, -11.17]

Cosa fa il codice. La funzione discesa_1d applica l'aggiornamento per 8 passi partendo da x = −5. Con η = 0.1 il fattore è 0.8: avvicinamento lento e monotono. Con η = 0.5 il fattore è 0: dopo il primo passo si è già in 3. Con η = 0.9 il fattore è −0.8: si oscilla sopra e sotto 3 avvicinandosi. Con η = 1.05 il fattore è −1.1: l'oscillazione si amplia, i valori divergono. Il numero stampato come «fattore» predice esattamente ciò che si vede.

Valli strette e il condizionamento

Nel mondo reale le variabili hanno curvature molto diverse. Prendiamo f(x, y) = x² + 10y²: una «valle» stretta e lunga, ripida lungo y (curvatura 20) e dolce lungo x (curvatura 2). Il gradiente è (2x, 20y), e un passo dà

x ← (1 − 2η)·x, y ← (1 − 20η)·y.

Per convergere servono entrambe |1 − 2η| < 1 e |1 − 20η| < 1: la seconda è la più restrittiva e impone η < 0.1. In generale, il passo massimo è 2 / (curvatura più grande). Ma con η vicino a quel limite la direzione dolce converge molto lentamente: il fattore (1 − 2η) resta vicino a 1. Il rapporto tra curvatura massima e minima (qui 10) è il numero di condizionamento: più è grande, più la discesa è lenta o a zig-zag.

def discesa_2d(eta, passi=200, p0=(-2.6, 0.8), tol=1e-6):
    p = np.array(p0)
    for k in range(passi):
        p = p - eta * np.array([2 * p[0], 20 * p[1]])
        if np.linalg.norm(p) < tol:
            return f"converge in {k + 1} passi"
        if np.linalg.norm(p) > 1e6:
            return f"DIVERGE dopo {k + 1} passi"
    return f"dopo {passi} passi |p| = {np.linalg.norm(p):.2e}"

for eta in (0.01, 0.05, 0.09, 0.099, 0.11):
    print(f"η = {eta:<5}: {discesa_2d(eta)}")

Output ottenuto:

η = 0.01 : dopo 200 passi |p| = 4.57e-02
η = 0.05 : converge in 141 passi
η = 0.09 : converge in 75 passi
η = 0.099: dopo 200 passi |p| = 1.41e-02
η = 0.11 : DIVERGE dopo 77 passi

Cosa fa il codice. Applica la discesa alla funzione x² + 10y² con vari passi e riporta se converge (norma sotto 10⁻⁶) o diverge (norma sopra 10⁶). Il passo ottimale per questo problema è un compromesso: η = 0.01 è sicuro ma dopo 200 passi non ha ancora raggiunto la tolleranza (direzione dolce lenta); con 0.05 e 0.09 converge in 141 e 75 passi, con zig-zag lungo y; a ridosso del limite (0.099) il fattore 1 − 20η ≈ −0.98 rende l'oscillazione quasi non smorzata e non si converge in 200 passi; 0.11 supera il limite 0.1 e diverge. È la stessa matematica che regola tutte le reti neurali, solo in milioni di dimensioni.

Conseguenza pratica: scalare le variabili. Se due caratteristiche hanno scale molto diverse (una in unità, l'altra in migliaia), la funzione di perdita ha una valle stretta come questa, e la discesa è lenta o instabile. Standardizzare le caratteristiche (lezione su NumPy) rende le curvature simili. È la ragione, più profonda, per cui nella lezione su scikit-learn abbiamo scalato prima di usare modelli lineari addestrati con questi metodi.

Regressione lineare con la discesa del gradiente

Riprendiamo la retta ai minimi quadrati della lezione di statistica, ma risolviamola senza la formula chiusa. La perdita è l'errore quadratico medio

L(a, b) = (1/n) Σ (yᵢ − a − b·xᵢ)².

Con la regola della catena (la derivata di un quadrato è 2·(...) per la derivata interna), indicando il residuo rᵢ = yᵢ − a − b·xᵢ:

∂L/∂a = −(2/n) Σ rᵢ, ∂L/∂b = −(2/n) Σ xᵢ·rᵢ.

rng = np.random.default_rng(0)
n = 200
x_grezza = rng.uniform(0, 100, n)                           # caratteristica con scala grande
y = 5.0 + 0.3 * x_grezza + rng.normal(0, 2.0, n)            # retta vera: a = 5, b = 0.3

def gd_regressione(x, y, eta, passi):
    a = b = 0.0
    for _ in range(passi):
        r = y - (a + b * x)
        a -= eta * (-2 * r.mean())
        b -= eta * (-2 * (x * r).mean())
    return a, b

# formula chiusa della lezione di statistica
b_chiusa = ((x_grezza - x_grezza.mean()) * (y - y.mean())).sum() / ((x_grezza - x_grezza.mean()) ** 2).sum()
a_chiusa = y.mean() - b_chiusa * x_grezza.mean()
print(f"formula chiusa          : a = {a_chiusa:.4f}, b = {b_chiusa:.4f}")

x_std = (x_grezza - x_grezza.mean()) / x_grezza.std()
a_s, b_s = gd_regressione(x_std, y, eta=0.1, passi=200)
# si riconverte alla scala originale: y = a_s + b_s*(x - m)/s  =>  b = b_s/s, a = a_s - b_s*m/s
print(f"GD su x standardizzata  : a = {a_s - b_s * x_grezza.mean() / x_grezza.std():.4f}, b = {b_s / x_grezza.std():.4f}  (200 passi, η = 0.1)")

for eta in (1e-3, 1e-4):
    a_g, b_g = gd_regressione(x_grezza, y, eta=eta, passi=200)
    esito = f"a = {a_g:.4f}, b = {b_g:.4f}" if abs(a_g) < 1e6 else f"DIVERGE (a = {a_g:.1e})"
    print(f"GD su x grezza, η = {eta:g}: {esito}")

Output ottenuto:

formula chiusa          : a = 5.0942, b = 0.2949
GD su x standardizzata  : a = 5.0942, b = 0.2949  (200 passi, η = 0.1)
GD su x grezza, η = 0.001: DIVERGE (a = -1.5e+162)
GD su x grezza, η = 0.0001: a = 0.0533, b = 0.3662

Cosa fa il codice. gd_regressione calcola ad ogni passo i residui e aggiorna a e b con le due derivate parziali (le medie corrispondono ai 1/n Σ). Confronta tre risultati con la formula chiusa: la discesa su x standardizzata converge in 200 passi con η = 0.1 e, riconvertita alla scala originale, coincide con la formula. Sulla variabile grezza (valori fino a 100) il passo deve essere piccolo: con η = 10⁻³ la discesa diverge (il codice lo segnala), con 10⁻⁴ si avvicina ma dopo 200 passi è ancora lontana dalla soluzione, in particolare sull'intercetta. Il motivo è il condizionamento di prima.

La versione «a pezzetti»: la discesa stocastica

Con milioni di esempi calcolare il gradiente su tutti i dati a ogni passo è troppo costoso. Si usa invece un mini-batch: un piccolo sottoinsieme casuale (diciamo 32 esempi) per stimare il gradiente. È una stima rumorosa ma molto più economica, e in media punta nella direzione giusta (la media dei gradienti su sottoinsiemi casuali è il gradiente completo). È la discesa stocastica del gradiente (SGD), legata all'idea di approssimazione stocastica di Robbins e Monro (1951).

def sgd_regressione(x, y, eta, epoche, batch, seme=1):
    rng = np.random.default_rng(seme)
    a = b = 0.0
    for _ in range(epoche):
        ordine = rng.permutation(len(x))                  # una «epoca» = un passaggio su tutti i dati, in ordine casuale
        for inizio in range(0, len(x), batch):
            idx = ordine[inizio:inizio + batch]
            r = y[idx] - (a + b * x[idx])
            a -= eta * (-2 * r.mean())
            b -= eta * (-2 * (x[idx] * r).mean())
    return a, b

a_sgd, b_sgd = sgd_regressione(x_std, y, eta=0.05, epoche=30, batch=16)
print("SGD (batch 16, 30 epoche):", round(a_sgd - b_sgd * x_grezza.mean() / x_grezza.std(), 3), round(b_sgd / x_grezza.std(), 4))

Output ottenuto:

SGD (batch 16, 30 epoche): 4.931 0.2985

Cosa fa il codice. Ad ogni epoca (un passaggio completo sui dati) mescola l'ordine degli esempi e li prende a gruppi di 16; per ogni gruppo calcola il gradiente sul solo gruppo e fa un passo. Il risultato si avvicina molto alla formula chiusa, ma con un po' di rumore: è il prezzo dell'economia. I valori esatti dipendono da seme e dimensione del batch.

La retropropagazione: la regola della catena su un neurone

Un neurone di classificazione calcola z = w·x + b, poi p = σ(z) (la sigmoide, lezione precedente), che si interpreta come probabilità della classe 1. L'errore per un esempio con etichetta y ∈ {0, 1} è la log-perdita (cross-entropy):

L = −[ y·ln p + (1 − y)·ln(1 − p) ].

Per addestrare servono ∂L/∂w e ∂L/∂b. La catena è w → z → p → L, e per la regola della catena le derivate si moltiplicano lungo il percorso:

∂L/∂w = (∂L/∂p) · (∂p/∂z) · (∂z/∂w).

I tre fattori: ∂L/∂p = −y/p + (1 − y)/(1 − p) = (p − y)/(p(1 − p)); ∂p/∂z = p(1 − p) (derivata della sigmoide); ∂z/∂w = x. Il prodotto semplifica il fattore p(1 − p):

∂L/∂z = p − y, ∂L/∂w = (p − y)·x, ∂L/∂b = p − y.

Un risultato elegante: il gradiente rispetto all'ingresso del neurone è semplicemente «previsto meno vero». La retropropagazione (backpropagation, Rumelhart, Hinton e Williams 1986) è questo: si calcolano i valori in avanti (passaggio forward), poi si moltiplicano le derivate all'indietro lungo il grafo, riutilizzando i valori già calcolati.

rng = np.random.default_rng(3)
# dati sintetici: due gruppi di punti nel piano
n = 300
X0 = rng.normal([-1.0, -1.0], 1.0, size=(n // 2, 2))
X1 = rng.normal([1.2, 1.0], 1.0, size=(n // 2, 2))
X = np.vstack([X0, X1])
y = np.concatenate([np.zeros(n // 2), np.ones(n // 2)])

sigmoide = lambda z: 1 / (1 + np.exp(-z))

def perdita(w, b):
    p = np.clip(sigmoide(X @ w + b), 1e-12, 1 - 1e-12)
    return float(-(y * np.log(p) + (1 - y) * np.log(1 - p)).mean())

def gradiente(w, b):
    p = sigmoide(X @ w + b)
    err = p - y                                   # ∂L/∂z per ogni esempio
    return X.T @ err / len(y), err.mean()         # medie dei (p - y)·x  e  dei (p - y)

# gradient check: confronto con le differenze finite
w0, b0 = np.array([0.3, -0.2]), 0.1
gw, gb = gradiente(w0, b0)
h = 1e-6
num_w = np.array([(perdita(w0 + h * e, b0) - perdita(w0 - h * e, b0)) / (2 * h) for e in np.eye(2)])
num_b = (perdita(w0, b0 + h) - perdita(w0, b0 - h)) / (2 * h)
print("gradiente retropropagato :", np.round(gw, 6), round(float(gb), 6))
print("gradiente numerico       :", np.round(num_w, 6), round(float(num_b), 6))

# addestramento con la discesa del gradiente
w, b = np.zeros(2), 0.0
for passo in range(401):
    gw, gb = gradiente(w, b)
    w, b = w - 0.5 * gw, b - 0.5 * gb
    if passo in (0, 10, 100, 400):
        accuratezza = ((sigmoide(X @ w + b) > 0.5) == y).mean()
        print(f"passo {passo:>3}: perdita {perdita(w, b):.4f} | accuratezza {accuratezza:.3f}")
print("pesi appresi:", np.round(w, 3), "| bias:", round(float(b), 3))

Output ottenuto:

gradiente retropropagato : [-0.428585 -0.549415] 0.025666
gradiente numerico       : [-0.428585 -0.549415] 0.025666
passo   0: perdita 0.4656 | accuratezza 0.917
passo  10: perdita 0.2291 | accuratezza 0.923
passo 100: perdita 0.1845 | accuratezza 0.937
passo 400: perdita 0.1834 | accuratezza 0.937
pesi appresi: [2.099 1.77 ] | bias: -0.457

Cosa fa il codice. perdita calcola la log-perdita media sull'intero insieme (con np.clip per evitare ln(0): lezione sui float). gradiente implementa esattamente le formule ricavate sopra, in forma vettoriale: X.T @ err / n è la media dei vettori (p − y)·x. Il gradient check confronta il gradiente retropropagato con quello numerico: le due righe coincidono fino alla sesta cifra, quindi la derivata è giusta. Poi il ciclo applica la discesa (η = 0.5) e stampa perdita e accuratezza in quattro momenti: la perdita scende e l'accuratezza sale verso il valore ottimo per due gruppi sovrapposti. Il modello è, in pratica, una regressione logistica: lo stesso algoritmo di LogisticRegression di scikit-learn, scritto a mano.

Convessità e minimi locali

Una funzione è convessa se il segmento che congiunge due punti qualsiasi del suo grafico sta sopra il grafico (una coppa). Per una funzione convessa ogni minimo locale è un minimo globale e la discesa del gradiente, con un passo adeguato, lo trova. L'errore quadratico della regressione lineare e la log-perdita della regressione logistica sono convessi (lo studieresti con la matrice delle derivate seconde; qui lo prendiamo come fatto: si veda Boyd e Vandenberghe).

Le reti neurali profonde non sono convesse: ci sono più valli. A seconda del punto di partenza si finisce in minimi diversi. Mostriamolo con una funzione a una variabile con due minimi, f(x) = x⁴ − 3x² + x.

f4 = lambda x: x ** 4 - 3 * x ** 2 + x
df4 = lambda x: 4 * x ** 3 - 6 * x + 1

for partenza in (-2.0, -0.2, 0.2, 2.0):
    x = partenza
    for _ in range(2000):
        x -= 0.01 * df4(x)
    print(f"partenza {partenza:>5}: arriva in x = {x:.4f}, f(x) = {f4(x):.4f}")

Output ottenuto:

partenza  -2.0: arriva in x = -1.3008, f(x) = -3.5139
partenza  -0.2: arriva in x = -1.3008, f(x) = -3.5139
partenza   0.2: arriva in x = 1.1309, f(x) = -1.0702
partenza   2.0: arriva in x = 1.1309, f(x) = -1.0702

Cosa fa il codice. Applica la discesa da quattro punti di partenza. Le partenze a sinistra di circa 0.17 finiscono nel minimo profondo (a sinistra, f ≈ −3.5); quelle a destra finiscono nel minimo più alto (a destra, f ≈ −1.1). Il minimo trovato dipende da dove inizi. Per questo l'inizializzazione dei parametri e l'uso di più partenze o di rumore (come quello dell'SGD) sono importanti.

Gli integrali come aree

L'analogia: percorso dalla velocità

Se conosci la velocità dell'auto in ogni istante, come ricostruisci quanta strada ha fatto? Dividi il tempo in intervalli piccolissimi Δt, in ciascuno la velocità è quasi costante: la strada percorsa è circa velocità × Δt. Sommi tutti i pezzi. È l'idea dell'integrale: l'area sotto una curva f tra a e b si approssima con n rettangoli sottili di base Δx = (b − a)/n e altezza f(xᵢ):

∫ₐᵇ f(x) dx = lim_{n→∞} Σ f(xᵢ)·Δx (somma di Riemann).

Il teorema fondamentale

Derivata e integrale sono operazioni inverse. Se F' = f, allora

∫ₐᵇ f(x) dx = F(b) − F(a).

L'intuizione è nella somma: dividi [a, b] in piccoli intervalli con estremi x₀ = a, x₁, …, xₙ = b. Per ogni intervallo, F(xᵢ₊₁) − F(xᵢ) ≈ F'(xᵢ)·Δx = f(xᵢ)·Δx (definizione di derivata). Sommando su tutti gli intervalli, i termini intermedi si elidono (F(x₁) compare con il segno meno e poi più...): resta F(b) − F(a), la somma telescopica. Quindi l'area è la «variazione totale» di F.

from scipy import integrate
import math

f = lambda x: x ** 2
a, b = 0.0, 3.0                                  # ∫ x² dx da 0 a 3 = 3³/3 = 9
print("valore esatto F(b) - F(a) = 9")
for n in (4, 16, 64, 256, 1024):
    xs = np.linspace(a, b, n, endpoint=False)    # estremi sinistri dei sottointervalli
    riemann = (f(xs) * (b - a) / n).sum()
    trapezi = np.trapezoid(f(np.linspace(a, b, n + 1)), dx=(b - a) / n)
    print(f"n = {n:>4}: somma di Riemann {riemann:.5f} (errore {abs(riemann - 9):.2e}) | trapezi {trapezi:.5f} (errore {abs(trapezi - 9):.2e})")

valore, errore = integrate.quad(f, a, b)
print("scipy.integrate.quad:", valore, "(stima dell'errore:", f"{errore:.1e})")

Output ottenuto:

valore esatto F(b) - F(a) = 9
n =    4: somma di Riemann 5.90625 (errore 3.09e+00) | trapezi 9.28125 (errore 2.81e-01)
n =   16: somma di Riemann 8.17383 (errore 8.26e-01) | trapezi 9.01758 (errore 1.76e-02)
n =   64: somma di Riemann 8.79016 (errore 2.10e-01) | trapezi 9.00110 (errore 1.10e-03)
n =  256: somma di Riemann 8.94733 (errore 5.27e-02) | trapezi 9.00007 (errore 6.87e-05)
n = 1024: somma di Riemann 8.98682 (errore 1.32e-02) | trapezi 9.00000 (errore 4.29e-06)
scipy.integrate.quad: 9.000000000000002 (stima dell'errore: 1.0e-13)

Cosa fa il codice. Approssima l'area sotto x² con rettangoli (somma di Riemann, usando l'estremo sinistro) e con trapezi (che «tagliano» meglio la curva), per n crescente. L'errore della somma di Riemann scende come 1/n, quello dei trapezi come 1/n²: i trapezi sono molto più efficienti. scipy.integrate.quad usa un algoritmo adattivo e restituisce anche una stima dell'errore.

Integrali e probabilità

La probabilità, per una variabile continua, è un'area (lezione sulle variabili aleatorie): P(a ≤ X ≤ b) = ∫ₐᵇ f(x) dx, e il valore atteso è E[X] = ∫ x·f(x) dx.

densita_normale = lambda x: math.exp(-x * x / 2) / math.sqrt(2 * math.pi)
area, _ = integrate.quad(densita_normale, -1.96, 1.96)
print("∫ densità normale standard tra -1.96 e 1.96 =", round(area, 6))
print("area totale su (-∞, +∞)                    =", round(integrate.quad(densita_normale, -np.inf, np.inf)[0], 6))

lam = 0.5
valore_atteso, _ = integrate.quad(lambda x: x * lam * math.exp(-lam * x), 0, np.inf)
print("E[X] per l'esponenziale con λ = 0.5 (teoria 1/λ = 2):", round(valore_atteso, 6))

Output ottenuto:

∫ densità normale standard tra -1.96 e 1.96 = 0.950004
area totale su (-∞, +∞)                    = 1.0
E[X] per l'esponenziale con λ = 0.5 (teoria 1/λ = 2): 2.0

Cosa fa il codice. Calcola come integrale ciò che nella lezione 4 abbiamo ricavato con la funzione erf: l'area della densità normale tra ±1.96 vale 0.95, e l'area totale 1. Poi verifica che il valore atteso dell'esponenziale, definito come integrale ∫ x·λe^(−λx) dx, vale 1/λ = 2. np.inf indica un estremo infinito.

Errori tipici (sintomo, causa, rimedio)

• Perdita che sale o nan. Causa: η troppo grande, o scala delle variabili. Rimedio: ridurre η (di un fattore 10), standardizzare le caratteristiche.
• Perdita che scende lentissimo. Causa: η troppo piccolo o valle stretta. Rimedio: aumentare η fin quasi al limite, ottimizzatori adattivi, scalare.
• Gradiente sbagliato senza accorgersene. Causa: una derivata scritta male. Rimedio: gradient check con differenze finite su un esempio piccolo.
• Dimenticare il segno meno (si sale invece di scendere). Rimedio: controllare che la perdita cali nei primi passi.
• Confondere una funzione non convessa con una convessa. Rimedio: più partenze; guardare la perdita di validazione.
• Usare ln(0) nella log-perdita. Rimedio: np.clip o funzioni stabili.
• Integrare numericamente con troppo pochi punti o con singolarità. Rimedio: aumentare n, usare quad e controllare la stima dell'errore.
• Dimenticare di dividere per n nella media dei gradienti: cambia la scala del passo ottimale.

Le prime due voci sono i due volti dello stesso parametro, η, e si distinguono guardando la curva della perdita: sale o dà nan nel primo caso, scende appena nel secondo. La terza, il gradiente scritto male, è la più insidiosa perché l'algoritmo gira senza errori e produce soltanto risultati peggiori: per questo il gradient check va fatto prima di ogni addestramento.

Come lavorare nella pratica

Parti dal gradient check ogni volta che scrivi un gradiente a mano, registra la perdita a ogni epoca e disegnala (la forma dice se η è adatto), standardizza le variabili prima di addestrare e fissa il seme per rendere riproducibili mescolamento e inizializzazione. In produzione usa le librerie (scikit-learn, PyTorch): scrivere il gradiente a mano serve a capire, non a sostituirle.

I termini nuovi sono pochi. La derivata parziale è la pendenza rispetto a una variabile con le altre ferme, il gradiente ∇f è il vettore delle derivate parziali e la derivata direzionale misura la variazione per unità di cammino in una direzione data. Il tasso di apprendimento η è la lunghezza del passo; un'epoca è un passaggio completo sui dati e un mini-batch il sottoinsieme che stima il gradiente. Una funzione convessa ha un solo minimo; la retropropagazione calcola i gradienti con la regola della catena all'indietro lungo il grafo; la somma di Riemann approssima l'area sotto una curva con rettangoli.

Mini-progetto

1. Implementa la discesa del gradiente per la regressione logistica sul dataset del tumore al seno di scikit-learn (standardizzato), confronta accuratezza e coefficienti con LogisticRegression(penalty=None).
2. Per la funzione f(x, y) = x² + 10y² verifica per più valori di η che il limite di convergenza sia 0.1 e che con η = 0.1 esatto il comportamento sia oscillatorio.
3. Disegna (con Matplotlib) la perdita in funzione del passo per η = 0.01, 0.1, 0.5, 1.0 nella regressione logistica sintetica.

Il primo punto usa un dataset incluso in scikit-learn, quindi non richiede download; il secondo ripete per due variabili la verifica che la lezione ha fatto con una; il terzo rende visibili le curve di perdita dei quattro regimi (lento, buono, oscillante, divergente) che finora hai letto come numeri.

Challenge

• Easy. Calcola a mano il gradiente di f(x, y) = x²y + sin(y) e verificalo con grad_numerico.
• Medium. Dimostra che, per f(x) = (a/2)(x − m)², la discesa del gradiente converge se e solo se 0 < η < 2/a, e trova l'η che converge in un passo.
• Hard. Aggiungi alla discesa l'inerzia (momentum): v ← β·v + ∇f, x ← x − η·v, con β = 0.9. Mostra su x² + 10y² che riduce lo zig-zag e confronta il numero di passi.
• Expert. Estendi la retropropagazione a una rete con uno strato nascosto da 3 neuroni (sigmoide) e un'uscita; ricava a mano i gradienti con la regola della catena, implementali in NumPy e verificali con il gradient check.
• Enterprise. Scegli un problema reale di regressione della tua organizzazione. Addestralo con la discesa del gradiente scritta a mano, con e senza standardizzazione, e documenta tasso di apprendimento scelto, numero di epoche e perdita finale.

La sfida Medium è la generalizzazione algebrica del caso risolto nel testo, la Hard mostra il rimedio più semplice al zig-zag delle valli strette, la Expert è il passaggio logico dalla regressione logistica alle reti profonde. Se ne scegli una sola, scegli la Expert: prepara direttamente la lezione sulla retropropagazione completa.

Che cosa sai fare adesso

Sai calcolare un gradiente e verificarlo con le differenze finite, dimostrare perché indica la salita più ripida, prevedere se la discesa del gradiente converge dal fattore 1 − 2η (o dalla curvatura massima nelle valli strette) e addestrare a mano una regressione lineare e un neurone logistico. Sai anche leggere la retropropagazione come regola della catena e calcolare un'area con le somme di Riemann e con scipy.integrate.quad. Il passo successivo è l'algebra lineare, il linguaggio con cui X @ w sostituisce i cicli.

Nella stessa serie trovi anche: 19 - Variabili aleatorie e distribuzioni, 25 - Calcolo matriciale e differenziazione automatica: dal Jacobiano a un autograd scritto da zero, 26 - Ottimizzazione vincolata: moltiplicatori di Lagrange, dualità e condizioni KKT (che richiamano questa lezione).