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 .$$

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$. Il sistema ha soluzione esattamente quando $b \in \operatorname{Im} A$.

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 .$$

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 .$$

Residuo: $b - A\hat X = \left(1 - \frac13,\ 1 - \frac13,\ 0 - \frac23\right) = \left(\frac23, \frac23, -\frac23\right)$. Controllo: è ortogonale alle colonne $(1, 0, 1)$ e $(0, 1, 1)$: $\frac23 - \frac23 = 0$ ✓. 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 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.

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

$$\sum x_i^2 = 2, \quad \sum x_i = 0, \quad m = 3, \quad \sum x_i y_i = 0 + 0 + 3 = 3, \quad \sum y_i = 4 .$$

$$\begin{cases} 2\beta + 0\gamma = 3 \\ 0\beta + 3\gamma = 4 \end{cases} \Longrightarrow \beta = \frac32, \quad \gamma = \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$).