La più semplice rete neurale, per trovare la retta che meglio modella un insieme di punti $(x,y)$¶
Luca Mari, aggiornamento maggio 2024
Quest'opera è distribuita con Licenza Creative Commons Attribuzione - Non commerciale - Condividi allo stesso modo 4.0 Internazionale.
[i file di questa attività: linear.ipynb, linearutils.py]
Obiettivi: comprendere il funzionamento di base di una (molto) semplice rete neurale, per come specificato in Python mediante un modulo di alto livello come PyTorch
, e poi anche comprendendo come la rete è addestrata e poi funziona "a basso livello".
Precompetenze: basi di Python; almeno qualche idea di analisi matematica.
Per eseguire questo notebook, supponiamo con VSCode, occorre:
- installare un interprete Python
- scaricare da https://code.visualstudio.com/download e installare VSCode
- eseguire VSCode e attivare le estensioni per Python e Jupyter
- ancora in VSCode:
- creare una cartella di lavoro e renderla la cartella corrente
- copiare nella cartella questo notebook e il file
linearutils.py
e aprire il notebook- creare un ambiente virtuale locale Python (Select Kernel | Python Environments | Create Python Environment | Venv, e scegliere un interprete Python):
- installare i moduli Python richiesti, eseguendo dal terminale:
pip install torch matplotlib ipympl
Supponiamo di disporre di un insieme di dati $(x_i, y_i)$, e da questo di voler costruire una funzione $f$ in grado di prevedere un $y$ per ogni $x$ dato, dunque tale che $y = f(x)$. Per ogni $x_i$, indicheremo dunque con $y_i$ il dato nel training set e con $\hat y_i$ il valore previsto dal modello.
Ipotizziamo che il modello $f$ da costruire sia una funzione lineare, dunque geometricamente una retta $y = k_0 + k_1 x$, così che il modello è una funzione $y = f_{k_0, k_1}(x)$ a due parametri, $k_0$ e $k_1$, di cui occorre stimare i valori: si tratta di un problema noto in statistica come regressione lineare.
Il problema è così semplice che, come è noto, potrebbe essere risolto per via analitica, con il metodo dei minimi quadrati, dunque assumendo che la retta migliore sia quella che minimizza la sua distanza euclidea con i dati $(x_i, y_i)$. Noi risolveremo però il problema per via numerica, costruendo una semplice rete neurale, da addestrare mediante i dati $(x_i, y_i)$, considerati come il training set.
Importiamo prima di tutto i moduli Python necessari.
%matplotlib widget
from linearutils import get_params, plot, polyfit, train
import torch
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
Costruiamo sinteticamente il training set, che supponiamo generato dalla sovrapposizione di valori disposti su una retta e di valori casuali estratti da una distribuzione normale standard, dunque a media 0 e deviazione standard 1, e quindi rappresentiamolo graficamente.
z0, z1 = 10, 1.5 # i "valori veri", e supposti ignoti e da stimare, dei parametri della retta e dunque del modello
x = torch.arange(-5, 5, 0.1).view(-1, 1)
y = z0 + z1 * x + 2 * torch.randn(x.size())
plot(x, y, 'Dati da cui costruire un modello di previsione')
k1 = 1 # i parametri stimati "a occhio" della retta e dunque del modello
k0 = 5
plot(x, y, f'Dati e modello con parametri stimati "a occhio" (k1={k1:.2f}, k0={k0:.2f})', k0=k0, k1=k1)
Ovviamente, potremmo provare a stimare "a occhio" i valori dei parametri...
Come ricordato, sotto l'ipotesi dei minimi quadrati il problema ha una soluzione algebrica, per cui possiamo trovare immediatamente i valori dei due parametri $k_0$ e $k_1$ della retta, e quindi disegnare la retta.
k = polyfit(x, y)
plot(x, y, f'Dati e modello analitico (k1={k[0]:.2f}, k0={k[1]:.2f})', k0=k[1], k1=k[0])
Costruiamo allora una rete neurale con un neurone di input e uno di output, per $x$ e $y$ rispettivamente, che avrà dunque due parametri, $k_0$ e $k_1$, corrispondenti all'intercetta della retta (il bias) e alla pendenza della retta (il "peso" della connessione).
Questa è dunque la struttura della rete:
I valori iniziali di $k_0$ e $k_1$ sono casuali, e quindi la retta di regressione non ha alcuna affidabilità.
model = torch.nn.Linear(1, 1) # definizione della rete visualizzata nella figura sopra, con un neurone di input e uno di output
k = get_params(model, 2)
plot(x, y, f'Dati e modello con parametri casuali (k1={k[1]:.2f}, k0={k[0]:.2f})', model)
Scegliamo come addestrare la rete specificando:
- l'iperparametro tasso di apprendimento (learning rate);
- l'iperparametro numero di epoche;
- il criterio per valutare la qualità delle previsioni, in questo caso l'errore quadratico medio (mean squared error, MSE);
- il metodo di ottimizzazione, in questo caso la back propagation mediante discesa stocastica lungo il gradiente (stochastic gradient descent, SGD).
learning_rate = 0.1
num_epochs = 20
criterion = torch.nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr = learning_rate)
Addestriamo la rete, visualizzando per ogni epoca le stime dei parametri e il valore della funzione di errore, e alla fine il grafico con la retta di regressione risultante.
Per ora operiamo "a scatola chiusa", dunque senza interessarci di come l'addestramento si realizza.
train(model, criterion, optimizer, x, y, num_epochs)
plot(x, y, "Dati e retta di regressione dopo l'addestramento", model)
epoca parametri MSE 1 [1.58, 1.94] 116.43 2 [3.25, 1.23] 75.13 3 [4.59, 1.72] 49.07 4 [5.66, 1.41] 32.55 5 [6.52, 1.62] 22.05 6 [7.21, 1.49] 15.37 7 [7.75, 1.59] 11.1 8 [8.19, 1.53] 8.38 9 [8.54, 1.57] 6.64 10 [8.83, 1.54] 5.52 11 [9.05, 1.56] 4.81 12 [9.23, 1.55] 4.36 13 [9.37, 1.56] 4.07 14 [9.49, 1.56] 3.88 15 [9.58, 1.56] 3.76 16 [9.65, 1.56] 3.68 17 [9.71, 1.56] 3.64 18 [9.76, 1.56] 3.6 19 [9.8, 1.56] 3.58 20 [9.83, 1.56] 3.57
Possiamo ora "aprire la scatola" per vedere come l'addestramento si realizza, e visualizzando la retta di regressione alla fine di ogni epoca, per mostrare la convergenza progressiva verso la retta che meglio stima quella da cui i dati nel training set sono stati generati.
Il codice qui sotto, che tra l'altro contiene vari dettagli concettualmente non così rilevanti, è basato su PyTorch
.
model = torch.nn.Linear(1, 1)
learning_rate = 0.05
num_epochs = 50
criterion = torch.nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr = learning_rate)
xx, yy = x.numpy(), y.numpy()
plt.close('all')
fig, ax = plt.subplots()
ax.scatter(xx, yy, marker='o')
[k0, k1] = get_params(model)
line, = ax.plot(xx, k0 + k1 * xx, color='blue')
def update(epoch):
y1 = model(x) # inferenza: calcola i valori y previsti dal modello
loss = criterion(y1, y) # calcola il valore della funzione di errore
optimizer.zero_grad()
loss.backward() # backpropagation: calcola il gradiente della funzione di errore
optimizer.step() # aggiorna i parametri del modello
mse = torch.mean((y1 - y)**2) # calcola il nuovo valore della funzione di errore
line.set_ydata(y1.detach().numpy())
[k0, k1] = get_params(model)
plt.title(f'epoca={epoch+1}, k0={k0:.2f}, k1={k1:.2f}, MSE={mse:.2f}')
return line,
_ = FuncAnimation(fig, update, frames=num_epochs, interval=100, repeat=False)
plt.show()
Come ultimo passo, possiamo "aprire la scatola" ancora di più, evitare di usare PyTorch
e scrivendo tutto "a mano".
Ricordiamo che la funzione di errore che abbiamo scelto, lo scarto quadratico medio (mean square error, MSE) tra training set e retta / modello, è definita così:
$\sum\limits_{i=1}^n (\hat y_i - y_i)^2 / n = \sum\limits_{i=1}^n (k_0 + k_1 x_i - y_i)^2 / n$
in cui $\hat y_i$ è dunque il valore previsto dal modello per $x_i$.
Per risolvere il problema, consideriamo l'MSE come una funzione dei parametri $k_0$ e $k_1$:
$g(k_0, k_1) = \sum\limits_{i=1}^n (k_0 + k_1 x_i - y_i)^2 / n$.
Trovando i valori di $k_0$ e $k_1$ che rendono minimo il valore della funzione di errore $g$, avremo risolto il problema.
Per via analitica, questo corrisponderebbe a calcolare il gradiente $\nabla g(k_0, k_1)$ di $g(k_0, k_1)$ e quindi trovare i valori dei suoi argomenti per cui $\nabla g(k_0, k_1) = 0$ (per come è definita, è chiaro che la funzione $g$ non ha un massimo, e dunque un suo zero sarebbe un minimo: l'errore minimo, appunto).
A questa idea è anche ispirato l'algoritmo della "discesa lungo il gradiente" della funzione $g$.
Dato il fatto che le componenti di $\nabla g$, cioè le derivate parziali della funzione $g(k_0, k_1)$, sono
$\partial g(k_0, k_1) / \partial k_0 = 2 \sum_{i=1}^n (k_0 + k_1 x_i - y_i) / n$
$\partial g(k_0, k_1) / \partial k_1 = 2 \sum_{i=1}^n (x_i (k_0 + k_1 x_i - y_i) / n)$
e dopo aver assegnato un valore iniziale a $k_0$ e $k_1$ (passo 0):
- per ogni $x_i$ nel training set, si calcola $\hat y_i = k_0 + k_1 x_i$ (passo 1);
- si calcolano le derivate parziali $\partial g / \partial k_0$ e $\partial g / \partial k_1$ (passo 2);
- si modificano i valori di $k_0$ e $k_1$ "verso lo zero" del gradiente (passo 3), dunque
$k_0 = k_0 - z \partial g / \partial k_0$
$k_1 = k_1 - z \partial g / \partial k_1$
dove $z$ determina la rapidità dell'apprendimento (learning rate),
ripetendo il processo per il numero di volte specificato.
L'animazione mostra come, ripetendo il passo dell'addestramento, il valore dell'MSE si riduce, e i valori dei parametri si avvicinano progressivamente ai valori dei parametri che erano stati usati per generare il training set.
from random import random
import numpy as np
xx, yy = x.numpy(), y.numpy()
plt.close('all')
fig2, ax2 = plt.subplots()
ax2.scatter(xx, yy, marker='o')
k0, k1 = random(), random() # passo 0: inizializzazione casuale dei parametri del modello
line2, = ax2.plot(xx, k0 + k1 * xx, color='blue')
lr = 0.05
num_epochs = 50
def update2(epoch):
global k0, k1
y1 = k0 + k1 * xx # passo 1: inferenza: calcola i valori y previsti dal modello
d0 = np.mean(2 * (y1 - yy)) # passo 2: backpropagation: calcola il gradiente della funzione di errore
d1 = np.mean(2 * xx * (y1 - yy))
k0 = k0 - lr * d0 # passo 3: aggiorna i parametri del modello
k1 = k1 - lr * d1
mse = np.mean((y1 - yy)**2) # calcola il nuovo valore della funzione di errore
line2.set_ydata(y1)
plt.title(f'epoca={epoch+1}, k0={k0:.2f}, k1={k1:.2f}, MSE={mse:.2f}')
return line2,
_ = FuncAnimation(fig2, update2, frames=num_epochs, interval=100, repeat=False)
plt.show()