Metodo dei minimi quadrati

Lezione 34, prima parte (videolezione n. 34). Esercizio svolto: Esercizio 87 · retta dei minimi quadrati per quattro punti. Prerequisiti: Complemento ortogonale e proiezioni ortogonaliL'ortogonale U⊥ di un sottospazio è un sottospazio di dimensione n − dim U, e R^n = U ⊕ U⊥; ogni vettore si scompone in proiezione su U più componente ortogonale; la proiezione è il punto di U più vicino e si calcola con un sistema o con la matrice di proiezione A(AᵀA)⁻¹Aᵀ.Complemento ortogonale e proiezioni ortogonali →, Sistemi lineari e teorema di Rouché-CapelliUn sistema lineare si scrive AX = B; ha soluzioni se e solo se B sta nell'immagine di A, cioè se rango(A) = rango(A|B) (Rouché-Capelli); le soluzioni sono una soluzione particolare più il nucleo e dipendono da n − r parametri.Sistemi lineari e teorema di Rouché-Capelli →. Seguito: Basi ortonormali e Gram-SchmidtUna base ortonormale è fatta di vettori di norma 1 a due a due ortogonali: le coordinate si calcolano con prodotti scalari e le proiezioni con una formula diretta. Il procedimento di Gram-Schmidt trasforma una base qualsiasi in una base ortogonale (e poi ortonormale) dello stesso sottospazio.Basi ortonormali e Gram-Schmidt →.

Il problema

Negli esperimenti si raccolgono molti dati e si cerca un modello con pochi parametri che li descriva. Esempio tipico: si misurano $m$ coppie $(x_i, y_i)$ e si sospetta una legge lineare $y = \beta x + \gamma$. Imponendo che ogni punto stia sulla retta si ottiene un sistema con $m$ equazioni (tante quante le misure) e solo $2$ incognite ($\beta$, $\gamma$):

$$\beta x_i + \gamma = y_i, \qquad i = 1, \dots, m .$$

Per gli errori di misura i punti non sono mai perfettamente allineati, e il sistema non ha soluzioni. Rinunciare non serve: si cerca la retta che "sbaglia di meno".

In generale: dato un sistema $AX = b$ con $A \in M_{m \times n}(\mathbb{R})$ (tipicamente $m > n$, sistema sovradeterminato) e senza soluzioni, si cerca $X$ tale che $AX$ sia il più vicino possibile a $b$.

Definizione. Una soluzione ai minimi quadrati di $AX = b$ è un vettore $\hat X \in \mathbb{R}^n$ che rende minima la norma $$\|A X - b\|^2 = \sum_{i=1}^m \big( (AX)_i - b_i \big)^2 .$$

Esempio. $A = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$, $b = \begin{pmatrix} 1 \\ 3 \end{pmatrix}$, cioè il sistema $x = 1$, $x = 3$ (incompatibile): $\|AX - b\|^2 = (x - 1)^2 + (x - 3)^2$ è minima per $\hat{x} = 2$. Passaggi: sviluppando, $(x - 1)^2 + (x - 3)^2 = 2x^2 - 8x + 10$, una parabola con il vertice in $x = -\frac{-8}{2 \cdot 2} = 2$; oppure, con le derivate, $\frac{d}{dx}(2x^2 - 8x + 10) = 4x - 8 = 0$ per $x = 2$ (teorema di Fermatx0 è punto di minimo (massimo) relativo se f(x0) ≤ f(x) (≥) per gli x del dominio vicini a x0. I candidati sono gli estremi del dominio, i punti dove f non è derivabile e i punti interni con f'(x0) = 0 (punti critici o stazionari). Teorema di Fermat: in un punto interno di minimo o massimo relativo dove f è derivabile, f'(x0) = 0. È solo una condizione necessaria: x³ in 0.Massimi e minimi relativi e teorema di Fermat →). Il valore minimo è $1^2 + 1^2 = 2$: il compromesso sbaglia di $1$ in ciascuna equazione.

Il nome viene dal fatto che si minimizza una somma di quadrati degli scarti.

La soluzione geometrica

Il vettore $AX = x_1 C_1 + \dots + x_n C_n$ (dove $C_j$ sono le colonne di $A$) è una combinazione lineare delle colonne, cioè un vettore di $\operatorname{Im} A$. Al variare di $X$, $AX$ descrive tutto il sottospazio $\operatorname{Im} A \subseteq \mathbb{R}^m$ (l'immagineIl nucleo (vettori mandati in 0) e l'immagine (vettori raggiunti) di una funzione lineare sono sottospazi; f è iniettiva se e solo se Ker f = {0}; dim Ker f + dim Im f = dim V (nullità + rango); l'antimmagine di un vettore è una soluzione particolare più il nucleo.Nucleo e immagine → di $X \mapsto AX$ è generata dalle colonne). Il sistema ha soluzione esattamente quando $b \in \operatorname{Im} A$: è la stessa cosa che dice il teorema di Rouché-CapelliUn sistema lineare si scrive AX = B; ha soluzioni se e solo se B sta nell'immagine di A, cioè se rango(A) = rango(A|B) (Rouché-Capelli); le soluzioni sono una soluzione particolare più il nucleo e dipendono da n − r parametri.Sistemi lineari e teorema di Rouché-Capelli →, perché aggiungere la colonna $b$ non fa salire il rango proprio quando $b$ è combinazione delle colonne.

Quindi il problema diventa: trovare il vettore di $\operatorname{Im} A$ più vicino a $b$. Per il teorema della migliore approssimazioneIl vettore di un sottospazio U più vicino a b è la proiezione ortogonale di b su U.Complemento ortogonale e proiezioni ortogonali → la risposta è la proiezione ortogonale:

$$A \hat X = p_{\operatorname{Im} A}(b) .$$

Le equazioni normali

La proiezione è caratterizzata dal fatto che il residuo $b - A\hat X$ è ortogonale a $\operatorname{Im} A$, cioè a ogni colonna $C_j$ di $A$:

$$C_j \cdot (b - A \hat X) = 0 \quad (j = 1, \dots, n) \qquad \Longleftrightarrow \qquad A^T (b - A \hat X) = 0$$

(la $j$-esima componente di $A^T w$ è proprio $C_j \cdot w$, perché le righe di $A^T$ sono le colonne di $A$).

Teorema. $\hat X$ è una soluzione ai minimi quadrati di $AX = b$ se e solo se risolve le equazioni normali $$A^T A\, \hat X = A^T b .$$ Se le colonne di $A$ sono linearmente indipendenti, $A^T A$ è invertibile e la soluzione è unica: $$\hat X = (A^T A)^{-1} A^T b .$$

Esempio. Per lo stesso sistema $A^T A = (2)$ e $A^T b = (4)$, quindi $\hat{x} = (A^T A)^{-1} A^T b = 2$.

Osservazioni.

• Il sistema normale è quadrato ($n \times n$) e ha sempre soluzione (la proiezione esiste sempre).
• Se $AX = b$ ha soluzioni vere, le soluzioni ai minimi quadrati sono esattamente quelle (residuo nullo).
• $A^T A$ è simmetrica; è invertibile quando le colonne di $A$ sono indipendenti: se $A^T A x = 0$ allora $\|Ax\|^2 = x^T A^T A x = 0$, quindi $Ax = 0$ e $x = 0$.
• Il vettore $A\hat X = A (A^TA)^{-1} A^T b$ contiene la matrice di proiezioneL'ortogonale U⊥ di un sottospazio è un sottospazio di dimensione n − dim U, e R^n = U ⊕ U⊥; ogni vettore si scompone in proiezione su U più componente ortogonale; la proiezione è il punto di U più vicino e si calcola con un sistema o con la matrice di proiezione A(AᵀA)⁻¹Aᵀ.Complemento ortogonale e proiezioni ortogonali → $P = A (A^TA)^{-1} A^T$ su $\operatorname{Im} A$.

Esempio: un sistema impossibile

$$\begin{cases} x = 1 \\ y = 1 \\ x + y = 0 \end{cases} \qquad A = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}, \quad b = \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}.$$

È impossibile: le prime due danno $x + y = 2 \ne 0$. Equazioni normali:

$$A^T A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}, \qquad A^T b = \begin{pmatrix} 1 \\ 1 \end{pmatrix}, \qquad \begin{cases} 2x + y = 1 \\ x + 2y = 1 \end{cases} \Longrightarrow\ x = y = \frac13 .$$

Da dove vengono i numeri: l'elemento $(i, j)$ di $A^T A$ è il prodotto scalare tra la colonna $i$ e la colonna $j$ di $A$. Con $C_1 = (1, 0, 1)$ e $C_2 = (0, 1, 1)$: $C_1 \cdot C_1 = 2$, $C_1 \cdot C_2 = 1$, $C_2 \cdot C_2 = 2$; e $A^T b = (C_1 \cdot b,\ C_2 \cdot b) = (1 + 0 + 0,\ 0 + 1 + 0) = (1, 1)$. Per risolvere si sottrae la seconda equazione dalla prima: $x - y = 0$, quindi $x = y$, e la prima diventa $3x = 1$.

Residuo: $b - A\hat X = \left(1 - \frac13,\ 1 - \frac13,\ 0 - \frac23\right) = \left(\frac23, \frac23, -\frac23\right)$ (qui $A\hat X = (x, y, x + y) = \left(\frac13, \frac13, \frac23\right)$). Controllo: è ortogonale alle colonne $(1, 0, 1)$ e $(0, 1, 1)$: $\frac23 - \frac23 = 0$ ✓. L'errore minimo è $\|b - A\hat X\|^2 = 3 \cdot \frac49 = \frac43$: nessuna scelta di $x, y$ fa meglio. Il compromesso "spartisce" l'errore in parti uguali tra le tre equazioni, come ci si aspetta per simmetria.

La retta di regressione

Per i punti $(x_1, y_1), \dots, (x_m, y_m)$ e il modello $y = \beta x + \gamma$:

$$A = \begin{pmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_m & 1 \end{pmatrix}, \qquad X = \begin{pmatrix} \beta \\ \gamma \end{pmatrix}, \qquad b = \begin{pmatrix} y_1 \\ \vdots \\ y_m \end{pmatrix}.$$

Le colonne di $A$ sono $C_1 = (x_1, \dots, x_m)$ e $C_2 = (1, \dots, 1)$, quindi $C_1 \cdot C_1 = \sum x_i^2$, $C_1 \cdot C_2 = \sum x_i$, $C_2 \cdot C_2 = 1 + \dots + 1 = m$, $C_1 \cdot b = \sum x_i y_i$, $C_2 \cdot b = \sum y_i$ (le SommatorieIl simbolo di sommatoria, le sue proprietà (linearità, additività, cambio di indice) e le somme notevoli di Gauss e geometrica.Sommatorie → vanno da $i = 1$ a $m$). Le equazioni normali, scritte con le somme, sono

$$\begin{pmatrix} \sum x_i^2 & \sum x_i \\ \sum x_i & m \end{pmatrix} \begin{pmatrix} \beta \\ \gamma \end{pmatrix} = \begin{pmatrix} \sum x_i y_i \\ \sum y_i \end{pmatrix} .$$

La seconda equazione, divisa per $m$, dice $\beta \bar x + \gamma = \bar y$ (con $\bar x$, $\bar y$ le medie): la retta passa sempre per il baricentro dei dati. Le colonne di $A$ sono indipendenti (quindi la soluzione è unica) purché non tutti gli $x_i$ siano uguali: due vettori sono dipendenti solo se uno è multiplo dell'altro, e $C_1 = c\, C_2 = (c, \dots, c)$ vuol dire proprio $x_1 = \dots = x_m = c$ (punti tutti sulla stessa verticale, per cui nessuna retta $y = \beta x + \gamma$ va bene).

Esempio. Punti $(-1, 0)$, $(0, 1)$, $(1, 3)$:

$$\sum x_i^2 = 1 + 0 + 1 = 2, \quad \sum x_i = -1 + 0 + 1 = 0, \quad m = 3, \quad \sum x_i y_i = (-1) \cdot 0 + 0 \cdot 1 + 1 \cdot 3 = 3, \quad \sum y_i = 0 + 1 + 3 = 4 .$$

$$\begin{cases} 2\beta + 0\gamma = 3 \\ 0\beta + 3\gamma = 4 \end{cases} \Longrightarrow \beta = \frac32, \quad \gamma = \frac43 .$$

Il baricentro è $(\bar x, \bar y) = \left(0, \frac43\right)$ e infatti $\frac32 \cdot 0 + \frac43 = \frac43$ ✓.

Retta: $y = \frac32 x + \frac43$. Residui $y_i - (\beta x_i + \gamma)$: $0 - (-\frac32 + \frac43) = \frac16$, $1 - \frac43 = -\frac13$, $3 - (\frac32 + \frac43) = \frac16$; somma $0$ ✓ (è la seconda equazione normale).

(Il calcolo è stato facile perché $\sum x_i = 0$: le due equazioni si separano. Con dati simmetrici rispetto a $0$ succede sempre.)

Altri modelli. Lo stesso metodo funziona per qualunque modello lineare nei parametri, per esempio una parabola $y = a x^2 + b x + c$: la matrice $A$ ha righe $(x_i^2, x_i, 1)$ e le incognite sono $a, b, c$.

Procedura

1. Scrivere il sistema $AX = b$ (una equazione per dato).
2. Calcolare $A^T A$ e $A^T b$.
3. Risolvere $A^T A\, X = A^T b$.
4. Controllo: il residuo $b - A\hat X$ deve essere ortogonale a ogni colonna di $A$.

Errori comuni

• Calcolare $AA^T$ invece di $A^TA$: è una matrice $m \times m$, non serve.
• Pensare che la soluzione ai minimi quadrati risolva il sistema: risolve il sistema normale; $A\hat X$ in generale è diverso da $b$.
• Minimizzare le distanze perpendicolari dei punti dalla retta: il metodo minimizza gli scarti verticali $y_i - (\beta x_i + \gamma)$.
• Dimenticare la colonna di $1$ nella matrice della retta di regressione (corrisponde al termine noto $\gamma$).