Calcolo numerico/Approssimazione media dei minimi quadrati

Descrizione del problema

Supponiamo di avere le coppie di dati $(x_{k},y_{k})_{k=0}^{m}$ con nodi $x_{k}$ a due a due distinti, e numerosi. Cerco

$P_{n}^{*}(x)=\sum _{j=1}^{n}a_{j}^{*}\phi _{j}(x),$

dove $\{\phi _{j}(x)\}_{j=0}^{n}$ è una base, e

$\sum _{k=0}^{m}(y_{k}-p_{n}^{*}(x_{k}))^{2}\leq \sum _{k=0}^{m}(y_{k}-p_{n}(x_{k}))^{2}$

per ogni $p_{n}\in P_{N}$, cioè $p_{n}^{*}$ è il polinomio che minimizza lo scarto quadratico medio. Inoltre, introducendo nella sommatoria pesi $w_{k}$ positivi, in modo da dare più importanza a determinati addendi rispetto ad altri, si ha

$\sum _{k=0}^{m}\omega _{k}(y_{k}-p_{n}^{*}(x_{k}))^{2}\leq \sum _{k=0}^{m}(y_{k}-p_{n}(x_{k}))^{2}\quad {\hbox{formula 1}}$

lo scopo è quello di minimizzare lo scarto quadratico medio pesato.

Determinazione del sistema

Consideriamo per semplicità la base canonica dell'insieme dei polinomi, quindi $\phi _{j}=x^{j}$. Consideriamo

$p_{n}(x)=\sum _{j=0}^{n}a_{j}\phi _{j}(x)$

e inserendo l'espressione di $p_{n}$ nel lato destro della formula 1:

$\sum _{k=0}^{m}\omega _{k}(y_{k}-\sum _{j=0}^{n}a_{j}\phi _{j}(x_{k}))^{2}=\Phi (a_{0},\dots ,a_{n})$

che è un funzionale nei parametri da determinare. Cerchiamo il valore dei parametri per cui otteniamo la minimizzazione di $p_{n}$.

Il minimo rende stazionario il gradiente di $\Phi$, quindi calcoliamo le derivate parziali. Per ogni $i=0,\dots ,n$, si ha

${\frac {\partial \Phi }{\partial a_{i}}}=\sum _{k=0}^{m}\omega _{k}*2*(y_{k}-\sum _{j=0}^{n}a_{j}\phi _{j}(x_{k}))*(-\phi _{i}(x_{k}))=0$

$\sum _{k=0}^{m}\omega _{k}\sum _{j=0}^{n}a_{j}\phi _{j}(x_{k})*\phi _{i}(x_{k})=\sum _{k=0}^{m}w_{k}y_{k}\phi _{i}(x_{k})$

e invertendo le sommatorie

$\sum _{j=0}^{n}(\sum _{k=0}^{m}\omega _{k}\phi _{j}(x_{k})*\phi _{i}(x_{k}))a_{j}=\sum _{k=0}^{m}\omega _{k}y_{k}\phi _{i}(x_{k})$

Abbiamo un sistema lineare di $n+1$ equazioni nelle $n+1$ incognite a $a_{0},\dots ,a_{n}$. Riscrivendo il tutto in maniera compatta abbiamo il sistema:

$C{\boldsymbol {a}}={\boldsymbol {b}}$

dove

$C=(C_{ij})=\sum _{k=0}^{m}\omega _{k}\phi _{j}(x_{k})\phi _{i}(x_{k}),$

mentre

$b_{i}=\sum _{k=0}^{m}\omega _{k}y_{k}\phi _{i}(x_{k}),\quad i,j=0,\dots n$

Per dimostrare che la soluzione di tale sistema è un minimo, analizziamo la matrice $C=H_{\Phi }$. Dobbiamo verificare che l'hessiana è definita positiva e simmetrica.

Proposizione 5.1

La matrice $C=H_{\Phi }$ è definita positiva (e di conseguenza esiste unica la soluzione del problema, e a ${\boldsymbol {a}}$ è un minimo per $\Phi$.)

Dimostrazione

Fattorizziamo la matrice $C$ per decomporla in pezzi, nella forma

$C=P^{t}WP$

$P^{t}$ è una matrice rettangolare di $n+1$ righe e $m+1$ colonne, con $m>>n$, $W$ è una matrice quadrata $m+1\times m+1$, mentre $P$ è $(m+1)\times (n+1)$.

$P$ è definita nel modo seguente.

$P_{st}=\phi _{t}(x_{s})$

$W={\rm {{diag}(w_{1},\dots ,w_{m+1})}}$

La matrice della forma $P^{t}WP$ è detta matrice delle equazioni normali pesate. Il termine noto si può scrivere come $P^{t}WY$.

Dimostriamo che $C$ è simmetrica.

$C^{t}=(P^{t}WP)^{t}=P^{t}W^{t}P=P^{t}WP$

perché $W$ è diagonale.

Dimostriamo che $C$ è definita positiva.

Per ogni ${\boldsymbol {x}}\neq 0$, dobbiamo dimostrare che

${\boldsymbol {x}}^{t}C{\boldsymbol {x}}>0$

${\boldsymbol {x}}^{t}C{\boldsymbol {x}}={\boldsymbol {x}}^{t}P^{t}WP{\boldsymbol {x}}$

e definendo ${\boldsymbol {y}}=P{\boldsymbol {x}}$ si ha

${\boldsymbol {x}}^{t}C{\boldsymbol {x}}={\boldsymbol {y}}^{t}W{\boldsymbol {y}}$

Supponiamo che $P$ abbia rango massimo. Allora ${\boldsymbol {y}}\neq 0$, e siccome $W$ è definita positiva per definizione, ${\boldsymbol {y}}^{t}W{\boldsymbol {y}}>0$.

Dimostriamo che $P$ ha rango massimo.

Caso particolare: consideriamo $\phi _{j}(x)=x^{j}$. Allora

$P_{s,t}=(x_{k}^{j}),k=0,\dots ,m,j=0,\dots ,n$

Per esteso

$P=\left({\begin{array}{ccccc}1&x_{0}&x_{0}^{2}&\dots &x_{0}^{n}\\1&x_{1}&x_{1}^{2}&\dots &x_{1}^{n}\\\dots &\dots &\dots &\dots &\dots \\1&x_{m}&x_{m}^{2}&\dots &x_{m}^{n}\end{array}}\right)$

Questa è la sottomatrice principale della matrice di Vandermonde. Siccome i nodi sono a due a due distinti, la Vandermonde è non singolare, allora ogni sua sottomatrice ha colonne linearmente indipendenti. $P$ ha rango massimo.

Considerando per lo spazio $P_{N}$ una qualsiasi altra base invece di quella canonica, vale la medesima proprietà.

Infatti la matrice del cambiamento di base è non singolare e quindi il rango è conservato.

Nel caso della base canonica dei polinomi,

$C=\sum _{k=0}^{m}\omega _{k}x_{k}^{i+j}$

In particolare:

$C_{00}=\sum _{k=0}^{m}\omega _{k}$

$C_{10}=\sum _{k=0}^{n}\omega _{k}x_{k}=C_{01}$

$C_{20}=\sum _{k=0}^{m}\omega _{k}x_{k}^{2}=C_{11}=C_{02}$

Abbiamo una matrice con elementi costanti lungo le antidiagonali. Una matrice con questa struttura viene chiamata matrice di Henkel.

Ricerca della soluzione del sistema

Data la scrittura $C=P^{t}WP$, l'approssimazione ai minimi quadrati, come impostato inizialmente, consiste nella risoluzione del sistema sovradeterminato $P{\boldsymbol {a}}={\boldsymbol {y}}$ di $m+1$ equazioni in $n+1$ incognite. Infatti, per minimizzare la norma del residuo, risolviamo

$\sum _{j=0}^{n}a_{j}\phi _{j}(x_{k})=y_{k},\quad k=0,\dots ,m$

Non c'è una soluzione classica di questo sistema, ma dobbiamo cercare una soluzione più debole.

Per minimizzare la norma euclidea del residuo, vogliamo determinare il vettore ${\boldsymbol {a}}^{*}$ tale che

$\phi ({\boldsymbol {a}}^{*})=(\vert {\boldsymbol {y}}-P{\boldsymbol {a}}\vert )^{2}=\min _{{\boldsymbol {a}}\in \mathbb {R} ^{n+1}}(\vert {\boldsymbol {y}}-P{\boldsymbol {a}}\vert )^{2}$

Invece di determinare la soluzione di $P{\boldsymbol {a}}={\boldsymbol {y}}$ cerchiamo la soluzione di

$P^{t}P{\boldsymbol {a}}=P^{t}{\boldsymbol {y}}={\boldsymbol {b}}$

(sistema delle equazioni normali). Se supponiamo di avere pesi $\omega _{k}$ generici, possiamo definire

$w(R)=R^{t}WR=\sum _{k}\omega _{k}r_{k}^{2}$

Possiamo pensare questa norma come la norma euclidea al quadrato di $W^{1/2}R$.

Possibilità per la risoluzione del sistema dei minimi quadrati

Prima possibilità: consideriamo la fattorizzazione di Cholesky, per matrici simmetriche e definite positive come $C$.

Risolviamo quindi i sistemi

$R^{t}{\boldsymbol {z}}={\boldsymbol {b}}$

$R{\boldsymbol {a}}={\boldsymbol {z}}$

Consideriamo per semplicità $w_{k}=1$.

Analizziamo il costo computazionale nella costruzione di $C$:

$C_{ij}=\sum _{k=0}^{m}\omega _{k}\phi _{j}(x_{k})\phi _{i}(x_{k})$

$C$ ha $n^{2}$ elementi, e ognuno richiede $m$ prodotti. $C$ è simmetrica quindi basta calcolare metà degli elementi, per un totale di $n^{2}*m/2$ operazioni.

Siccome il costo della fattorizzazione di Cholesky di $C$ è $n^{3}/6$, in totale abbiamo: $n^{2}/2*(m+n/3)$ operazioni. Quando i dati sono molti ($m$ grande) si hanno problemi di accuratezza e di rappresentabilità nella fattorizzazione di Cholesky.

Fattorizzazione QR per la risoluzione del sistema dei minimi quadrati

Supponiamo di avere $A\in C_{m\times n}$ di rango massimo (nel problema corrisponde a $P$). Calcoliamo la fattorizzazione $QR$ di $A$.

Considerata $Q_{h}=p_{1}p_{2}\dots p_{n}$ (nel caso di matrici rettangolari abbiamo $n$ colonne da annullare), otteniamo una matrice $R$ con $m$ righe e $n$ colonne, suddivisa in due blocchi: $R_{1}$ triangolare superiore e $r_{2}$ blocco trapezoidale di zeri delle sottocolonne annullate.

Il blocco superiore è sicuramente non singolare, perché $A$ ha rango massimo.

$\vert A{\boldsymbol {x}}-{\boldsymbol {b}}\vert _{r}=\vert QR{\boldsymbol {x}}-{\boldsymbol {b}}\vert$

$=\vert Q(R{\boldsymbol {x}}-QH{\boldsymbol {b}})\vert$

ma siccome $Q$ è unitaria conserva la norma 2:

$=\vert R{\boldsymbol {x}}-Q_{h}{\boldsymbol {b}}\vert$

e ponendo $Q_{h}{\boldsymbol {b}}={\boldsymbol {c}}$, siccome $R$ ha due blocchi e il secondo è nullo, e siccome ${\boldsymbol {c}}=({\boldsymbol {c}}_{1},{\boldsymbol {c}}_{2})$, si ha

$\vert R_{1}{\boldsymbol {x}}-{\boldsymbol {c}}_{1}-{\boldsymbol {c}}_{2}\vert _{2}$

Cerchiamo

$\min _{x\in \mathbb {C} ^{n}}\vert A{\boldsymbol {x}}-{\boldsymbol {b}}\vert ^{2}$

che per quanto appena osservato è

$=\min \vert R_{1}{\boldsymbol {x}}-{\boldsymbol {c}}\vert ^{2}$

$(\vert R_{1}{\boldsymbol {x}}-{\boldsymbol {c}}_{1}-{\boldsymbol {c}}_{2}\vert )^{2}\leq \vert R_{1}{\boldsymbol {x}}-{\boldsymbol {c}}_{1}\vert ^{2}+\vert {\boldsymbol {c}}_{2}\vert _{2}^{2}$

$\min _{{\boldsymbol {x}}\in \mathbb {C} ^{n}}(\vert R_{1}{\boldsymbol {x}}-{\boldsymbol {c}}_{1}\vert )^{2}+\vert {\boldsymbol {c}}_{2}\vert ^{2}$

Il primo addendo è minimo quando

$R_{1}{\boldsymbol {x}}={\boldsymbol {c}}_{1}$

e rimane

$\min \vert {\boldsymbol {c}}_{2}\vert ^{2}$

Si conclude che, data la matrice $A$, si calcola la fattorizzazione $QR$ e si risolve quindi

$R_{1}{\boldsymbol {x}}={\boldsymbol {c}}_{1}$

Il costo è $n^{3}/3$ per la fattorizzazione e $mn^{2}$ per la risoluzione del sistema lineare, e in totale

$n^{2}*(m-n/3)$

che è maggiore di quello della fattorizzazione di Cholesky. Questo metodo però è più sicuro della fattorizzazione di Cholesky.

Relazione tra i due metodi

Tornando al problema di approssimazione, ci si chiede qual è la relazione tra $R_{1}$ che viene dall'applicazione della $QR$ ad $A$ generica e ${\bar {R}}$ che otteniamo facendo la Cholesky. Queste matrici sono uguali a meno di matrici di fase (diagonali unitarie).

$P^{t}WP={\bar {R}}_{h}{\bar {R}}=R^{h}Q_{h}QR=R_{1}^{h}R_{1}$

e $R_{1}$ differisce per una matrice di fase da ${\bar {R}}$.

Nel caso di pesi diversi da 1, si deve considerare la norma euclidea pesata e si può ripetere il procedimento analogo.

Consideriamo la norma 2 di

$W^{1/2}(A{\boldsymbol {x}}-{\boldsymbol {b}})$

quindi consideriamo $QR$ di ${\bar {A}}=W^{1/2}*A$.

Condizionamento

Applicando la $QR$ ad $A$, $K(R_{1})=K(R)=K(A)$.

Costruendo $C=A_{h}A$ con la Cholesky, si ha $K(A_{h}A)=K(A)^{2}$.

Quindi le equazioni normali hanno lo svantaggio di un indice di condizionamento quadrato rispetto a quello di $A$.