"""Corso Machine Learning Engineer — Modulo 02, Lezione 10 Calcolo matriciale e differenziazione automatica. Quattro esercizi di codice. Completa le funzioni marcate con TODO e rilancia: lo script gira anche se non hai completato niente e ti dice che cosa manca. python3 esercizi.py Richiede Python 3.12+ e NumPy 2.3 (testato con NumPy 2.3.5). Nessun'altra libreria: l'autograd lo scrivi tu. Le soluzioni commentate sono in fondo al file. CONVENZIONE DI LAYOUT (dichiararla e` meta` del lavoro): - il gradiente di uno scalare ha la stessa forma dell'oggetto rispetto a cui si deriva, quindi d(f)/d(W) ha la forma di W; - il Jacobiano di F: R^n -> R^m ha forma (m, n), una riga per uscita. """ import numpy as np TODO = "__TODO__" # ====================================================================== # Strumenti già pronti: non c'è niente da completare qui. # ====================================================================== 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. Vale 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 def errore_relativo(a, b, eta=1e-12): """Criterio di accettazione del gradient check: mai l'errore assoluto.""" a, b = np.ravel(np.asarray(a, dtype=float)), np.ravel(np.asarray(b, dtype=float)) scala = np.maximum(np.maximum(np.abs(a), np.abs(b)), eta) return float(np.max(np.abs(a - b) / scala)) # ====================================================================== # ESERCIZIO 1 — le cinque identità, verificate contro le differenze centrate # ====================================================================== def grad_prodotto_scalare(a, x): """d(a^T x)/dx. Suggerimento: una riga, e non dipende da x.""" return TODO def grad_forma_quadratica(A, x): """d(x^T A x)/dx, con A NON necessariamente simmetrica. Trappola: 2Ax vale solo se A == A.T. Qui serve la formula generale. """ return TODO def grad_minimi_quadrati(A, b, x): """d(||Ax - b||^2)/dx. Annullarla da' l'equazione normale A^T A x = A^T b.""" return TODO def grad_pesi_strato_lineare(delta, x): """d(L)/dW per z = Wx, dato delta = d(L)/dz. È un prodotto esterno. Controlla la forma PRIMA del valore: deve essere (len(delta), len(x)). """ return TODO def softmax(z): """Softmax stabile: sottrarre il massimo non cambia il risultato matematico.""" e = np.exp(z - z.max()) return e / e.sum() def grad_softmax_entropia_incrociata(z, y_onehot): """d(L)/dz per L = -sum(y * log(softmax(z))). Il risultato è sorprendentemente corto.""" return TODO def esegui_esercizio_1(): print("ESERCIZIO 1 — le cinque identità") a = np.array([1.0, 2.0, 3.0]) x = np.array([4.0, 5.0, 6.0]) A_ns = np.array([[2.0, 1.0], [0.0, 3.0]]) x2 = np.array([1.0, 2.0]) A_ls = np.array([[1.0, 0.0], [1.0, 1.0], [0.0, 2.0]]) b_ls = np.array([1.0, 1.0, 1.0]) W = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) x3 = np.array([1.0, 0.0, -1.0]) delta = np.array([10.0, 20.0]) z = np.array([1.0, 2.0, 0.0]) y = np.array([0.0, 1.0, 0.0]) prove = [ ("d(a^T x)/dx", grad_prodotto_scalare(a, x), gradiente_numerico(lambda v: a @ v, x)), ("d(x^T A x)/dx, A non simmetrica", grad_forma_quadratica(A_ns, x2), gradiente_numerico(lambda v: v @ A_ns @ v, x2)), ("d(||Ax - b||^2)/dx", grad_minimi_quadrati(A_ls, b_ls, x2), gradiente_numerico(lambda v: float((A_ls @ v - b_ls) @ (A_ls @ v - b_ls)), x2)), ("d(L)/dW con z = Wx", grad_pesi_strato_lineare(delta, x3), gradiente_numerico(lambda M: delta @ (M @ x3), W)), ("softmax + entropia incrociata", grad_softmax_entropia_incrociata(z, y), gradiente_numerico(lambda v: float(-np.log(softmax(v) @ y)), z)), ] for nome, tuo, numerico in prove: if tuo is TODO: print(f" [ da fare ] {nome}") continue tuo = np.asarray(tuo, dtype=float) if tuo.shape != numerico.shape: print(f" [ FORMA ] {nome}: tu {tuo.shape}, atteso {numerico.shape}") continue err = errore_relativo(tuo, numerico) esito = "ok" if err < 1e-6 else "NO" print(f" [ {esito:^8} ] {nome}: errore relativo {err:.2e} | {riga(tuo)}") print() # ====================================================================== # ESERCIZIO 2 — contare le operazioni: modo diretto contro modo inverso # ====================================================================== class Duale: """Numero duale a + b*eps con eps^2 = 0. Ogni operazione conta 2: valore e derivata.""" def __init__(self, valore, derivata, conta): self.valore, self.derivata, self.conta = valore, derivata, conta def __add__(self, altro): self.conta[0] += 2 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 gradiente_modo_diretto(valori, conta): """Gradiente COMPLETO di f(x) = sum(x_i^2) con i numeri duali. Una passata da' una sola derivata parziale: servono len(valori) passate. Restituisci un array NumPy con tutte le derivate. """ return TODO def gradiente_modo_inverso(valori, conta): """Lo stesso gradiente con UNA passata avanti e UNA all'indietro. Conta 1 per ogni operazione elementare in avanti e 1 per ogni operazione all'indietro. """ return TODO def esegui_esercizio_2(): print("ESERCIZIO 2 — il costo dei due modi") rng = np.random.default_rng(7) print(f" {'n':>4} | {'diretto':>9} | {'inverso':>9} | {'rapporto':>9} | esito") for n in (2, 5, 10, 20, 50, 100): valori = rng.normal(size=n) c_d, c_i = [0], [0] g_d = gradiente_modo_diretto(valori, c_d) g_i = gradiente_modo_inverso(valori, c_i) if g_d is TODO or g_i is TODO: print(f" {n:>4} | [ da fare: completa gradiente_modo_diretto / _inverso ]") continue atteso = 2.0 * valori ok = (errore_relativo(g_d, atteso) < 1e-12 and errore_relativo(g_i, atteso) < 1e-12) rapporto = c_d[0] / c_i[0] if c_i[0] else float("nan") print(f" {n:>4} | {c_d[0]:>9} | {c_i[0]:>9} | {rapporto:>9.1f} | {'ok' if ok else 'NO'}") print() # ====================================================================== # ESERCIZIO 3 — l'autograd: completa le tre regole mancanti # ====================================================================== def somma_su_broadcast(grad, forma): """Rimette il gradiente nella forma dell'operando sommando gli assi del 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 def __add__(self, altro): altro = altro if isinstance(altro, Nodo) else Nodo(altro) fuori = Nodo(self.dato + altro.dato, (self, altro), "+") def propaga(): 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(): 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(): # TODO (3a): due righe. Il gradiente sull'ingresso è delta @ W^T, # quello sui pesi è X^T @ delta. Ricordati di usare +=. raise NotImplementedError("TODO 3a: completa la regola all'indietro di @") fuori._indietro = propaga return fuori def tanh(self): t = np.tanh(self.dato) fuori = Nodo(t, (self,), "tanh") def propaga(): # TODO (3b): una riga. La derivata di tanh è 1 - tanh^2, e t ce l'hai già. raise NotImplementedError("TODO 3b: completa la regola all'indietro di tanh") fuori._indietro = propaga return fuori def relu(self): fuori = Nodo(np.maximum(self.dato, 0.0), (self,), "relu") def propaga(): self.grad += (self.dato > 0) * fuori.grad 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, seme a 1 sulla radice, grafo percorso all'inverso.""" 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) visita(self) self.grad = np.ones_like(self.dato) for nodo in reversed(ordine): nodo._indietro() def entropia_incrociata(logit, y_onehot): """Softmax ed entropia incrociata fuse in un nodo solo. La derivata è (p - y)/m.""" 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 def esegui_esercizio_3(): print("ESERCIZIO 3 — l'autograd su un caso verificabile a mano") X = Nodo([[1.0, -2.0, 3.0]]) W = Nodo([[1.0, 2.0], [0.0, 1.0], [1.0, 0.0]]) try: L = (X @ W).relu().sum() L.indietro() except NotImplementedError as mancante: print(f" [ da fare ] {mancante}") print() return atteso_W = np.array([[1.0, 0.0], [-2.0, 0.0], [3.0, 0.0]]) atteso_X = np.array([[1.0, 0.0, 1.0]]) print(f" L = {float(L.dato):.4f} (atteso 4.0000)") print(f" [ {'ok' if errore_relativo(W.grad, atteso_W) < 1e-12 else 'NO':^8} ] dL/dW: {riga(W.grad)}") print(f" [ {'ok' if errore_relativo(X.grad, atteso_X) < 1e-12 else 'NO':^8} ] dL/dX: {riga(X.grad)}") a = Nodo(3.0) b = a * a + a b.indietro() print(f" [ {'ok' if abs(float(a.grad) - 7.0) < 1e-12 else 'NO':^8} ] accumulo: d(a^2+a)/da = {float(a.grad):.4f} (atteso 7)") print() # ====================================================================== # ESERCIZIO 4 — gradient check su una rete intera, sulla spirale a due bracci # ====================================================================== def spirale_a_due_bracci(per_braccio=300, rumore=0.18, seme=42): """Due bracci intrecciati. Dati SIMULATI, non un dataset reale.""" rng = 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 px = t * np.cos(angolo) + rng.normal(0.0, rumore, per_braccio) py = t * np.sin(angolo) + rng.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) def gradient_check(perdita_di, parametri, quanti=20, h=1e-5, seme=1): """Campiona `quanti` parametri a caso e confronta con le differenze centrate. Restituisce (errore relativo peggiore, indice globale del parametro peggiore). Su una rete vera non si controllano MAI tutti i parametri: bastano venti. """ rng = np.random.default_rng(seme) piatti = [(p, i) for p in parametri for i in range(p.dato.size)] scelti = rng.choice(len(piatti), size=min(quanti, len(piatti)), replace=False) peggiore, colpevole = 0.0, -1 for indice in scelti: p, i = piatti[int(indice)] 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 numerico = (piu - meno) / (2 * h) err = errore_relativo(numerico, float(p.grad.flat[i])) if err > peggiore: peggiore, colpevole = err, int(indice) return peggiore, colpevole def addestra_spirale(epoche=3000, passo=0.5, seme=0): """Addestra una rete 2-16-16-2 con il SOLO autograd di questo file. TODO (4): completa il ciclo seguendo i tre passi commentati e restituisci (perdita_finale, accuratezza_finale, errore_del_gradient_check). """ return TODO def esegui_esercizio_4(): print("ESERCIZIO 4 — addestramento e gradient check sulla spirale") risultato = addestra_spirale() if risultato is TODO: print(" [ da fare ] completa addestra_spirale()") print() return perdita_finale, accuratezza, errore = risultato print(f" perdita finale : {perdita_finale:.4f}") print(f" accuratezza finale : {accuratezza:.4f}") print(f" [ {'ok' if errore < 1e-6 else 'NO':^8} ] gradient check su 20 parametri: {errore:.2e}") print() # ====================================================================== def main(): print("Lezione 10 — calcolo matriciale e differenziazione automatica") print(f"NumPy {np.__version__}\n") esegui_esercizio_1() esegui_esercizio_2() esegui_esercizio_3() esegui_esercizio_4() print("Le soluzioni commentate sono in fondo a questo file.") if __name__ == "__main__": main() # ====================================================================== # SOLUZIONI COMMENTATE # ====================================================================== # # --- ESERCIZIO 1 ------------------------------------------------------- # # def grad_prodotto_scalare(a, x): # return a # non dipende da x: la funzione è lineare # # def grad_forma_quadratica(A, x): # return (A + A.T) @ x # 2Ax SOLO se A è simmetrica # # def grad_minimi_quadrati(A, b, x): # return 2 * A.T @ (A @ x - b) # annullarla da' A^T A x = A^T b # # def grad_pesi_strato_lineare(delta, x): # return np.outer(delta, x) # prodotto esterno: forma (len(delta), len(x)) # # def grad_softmax_entropia_incrociata(z, y_onehot): # return softmax(z) - y_onehot # y_cappello meno y, niente Jacobiano K x K # # L'errore relativo atteso è dell'ordine di 1e-10 o meno su tutte e cinque: # le differenze centrate sono esatte sui polinomi di secondo grado e quasi # esatte sul resto, a patto di usare h = 1e-5 (vedi la sezione sul gradient check). # # --- ESERCIZIO 2 ------------------------------------------------------- # # def gradiente_modo_diretto(valori, conta): # g = np.zeros(len(valori)) # for i_seme in range(len(valori)): # 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 # g[i_seme] = s.derivata # return g # # def gradiente_modo_inverso(valori, conta): # s = 0.0 # for v in valori: # passata in avanti: 2 operazioni per ingresso # q = v * v; conta[0] += 1 # s = s + q; conta[0] += 1 # g = np.zeros(len(valori)) # aggiunto_s = 1.0 # il seme # for i, v in enumerate(valori): # passata all'indietro: altre 2 per ingresso # aggiunto_q = aggiunto_s; conta[0] += 1 # g[i] = aggiunto_q * (2 * v); conta[0] += 1 # return g # # Con questa contabilità il modo diretto costa 4n^2 e il modo inverso 4n: # il rapporto stampato è esattamente n, cioè 2, 5, 10, 20, 50, 100. # È tutta qui la ragione per cui le librerie usano il modo inverso: la # perdita è UNA uscita e i parametri sono milioni di ingressi. # # --- ESERCIZIO 3 ------------------------------------------------------- # # TODO 3a, dentro __matmul__: # self.grad += fuori.grad @ altro.dato.T # delta W^T # altro.grad += self.dato.T @ fuori.grad # X^T delta, il prodotto esterno # # TODO 3b, dentro tanh: # self.grad += (1.0 - t ** 2) * fuori.grad # # Attenzione a tre cose. (1) Il += non è uno stile: se un nodo compare due # volte nel grafo i contributi si sommano, e con un = ne perderesti uno. # (2) L'ordine dei fattori in @ non è libero: scrivere self.dato.T @ fuori.grad # al posto di fuori.grad @ altro.dato.T da' un errore di forma sugli strati # rettangolari e, sugli strati QUADRATI, nessun errore e un numero sbagliato. # (3) In t = 0 la ReLU non è derivabile e qui si sceglie 0: è una convenzione, # ed è il motivo per cui il gradient check fallisce su quel punto. # # Valori attesi dell'esercizio 3, verificabili a mano: # Z = X @ W = [4, 0], quindi delta = relu'(Z) = [1, 0] # dL/dW = X^T delta = [[1, 0], [-2, 0], [3, 0]] # dL/dX = delta W^T = [[1, 0, 1]] # d(a^2 + a)/da in a = 3 vale 2a + 1 = 7 # # --- ESERCIZIO 4 ------------------------------------------------------- # # def addestra_spirale(epoche=3000, passo=0.5, seme=0): # X_dati, y_dati = spirale_a_due_bracci() # Y = np.eye(2)[y_dati] # rng = np.random.default_rng(seme) # # def strato(n_in, n_out): # return (Nodo(rng.normal(0.0, np.sqrt(2.0 / n_in), (n_in, n_out))), # Nodo(np.zeros(n_out))) # # W1, b1 = strato(2, 16) # W2, b2 = strato(16, 16) # W3, b3 = strato(16, 2) # parametri = [W1, b1, W2, b2, W3, b3] # # def avanti(X): # H1 = (Nodo(X) @ W1 + b1).tanh() # H2 = (H1 @ W2 + b2).tanh() # return H2 @ W3 + b3 # # for _ in range(epoche): # perdita = entropia_incrociata(avanti(X_dati), Y) # for p in parametri: # (1) azzerare i gradienti # p.grad = np.zeros_like(p.dato) # perdita.indietro() # (2) una passata all'indietro # for p in parametri: # (3) un passo di discesa # p.dato -= passo * p.grad # # previste = avanti(X_dati).dato.argmax(axis=1) # accuratezza = float((previste == y_dati).mean()) # perdita_finale = float(entropia_incrociata(avanti(X_dati), Y).dato) # errore, _ = gradient_check( # lambda: float(entropia_incrociata(avanti(X_dati), Y).dato), parametri) # return perdita_finale, accuratezza, errore # # Nota sul gradient check finale: va fatto quando i gradienti in .grad sono # quelli dei pesi ATTUALI. Nel ciclo qui sopra l'ultima passata all'indietro # precede l'ultimo aggiornamento, quindi prima di controllare conviene rifare # una passata avanti e indietro sui pesi finali. # # Caso limite richiesto dal mini-progetto: se il parametro campionato ha # gradiente esattamente nullo, l'errore relativo ha 0 sia al numeratore sia # al denominatore; il termine eta = 1e-12 evita la divisione per zero e il # risultato è 0, cioè il controllo passa senza avere verificato nulla. # Per questo si campionano venti parametri e non uno.