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^{\ast }(x)=\sum _{j=1}^n{a}_j^{\ast }{\varphi }_j(x),$

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

$\sum _{k=0}^m(y_k-p_n^{\ast }(x_k){)}^2\le \sum _{k=0}^m(y_k-p_n(x_k){)}^2$

per ogni $p_n\in P_N$, cioè $p_n^{\ast }$ è 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^{\ast }(x_k){)}^2\le \sum _{k=0}^m(y_k-p_n(x_k){)}^2\text{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 ${\varphi }_j=x^j$. Consideriamo

$p_n(x)=\sum _{j=0}^n{a}_j{\varphi }_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{\varphi }_j(x_k){)}^2=\mathrm{\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 $\mathrm{\Phi }$, quindi calcoliamo le derivate parziali. Per ogni $i=0,\dots ,n$, si ha

$\frac{\partial \mathrm{\Phi }}{\partial a_i}=\sum _{k=0}^m{\omega }_k\ast 2\ast (y_k-\sum _{j=0}^n{a}_j{\varphi }_j(x_k))\ast (-{\varphi }_i(x_k))=0$

$\sum _{k=0}^m{\omega }_k\sum _{j=0}^n{a}_j{\varphi }_j(x_k)\ast {\varphi }_i(x_k)=\sum _{k=0}^m{w}_k{y}_k{\varphi }_i(x_k)$

e invertendo le sommatorie

$\sum _{j=0}^n(\sum _{k=0}^m{\omega }_k{\varphi }_j(x_k)\ast {\varphi }_i(x_k))a_j=\sum _{k=0}^m{\omega }_k{y}_k{\varphi }_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𝒂=𝒃$

dove

$C=(C_{ij})=\sum _{k=0}^m{\omega }_k{\varphi }_j(x_k){\varphi }_i(x_k),$

mentre

$b_i=\sum _{k=0}^m{\omega }_k{y}_k{\varphi }_i(x_k),i,j=0,\dots n$

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

Proposizione 5.1

La matrice $C=H_{\mathrm{\Phi }}$ è definita positiva (e di conseguenza esiste unica la soluzione del problema, e a $𝒂$ è un minimo per $\mathrm{\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}={\varphi }_t(x_s)$

$W=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathrm{(}{\mathrm{w}}_{\mathrm{1}},\mathrm{\dots },{\mathrm{w}}_{\mathrm{m}\mathrm{+}\mathrm{1}}\mathrm{)}$

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 $𝒙\ne 0$, dobbiamo dimostrare che

${𝒙}^{t}C𝒙>0$

${𝒙}^{t}C𝒙={𝒙}^t{P}^{t}WP𝒙$

e definendo $𝒚=P𝒙$ si ha

${𝒙}^{t}C𝒙={𝒚}^{t}W𝒚$

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

Dimostriamo che $P$ ha rango massimo.

Caso particolare: consideriamo ${\varphi }_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𝒂=𝒚$ di $m+1$ equazioni in $n+1$ incognite. Infatti, per minimizzare la norma del residuo, risolviamo

$\sum _{j=0}^n{a}_j{\varphi }_j(x_k)=y_k,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 ${𝒂}^{\ast }$ tale che

$\varphi ({𝒂}^{\ast })=(|𝒚-P𝒂|{)}^2=\underset{𝒂\in {ℝ}^{n+1}}{min}(|𝒚-P𝒂|{)}^2$

Invece di determinare la soluzione di $P𝒂=𝒚$ cerchiamo la soluzione di

$P^{t}P𝒂=P^t𝒚=𝒃$

(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𝒛=𝒃$

$R𝒂=𝒛$

Consideriamo per semplicità $w_k=1$.

Analizziamo il costo computazionale nella costruzione di $C$:

$C_{ij}=\sum _{k=0}^m{\omega }_k{\varphi }_j(x_k){\varphi }_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\ast m/2$ operazioni.

Siccome il costo della fattorizzazione di Cholesky di $C$ è $n^3/6$, in totale abbiamo: $n^2/2\ast (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.

$|A𝒙-𝒃{|}_r=|QR𝒙-𝒃|$

$=|Q(R𝒙-QH𝒃)|$

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

$=|R𝒙-Q_h𝒃|$

e ponendo $Q_h𝒃=𝒄$, siccome $R$ ha due blocchi e il secondo è nullo, e siccome $𝒄=({𝒄}_1,{𝒄}_2)$, si ha

$|R_1𝒙-{𝒄}_1-{𝒄}_2{|}_2$

Cerchiamo

$\underset{x\in {ℂ}^n}{min}|A𝒙-𝒃{|}^2$

che per quanto appena osservato è

$=min|R_1𝒙-𝒄{|}^2$

$(|R_1𝒙-{𝒄}_1-{𝒄}_2|{)}^2\le |R_1𝒙-{𝒄}_1{|}^2+|{𝒄}_2{|}_2^2$

$\underset{𝒙\in {ℂ}^n}{min}(|R_1𝒙-{𝒄}_1|{)}^2+|{𝒄}_2{|}^2$

Il primo addendo è minimo quando

$R_1𝒙={𝒄}_1$

e rimane

$min|{𝒄}_2{|}^2$

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

$R_1𝒙={𝒄}_1$

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

$n^2\ast (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 $\overline{R}$ che otteniamo facendo la Cholesky. Queste matrici sono uguali a meno di matrici di fase (diagonali unitarie).

$P^{t}WP={\overline{R}}_h\overline{R}=R^h{Q}_{h}QR=R_1^h{R}_1$

e $R_1$ differisce per una matrice di fase da $\overline{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𝒙-𝒃)$

quindi consideriamo $QR$ di $\overline{A}=W^{1/2}\ast 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$.