02 - Regressione Lineare from scratch in Python

Hai una tabella di esempi, ciascuno con alcune misure e un numero da prevedere, e vuoi una formula che dalle misure restituisca quel numero. La regressione lineare è la formula più semplice che si può scrivere: una somma pesata delle misure, più un termine costante. Scegliere i pesi è un problema di minimi quadrati e ha una soluzione esatta; il lavoro vero sta nel calcolarla in modo affidabile. Questo articolo la ricava, la calcola in tre modi che danno gli stessi coefficienti, mostra con una misura quando uno dei tre fallisce, e la ottiene anche per tentativi con la discesa del gradiente, controllando il gradiente contro un calcolo numerico. Serve il capitolo 01 per le parole «train», «test» e «baseline».

Tutto gira su numpy 2.4.6 e scikit-learn 1.9.0, con il dataset Diabetes incluso nella libreria: 442 pazienti, dieci misure di base (età, sesso, indice di massa corporea, pressione media e sei esami del sangue) e come valore da prevedere un indice quantitativo della progressione della malattia a un anno dalla misura (descrizione della documentazione di scikit-learn). Nessun download; ogni uscita riportata è identica in due esecuzioni.

Il modello e la perdita

Indichiamo con X la matrice dei dati, una riga per esempio e una colonna per misura, con n righe e d colonne, e con y il vettore dei valori da prevedere. Il modello prevede ŷ = Xβ, dove β contiene un peso per colonna. Per avere anche un termine costante (l'intercetta: la previsione quando tutte le misure valgono zero) si aggiunge a X una colonna di uni, il cui peso diventa appunto l'intercetta. Un esempio con numeri: per prevedere il prezzo di un appartamento da metri quadrati e piano, β avrà tre pesi, uno per i metri quadrati, uno per il piano e l'intercetta.

Si sceglie β minimizzando l'errore quadratico medio (MSE, mean squared error), la media dei quadrati dei residui, cioè delle differenze fra valore vero e previsto: MSE(β) = (1/n)·||y − Xβ||². Elevare al quadrato rende ogni errore positivo e punisce molto di più gli errori grandi, il che ha un lato utile (la funzione è liscia e ha un solo minimo) e uno scomodo (un singolo valore anomalo tira la retta verso di sé). Per giudicare il modello si usa anche R², uguale a 1 meno il rapporto fra la somma dei quadrati dei residui del modello e quella di un modello che prevede sempre la media di y: vale 1 per un adattamento perfetto, 0 se non si fa meglio della media, e può diventare negativo sul test.

Le equazioni normali e la geometria dei residui

Il gradiente di ||y − Xβ||² rispetto a β è −2Xᵀ(y − Xβ). Porlo uguale a zero dà le equazioni normali: XᵀX β = Xᵀy. Se le colonne di X sono linearmente indipendenti, XᵀX è invertibile e la soluzione è unica; la perdita è una parabola in più dimensioni e il punto in cui il gradiente si annulla è il suo minimo globale. Il manuale di James, Witten, Hastie e Tibshirani (An Introduction to Statistical Learning) ricava la stessa soluzione nel paragrafo 3.1.1, «Estimating the Coefficients».

Il significato geometrico aiuta a ricordare la formula: Xβ̂ è la proiezione ortogonale di y sullo spazio generato dalle colonne di X, quindi i residui sono perpendicolari a ogni colonna, cioè Xᵀ(y − Xβ̂) = 0. Con la colonna di uni, in particolare, i residui sommano a zero. Lo verifichiamo su un esercizio da fare a mano.

Esercizio svolto: tre punti

Dati i punti (1, 1), (2, 2), (3, 2), trova la retta y = b₀ + b₁x che minimizza il MSE. A mano: XᵀX vale la matrice con righe (3, 6) e (6, 14), perché ci sono tre punti, la somma delle x è 6 e la somma dei quadrati è 14; Xᵀy vale (5, 11). Il determinante è 3·14 − 6·6 = 6, quindi β̂ = (1/6)·(14·5 − 6·11, −6·5 + 3·11) = (4/6, 3/6), cioè intercetta 0,667 e pendenza 0,5. Il programma conferma il conto e la perpendicolarità dei residui.

import numpy as np

x = np.array([1.0, 2.0, 3.0])
y = np.array([1.0, 2.0, 2.0])
X = np.c_[np.ones(3), x]
print("X'X =", (X.T @ X).tolist())
print("X'y =", (X.T @ y).tolist())
print("determinante:", round(float(np.linalg.det(X.T @ X)), 6))
b = np.linalg.solve(X.T @ X, X.T @ y)
print("intercetta, pendenza:", b.round(4).tolist())
print("residui:", (y - X @ b).round(4).tolist(), " somma:", round(float((y - X @ b).sum()), 12))
print("X' residui (ortogonalita'):", (X.T @ (y - X @ b)).round(12).tolist())

# Output:
# X'X = [[3.0, 6.0], [6.0, 14.0]]
# X'y = [5.0, 11.0]
# determinante: 6.0
# intercetta, pendenza: [0.6667, 0.5]
# residui: [-0.1667, 0.3333, -0.1667]  somma: 0.0
# X' residui (ortogonalita'): [0.0, 0.0]

I residui sono −0,1667, 0,3333 e −0,1667: sommano a zero e, pesati con x, a zero anch'essi, che sono le due condizioni Xᵀr = 0 scritte per esteso. Se cambi un punto, per esempio il terzo in (3, 5), la pendenza passa da 0,5 a 2 e l'intercetta da 0,667 a −1,333 (stesso calcolo, con Xᵀy uguale a (8, 20)): è la sensibilità ai valori anomali di cui sopra.

Tre modi di calcolare la stessa soluzione

Il metodo di scuola risolve le equazioni normali. Un metodo numericamente più sicuro, numpy.linalg.lstsq, risolve il problema dei minimi quadrati senza formare XᵀX, usando una decomposizione ai valori singolari (SVD); scikit-learn, dice la sua documentazione nel paragrafo sulla complessità degli OLS, calcola la soluzione con la SVD di X, con un costo proporzionale a n·d² se n ≥ d. Il listato li confronta sullo stesso dataset, dopo aver messo da parte il 25% degli esempi.

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X, y = load_diabetes(return_X_y=True, scaled=False)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.25, random_state=0)
sc = StandardScaler().fit(Xtr)                       # media e scarto dal solo training
A = np.c_[np.ones(len(Xtr)), sc.transform(Xtr)]      # colonna di 1 = intercetta
A_te = np.c_[np.ones(len(Xte)), sc.transform(Xte)]

b_eq = np.linalg.solve(A.T @ A, A.T @ ytr)           # equazioni normali
b_ls = np.linalg.lstsq(A, ytr, rcond=None)[0]        # minimi quadrati via SVD
sk = LinearRegression().fit(sc.transform(Xtr), ytr)
b_sk = np.r_[sk.intercept_, sk.coef_]

print("righe train/test, feature:", len(Xtr), len(Xte), X.shape[1])
print("max |equazioni normali - lstsq|:", np.abs(b_eq - b_ls).max() < 1e-9)
print("max |lstsq - sklearn|:", np.abs(b_ls - b_sk).max() < 1e-9)
print("intercetta (= media di y nel train):", round(b_ls[0], 3), round(ytr.mean(), 3))
print("R2 train:", round(r2_score(ytr, A @ b_ls), 3), " R2 test:", round(r2_score(yte, A_te @ b_ls), 3))

# Output:
# righe train/test, feature: 331 111 10
# max |equazioni normali - lstsq|: True
# max |lstsq - sklearn|: True
# intercetta (= media di y nel train): 151.921 151.921
# R2 train: 0.555  R2 test: 0.359

Le tre soluzioni coincidono a meno di 1e-9: sono matematicamente lo stesso numero, calcolato con tre algoritmi. In np.linalg.solve(A.T @ A, A.T @ ytr) uso solve e non l'inversa esplicita, perché risolvere il sistema è più preciso e più economico che invertire e moltiplicare. Lo StandardScaler è calcolato sul solo training, e le stesse media e scarto sono applicati al test; con colonne centrate l'intercetta coincide con la media di y sul training (151,921 per entrambi i valori stampati). Il valore da leggere con più cura è l'R²: 0,555 sul training e 0,359 sul test. Il modello spiega meno sul test, ed è normale; con soli 111 esempi di test la stima è inoltre rumorosa, quindi non confrontare queste due cifre con quelle di un altro seme. Il capitolo 05 mostra come stimare la variabilità.

Quando le equazioni normali tradiscono

Il numero di condizionamento κ(X), il rapporto fra il più grande e il più piccolo valore singolare, dice quanto un piccolo errore sui dati si amplifica nella soluzione. Formare XᵀX eleva al quadrato quel numero: κ(XᵀX) = κ(X)². Con numeri in doppia precisione (circa 16 cifre decimali) un condizionamento di 1e18 lascia zero cifre sicure. Per vederlo costruiamo una situazione in cui la risposta giusta è nota: i dati sono un polinomio con tutti i coefficienti uguali a 1, calcolato senza rumore su 100 punti in [0, 1], e confrontiamo l'errore massimo sui coefficienti dei due metodi.

import numpy as np

x = np.linspace(0, 1, 100)
for grado in (5, 10, 14):
    X = np.vander(x, grado + 1, increasing=True)      # 1, x, x^2, ..., x^grado
    beta = np.ones(grado + 1)
    y = X @ beta                                      # dati senza rumore: la soluzione vera e' nota
    b_eq = np.linalg.solve(X.T @ X, X.T @ y)
    b_ls = np.linalg.lstsq(X, y, rcond=None)[0]
    print(f"grado {grado:2d}  cond(X)={np.linalg.cond(X):.1e}  cond(X'X)={np.linalg.cond(X.T @ X):.1e}"
          f"  errore eq.normali={np.abs(b_eq - beta).max():.1e}  errore lstsq={np.abs(b_ls - beta).max():.1e}")

# Output:
# grado  5  cond(X)=3.7e+03  cond(X'X)=1.4e+07  errore eq.normali=1.8e-10  errore lstsq=1.2e-13
# grado 10  cond(X)=2.1e+07  cond(X'X)=4.5e+14  errore eq.normali=2.4e-03  errore lstsq=2.8e-10
# grado 14  cond(X)=2.3e+10  cond(X'X)=2.0e+18  errore eq.normali=4.8e+01  errore lstsq=2.1e-06

Fino al grado 5 le equazioni normali funzionano (errore 1,8e-10). Al grado 10 l'errore sale a 2,4e-3 mentre lstsq resta a 2,8e-10; al grado 14 le equazioni normali sbagliano di 48 su coefficienti che valgono 1, cioè la risposta è insensata, e lstsq è ancora a 2e-6. Il condizionamento di X è 2,3e10 e quello di XᵀX è 2,0e18, cioè oltre la precisione della macchina. Le ultime cifre di queste uscite dipendono dalla libreria di algebra lineare con cui è compilato NumPy, ma gli ordini di grandezza e il punto in cui il metodo ingenuo cede sono robusti. Questa è la ragione per cui scikit-learn non risolve mai le equazioni normali, e la ragione per cui conviene standardizzare le misure: colonne di scale diverse peggiorano il condizionamento.

Discesa del gradiente

Se n è enorme, o i dati arrivano a flusso, non si forma XᵀX e si cerca il minimo per tentativi: a ogni passo si sposta β di un tratto η (il learning rate) nella direzione opposta al gradiente, β ← β − η·∇MSE(β), con ∇MSE = (2/n)·Xᵀ(Xβ − y). Prima di fidarsi di un gradiente scritto a mano si controlla contro le differenze finite centrate, (f(β+h·eⱼ) − f(β−h·eⱼ)) / 2h, che richiedono solo di calcolare la perdita. Il test cattura quasi tutti gli errori di segno e di fattore 2 o n.

Il passo massimo è fissato dalla curvatura: l'Hessiana del MSE è (2/n)·XᵀX, e il metodo converge solo se η < 2/λ_max, dove λ_max è il suo autovalore maggiore. Ogni direzione propria dell'errore viene moltiplicata ad ogni passo per (1 − η·λᵢ): se il modulo supera 1 l'errore cresce, e la direzione più ripida è la prima a uscire dal limite.

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.preprocessing import StandardScaler

X, y = load_diabetes(return_X_y=True, scaled=False)
A = np.c_[np.ones(len(X)), StandardScaler().fit_transform(X)]
n, d = A.shape

def mse(b):  return np.mean((A @ b - y) ** 2)
def grad(b): return 2 / n * A.T @ (A @ b - y)

# 1) il gradiente analitico contro le differenze finite centrate
rng = np.random.default_rng(0)
b = rng.standard_normal(d)
h, num = 1e-5, np.zeros(d)
for j in range(d):
    e = np.zeros(d); e[j] = h
    num[j] = (mse(b + e) - mse(b - e)) / (2 * h)
rel = np.linalg.norm(num - grad(b)) / np.linalg.norm(grad(b))
print("errore relativo gradiente analitico vs differenze finite < 1e-6:", rel < 1e-6)

# 2) soglia del passo: eta < 2 / lambda_max dell'Hessiana (2/n) A^T A
lam = np.linalg.eigvalsh(2 / n * A.T @ A).max()
soglia = 2 / lam
b_chiusa = np.linalg.lstsq(A, y, rcond=None)[0]
print("lambda_max Hessiana:", round(lam, 3), " soglia eta:", round(soglia, 4))
for frac in (0.9, 1.1):
    b = np.zeros(d)
    for _ in range(500):
        b = b - frac * soglia * grad(b)
    print(f"eta = {frac} x soglia, dopo 500 passi: MSE =", f"{mse(b):.4g}")
b = np.zeros(d)
for t in range(5000):
    b = b - 0.9 * soglia * grad(b)
print("eta = 0.9 x soglia, dopo 5000 passi: MSE", round(mse(b), 3), "contro forma chiusa", round(mse(b_chiusa), 3))
print("coefficienti uguali a 1e-6:", np.abs(b - b_chiusa).max() < 1e-6)

# Output:
# errore relativo gradiente analitico vs differenze finite < 1e-6: True
# lambda_max Hessiana: 8.048  soglia eta: 0.2485
# eta = 0.9 x soglia, dopo 500 passi: MSE = 2860
# eta = 1.1 x soglia, dopo 500 passi: MSE = 2.776e+82
# eta = 0.9 x soglia, dopo 5000 passi: MSE 2859.696 contro forma chiusa 2859.696
# coefficienti uguali a 1e-6: True

Il controllo passa: l'errore relativo fra gradiente analitico e numerico è sotto 1e-6. Per le misure standardizzate dell'esempio λ_max vale 8,048 e la soglia 0,2485. Con η pari al 90% della soglia, 500 passi portano il MSE a 2860, e dopo 5000 passi il risultato coincide con la forma chiusa (2859,696) con coefficienti uguali a meno di 1e-6; con il 110% della soglia il MSE dopo 500 passi vale 2,8e82, quindi il metodo diverge. La lentezza della convergenza dipende invece dall'autovalore minore: più sono sbilanciati gli autovalori (cioè più è alto il condizionamento), più passi servono. Per questo la discesa del gradiente su dati non standardizzati può richiedere molte più iterazioni. Le varianti a mini-batch e il momento non sono eseguite qui; il capitolo 10 le usa per le reti neurali.

Regolarizzazione: Ridge e Lasso

Quando le colonne sono quasi dipendenti, XᵀX è quasi singolare e i coefficienti diventano instabili. Si aggiunge allora una penalità sulla grandezza dei pesi. La regressione Ridge minimizza ||y − Xβ||² + α·||β||² (norma L2) e ha ancora una soluzione chiusa: (XᵀX + αI)β = Xᵀy. Per α > 0 la matrice è sempre invertibile. Il Lasso usa la norma L1, la somma dei valori assoluti, che non è derivabile in zero: non esiste una formula chiusa, ma la discesa per coordinate risolve un peso alla volta con la funzione soft-thresholding S(z, t) = segno(z)·max(|z| − t, 0), e porta alcuni pesi esattamente a zero.

Attenzione alla convenzione. Secondo la documentazione di scikit-learn, Ridge minimizza ||Xw − y||² + α||w||², mentre Lasso minimizza (1/(2n))·||Xw − y||² + α||w||₁: lo stesso α non ha lo stesso peso nei due casi, e un Lasso scritto a mano con l'altra convenzione dà risultati diversi da quelli di Lasso(alpha=...). Nel listato il Lasso a coordinate segue la convenzione di scikit-learn, e il confronto lo verifica.

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.linear_model import Lasso, Ridge
from sklearn.preprocessing import StandardScaler

X, y = load_diabetes(return_X_y=True, scaled=False)
X = StandardScaler().fit_transform(X)
yc = y - y.mean()                                    # dati centrati: niente intercetta
n, d = X.shape

# Ridge: min ||y - Xb||^2 + a ||b||^2  ->  (X'X + a I) b = X'y
a = 10.0
b_ridge = np.linalg.solve(X.T @ X + a * np.eye(d), X.T @ yc)
sk_r = Ridge(alpha=a, fit_intercept=False).fit(X, yc)
print("Ridge uguale a sklearn:", np.allclose(b_ridge, sk_r.coef_, atol=1e-8))

# Lasso: min (1/2n)||y - Xb||^2 + a |b|_1, discesa per coordinate con soft-thresholding
def soft(z, t): return np.sign(z) * np.maximum(np.abs(z) - t, 0)
def lasso_cd(X, y, a, iters=2000):
    b = np.zeros(X.shape[1]); n = len(y)
    for _ in range(iters):
        for j in range(X.shape[1]):
            r = y - X @ b + X[:, j] * b[j]            # residuo senza la feature j
            b[j] = soft(X[:, j] @ r / n, a) / (X[:, j] @ X[:, j] / n)
    return b
a_l = 5.0
b_l = lasso_cd(X, yc, a_l)
sk_l = Lasso(alpha=a_l, fit_intercept=False, tol=1e-10, max_iter=100000).fit(X, yc)
print("Lasso uguale a sklearn:", np.allclose(b_l, sk_l.coef_, atol=1e-6))
print("coefficienti a zero (Lasso):", int((b_l == 0).sum()), "su", d)
print("norma L1 OLS vs Ridge vs Lasso:",
      [round(float(np.abs(v).sum()), 1) for v in (np.linalg.lstsq(X, yc, rcond=None)[0], b_ridge, b_l)])

# Output:
# Ridge uguale a sklearn: True
# Lasso uguale a sklearn: True
# coefficienti a zero (Lasso): 5 su 10
# norma L1 OLS vs Ridge vs Lasso: [164.6, 105.0, 65.0]

I due confronti con scikit-learn passano (tolleranze 1e-8 e 1e-6). Il Lasso con alpha=5.0 azzera 5 coefficienti su 10, cioè sceglie da solo quali misure tenere; i pesi di OLS, Ridge e Lasso hanno somma dei valori assoluti rispettivamente 164,6, 105,0 e 65,0: la penalità restringe i pesi, ed è il fenomeno chiamato shrinkage. Si centra y e si standardizza X perché la penalità tratta tutti i pesi allo stesso modo: senza standardizzare, una misura espressa in unità piccole avrebbe un peso grande e verrebbe penalizzata più delle altre. Come scegliere α è il tema del capitolo 04, che lega la penalità al compromesso fra bias e varianza.

Quando non usare la regressione lineare

• La relazione fra misure e valore è chiaramente curva o ha soglie: i residui mostrano allora un andamento sistematico se li metti in grafico contro le previsioni, e servono termini polinomiali oppure alberi (capitolo 06).
• Il valore da prevedere è una categoria: si usa la regressione logistica, perché una retta può restituire «probabilità» minori di 0 o maggiori di 1.
• Le misure sono più numerose degli esempi (d > n): XᵀX non è invertibile, e serve una penalità.

Le tre voci non sono un elenco di difetti ma tre segnali diversi da cercare nei dati; il primo si vede nei residui, il secondo nel tipo di target, il terzo nella forma della tabella. In tutti e tre i casi il rimedio comincia dal confronto con una baseline e da una valutazione su dati non visti, come nel capitolo 01.

Che cosa sai fare adesso

Sai ricavare le equazioni normali e interpretarle come una proiezione, calcolare i coefficienti con tre algoritmi e sapere quando il primo non è affidabile, verificare un gradiente con le differenze finite e fissare il passo della discesa dal massimo autovalore, e usare Ridge e Lasso con la convenzione giusta. Nel capitolo successivo il bersaglio non è più un numero ma una probabilità: 03 - Regressione logistica e classificazione binaria.

Nella stessa serie trovi anche: 12 - Deployment ML in produzione: serving e monitoring.