Regressione lineare multipla: forma matriciale ed equazioni normali

Nell'articolo sulla regressione lineare abbiamo stimato il prezzo di un appartamento a partire dalla sola superficie. La retta funzionava bene, ma sbagliava comunque di qualche migliaio di euro: due appartamenti di metratura simile possono costare diversamente per tanti motivi.

Uno di questi è l'età dell'immobile. In questo articolo la aggiungiamo al modello e vediamo come cambiano i conti quando le variabili sono più di una: si passa dalle formule con le sommatorie alla forma matriciale, che vale per qualunque numero di variabili.

I dati

Agli stessi dieci appartamenti aggiungiamo gli anni dalla costruzione:

Guardando i residui della retta del primo articolo, l'appartamento da 70 m² costava 8,5 k€ più del previsto e quello da 65 m² 6,2 k€ meno: il primo ha 10 anni, il secondo 40. L'età sembra spiegare una parte di quegli errori.

La notazione

Con più variabili, $m$ e $q$ non bastano più. Si usa una notazione generale:

• $i$ indica l'osservazione: l'i-esimo appartamento, da 1 a $N$;
• $\varphi_1, \varphi_2, \dots$ sono i regressori (o feature): le variabili di input, qui superficie ed età;
• $\theta_0, \theta_1, \dots$ sono i parametri da stimare. $\theta_0$ è l'intercetta, quella che prima chiamavamo $q$;
• $\varepsilon(i)$ è l'errore: la parte di $y(i)$ che il modello non riesce a spiegare.

Il caso del primo articolo è quello con un solo regressore: $\theta_0 = q$, $\theta_1 = m$.

C'è un trucco che semplifica tutto: si aggiunge un regressore finto $\varphi_0 = 1$, sempre uguale a uno. Così anche l'intercetta diventa un parametro "moltiplicato per qualcosa", e l'intero modello si scrive come un prodotto scalare tra due vettori:

$d$ è il numero di parametri, intercetta compresa. Nel nostro esempio $d = 3$: intercetta, superficie, età.

Con un regressore il modello è una retta, con due è un piano nello spazio tridimensionale (superficie, età, prezzo). Oltre è un iperpiano, che non si può disegnare ma si calcola esattamente allo stesso modo.

Tutto in una matrice

Scrivere un'equazione per ogni appartamento diventa scomodo. Si impilano allora le $N$ osservazioni una sotto l'altra:

Ogni riga di $X$ è un appartamento, ogni colonna un regressore. La prima colonna è tutta di uni: è il $\varphi_0 = 1$ di prima. Le dieci equazioni diventano una sola:

dove $E$ è il vettore dei dieci errori.

La funzione di costo

È lo stesso errore quadratico medio del primo articolo, riscritto con i nuovi simboli. La previsione per l'appartamento $i$ prima era $m\,x_i + q$, ora è $\varphi(i)^\top\theta$; l'MSE prende il nome di funzione di costo e si indica con $J(\theta)$:

La stessa formula, con i vettori

La sommatoria si può scrivere in modo più compatto usando le matrici della sezione precedente:

• $X\theta$ è il vettore delle previsioni: moltiplicare la matrice per $\theta$ calcola $\varphi(i)^\top\theta$ per ogni riga, cioè per ogni appartamento;
• $Y - X\theta$ è quindi il vettore dei residui, uno per appartamento;
• $\lVert v \rVert^2$, la norma al quadrato di un vettore, è la somma dei quadrati dei suoi elementi: $v_1^2 + v_2^2 + \dots + v_N^2$.

Mettendo insieme le tre cose, la somma dei residui al quadrato diventa $\lVert Y - X\theta \rVert^2$:

Il θ migliore

La stima $\hat\theta$ è il vettore che rende $J$ più piccolo possibile, e si scrive così:

Attenzione alla differenza tra le due parole:

• $\min_\theta J(\theta)$ è il valore più basso che l'errore può raggiungere. Nel nostro esempio, 4,28 k€²;
• $\arg\min_\theta J(\theta)$ è il $\theta$ che lo raggiunge, cioè i parametri da usare: $\hat\theta = (32{,}59;\ 2{,}069;\ -0{,}392)$, come vedremo tra poco. "arg" sta per argomento, il valore da mettere dentro la funzione.

A noi interessa il secondo: non quanto vale l'errore minimo, ma quale modello lo produce.

Perché il minimo è unico

Nel primo articolo, tenendo fermo $q$ e muovendo solo $m$, l'MSE disegnava una curva a U. Qui succede lo stesso con qualunque parametro: ogni residuo $y(i) - \varphi(i)^\top\theta$ è una funzione lineare di $\theta$, il suo quadrato è una parabola rivolta verso l'alto, e una somma di parabole rivolte verso l'alto è ancora una parabola rivolta verso l'alto.

Se $J$ è una U in ogni direzione, nel suo insieme è una scodella. Una scodella ha un solo punto più basso: non ci sono "buche" secondarie in cui un algoritmo potrebbe fermarsi credendo di aver trovato il minimo. In matematica questa proprietà si chiama convessità.

C'è un'eccezione. Se due colonne di $X$ contengono la stessa informazione, la scodella ha il fondo piatto, come una grondaia, e i minimi sono infiniti. È il caso dei regressori ridondanti, che vedremo più avanti.

Le equazioni normali

Nel punto più basso della scodella, la funzione di costo non diminuisce più in nessuna direzione: il suo gradiente è nullo. Possiamo sfruttare questa proprietà per trovare direttamente i parametri che minimizzano l'errore, senza procedere per tentativi come nella discesa del gradiente.

Partiamo dalla funzione di costo della regressione lineare multipla:

Sviluppando il quadrato, otteniamo:

Deriviamo la funzione di costo

Abbiamo sviluppato la funzione di costo in modo da separare i termini costanti da quelli che dipendono dai parametri $\theta$:

Ora vogliamo trovare i valori di $\theta$ che rendono $J(\theta)$ minima. Per farlo, calcoliamo il gradiente, cioè il vettore delle derivate parziali rispetto a tutti i parametri del modello.

1. Due regole di derivazione matriciale

Ci servono due regole, che generalizzano quelle delle funzioni a una variabile:

Nella prima formula, $b$ è un vettore di costanti: la derivata di una funzione lineare restituisce i suoi coefficienti. È l'equivalente di:

La seconda formula riguarda invece una funzione quadratica. È l'equivalente della regola:

La differenza è che, nel caso matriciale, la derivata contiene $A + A^\top$, perché il vettore $x$ compare sia a sinistra sia a destra della matrice.

Se $A$ è simmetrica, cioè $A^\top = A$, la formula si semplifica:

2. Applichiamo le regole alla funzione di costo

Riprendiamo i tre termini di $J(\theta)$ e deriviamoli uno alla volta rispetto a $\theta$.

Primo termine: $Y^\top Y$

Questo termine non dipende da $\theta$: contiene soltanto i dati osservati $Y$, che consideriamo costanti durante l'addestramento. La sua derivata è quindi zero.

Secondo termine: $-2\,\theta^\top X^\top Y$

Applichiamo la prima regola, identificando $x = \theta$ e $b = X^\top Y$:

Il fattore $-2$ rimane davanti alla derivata:

Terzo termine: $\theta^\top X^\top X \theta$

Applichiamo la seconda regola, identificando $x = \theta$ e $A = X^\top X$:

La matrice $X^\top X$ è simmetrica, perché:

I due termini nella parentesi sono quindi uguali. Possiamo sommarli:

3. Mettiamo insieme le derivate

La derivata di una somma è la somma delle derivate. Inoltre, il fattore costante $1/N$ rimane davanti:

Abbiamo così ottenuto il gradiente della funzione di costo.

4. Poniamo il gradiente uguale a zero

Nel punto di minimo, il gradiente è nullo. Poniamo quindi:

Moltiplichiamo entrambi i membri per $N/2$, eliminando i fattori costanti:

Portiamo il primo termine a destra:

Queste sono le equazioni normali: abbiamo trasformato il problema di minimizzazione della funzione di costo in un sistema lineare di equazioni.

La soluzione in forma chiusa

Se $X^\top X$ è invertibile, possiamo moltiplicare entrambi i membri delle equazioni normali per $(X^\top X)^{-1}$:

Poiché $(X^\top X)^{-1} X^\top X = I$, dove $I$ è la matrice identità, rimane:

Questa è la soluzione dei minimi quadrati in forma chiusa: una formula che permette di calcolare direttamente i parametri ottimali, senza aggiornamenti iterativi.

Se $X^\top X$ non è invertibile, la formula con l'inversa non è applicabile direttamente. Il problema può comunque avere soluzioni, eventualmente non uniche.

Con un solo regressore e un'intercetta, questa formula restituisce gli stessi parametri $m$ e $q$ della regressione lineare semplice. Cambia la notazione, ma il principio è identico.

import numpy as np

superficie = np.array([40, 50, 55, 65, 70, 80, 90, 100, 110, 120])
eta = np.array([20, 45, 25, 40, 10, 30, 15, 35, 20, 40])
prezzo = np.array([105, 120, 140, 150, 175, 185, 210, 225, 255, 265])

# la colonna di uni è il regressore finto per l'intercetta
X = np.column_stack([np.ones(len(prezzo)), superficie, eta])

theta = np.linalg.solve(X.T @ X, X.T @ prezzo)
print(theta.round(3))                           # [32.595  2.069 -0.392]
print(round(np.array([1, 85, 20]) @ theta, 1))  # 200.6

Quando XᵀX non è invertibile

La formula funziona solo se $X^\top X$ si può invertire. Ci sono due casi tipici in cui non succede.

Regressori ridondanti

Supponiamo di aggiungere al dataset anche la superficie in piedi quadrati: è la stessa informazione della superficie in metri quadrati. La quarta colonna di $X$ è un multiplo esatto della seconda, quindi le colonne sono linearmente dipendenti. Le soluzioni sono infinite: un m² in più può essere "pagato" dal parametro dei metri, da quello dei piedi o da un mix dei due, e la previsione non cambia.

$X^\top X$ è singolare, o meglio, a causa degli arrotondamenti, quasi singolare. Ed è qui che sta il pericolo: np.linalg.inv spesso non dà errore, ma restituisce numeri insensati. Sul mio computer:

Con l'inversa "forzata" la funzione di costo esplode da 4,3 a oltre 11.000 k€². I valori esatti dipendono dalla libreria e dal processore, il che è un altro sintomo del problema.

Esiste anche una versione più subdola del problema: regressori non identici ma molto correlati, come la superficie e il numero di locali. $X^\top X$ è invertibile, ma i parametri diventano instabili e difficili da interpretare: il modello non sa a chi dei due attribuire l'effetto.

Troppi pochi dati

Se le osservazioni sono meno dei parametri, o tante quante ($N \le d$), i dati non bastano a determinare il modello. Con soli 3 appartamenti e 3 parametri il piano passa esattamente per tutti e tre i punti, con errore zero:

Errore zero non vuol dire modello perfetto: vuol dire che il modello ha memorizzato i dati invece di imparare la relazione. Con 3 appartamenti il valore di un m² sale a 2.460 €, un numero che non regge su nessun altro immobile. È il primo assaggio dell'overfitting.

Il limite delle equazioni normali

$X^\top X$ è una matrice $d \times d$, e risolvere il sistema costa all'incirca $d^3$ operazioni. Con 3 parametri è istantaneo, con qualche migliaio ancora fattibile. Con milioni di feature, come nei modelli per testo o immagini, diventa impraticabile, anche solo per la memoria necessaria a contenere la matrice.

In quei casi si rinuncia alla soluzione esatta e si usa un metodo iterativo: la discesa del gradiente, già vista nel primo articolo. Con più variabili però nasconde una trappola, e il prossimo articolo è dedicato a quella: discesa del gradiente e normalizzazione dei dati.

Riepilogo

• Con più regressori il modello si scrive $y(i) = \varphi(i)^\top \theta + \varepsilon(i)$, con un regressore finto $\varphi_0 = 1$ per l'intercetta.
• Impilando le osservazioni si ottiene $Y = X\theta + E$: ogni riga di $X$ è un'osservazione, ogni colonna un regressore.
• La funzione di costo $J(\theta) = \frac{1}{N} \lVert Y - X\theta \rVert^2$ è convessa: ha un unico minimo.
• Annullando il gradiente si ottengono le equazioni normali $X^\top X \theta = X^\top Y$, da cui $\hat\theta = (X^\top X)^{-1} X^\top Y$.
• $X^\top X$ non è invertibile con regressori ridondanti (si toglie la colonna, o si usa la pseudo-inversa) o con troppi pochi dati ($N \le d$): il modello memorizza invece di imparare.
• Con moltissimi parametri le equazioni normali diventano troppo costose e si passa alla discesa del gradiente.