Implementazione del problema della approssimazione ai minimi quadrati

1. Corso di Analisi Numerica A.A. 2004-2005 Implementazione del problema della approssimazione ai minimi quadrati Camillo Bosco

2. Definizione del problema (richiamo) CASO DISCRETO: siano dati m punti (x1,y1),…,(xm,ym). Si vuole determinare un polinomio p*(x) appartenente a Pn, con m>n ed m,n appartenenti a N, tale che: risulti minima rispetto ai coefficienti del polinomio. I valori w1, . . . , wm sono delle costanti positive dette pesi. A scopo didattico tali costanti assumono tutte valore 1 (peso uniforme).

3. Fitting Lineare: una istanza del problema Nel caso in cui n=1 vogliamo determinare il polinomio di primo grado p(x)=ax+b (geometricamente una retta) che costituisce la migliore approssimazione ai minimi quadrati. Vogliamo cioè, noti m punti (xi,yi) i=1…m trovare i valori di a e b che minimizzano: Come minimizzare tale funzione ? Calcolando le derivate prime di F rispetto ad a e rispetto a b e ponendole uguali a zero. Si ottiene così un sistema di due equazioni in due incognite che risolto consente di esprimere a e b (le incognite !) in funzione delle coordinate degli n punti.

4. Le soluzioni del sistema

5. Esempio 1 : fitting_lineare.m

6. Ricordiamo la teoria ! Risolvere il problema discreto ai minimi quadrati equivale a determinare la soluzione di un sistema lineare sovradeterminato (m>n) ottenuto calcolando i valori p*(xi). Qualora tale sistema Ax=b non abbia soluzione si determina il vettore soluzione Risolvere il problema discreto ai minimi quadrati equivale a determinare la soluzione di un sistema lineare sovradeterminato (m>n) ottenuto calcolando i valori p*(xi). Qualora tale sistema Ax=b non abbia soluzione si determina il vettore soluzione : Tale vettore, che minimizza la somma dei quadrati delle componenti del vettore resto R = Ax-b, è soluzione del sistema: ATAx = ATb Il problema viene quindi ricondotto alla risoluzione di un sistema lineare “classico”

7. Ripassiamo il metodo di Jacobi ! E’ un metodo iterativo per la risoluzione di un sistema lineare quadrato (m=n). Al passo k+1-esimo le componenti del vettore soluzione sono così definite: i=1…n MATLAB risolve un sistema lineare utilizzando l’operatore \ nella forma x=A\b. Tale operatore corrisponde alla funzione mldivide. Matlab utilizza un metodo diretto: il metodo di eliminazione di Gauss

8. Esempio 2 : min_quad.m

9. Definizione del problema (caso continuo) CASO CONTINUO: sia w(x) una funzione positiva, continua ed integrabile in [a,b]. Si vuole determinare p*(x) appartenente a Pn in modo tale da minimizzare: La funzione w(x) è detta funzione peso. A scopo didattico si sceglie w(x)=1 cioè peso uniforme ed unitario.

10. Una istanza del problema Nel caso in cui w(x)=1 e consideriamo lo spazio V=C[a,b], ovvero lo spazio delle funzioni continue in [a,b], vogliamo trovare i coefficienti a0,…,an del polinomio p*(x) che approssima f appartenente a V. Dobbiamo minimizzare: Come minimizzare tale funzione ? Calcolando le derivate prime di F rispetto ai valori a0,…,an e ponendole uguali a zero. Otteniamo il sistema:

11. Esempio 3 : min_quad_continuo1.m In tal caso [a,b] = [0,1] il sistema diventa: PROBLEMA: La matrice dei coefficienti di tale sistema è: i,j= 0…n Si tratta della matrice di Hilbert. E’ una matrice malcondizionata !

12. Soluzione ! Per evitare di ottenere una matrice malcondizionata e difficilmente invertibile si sceglie una base differente da quella canonica in Pn. L’insieme S={1,x,x2,…,xn}, base canonica di Pn, è costituito da polinomi non ortogonali rispetto a nessun prodotto interno. Vogliamo trasformare la base canonica in una famiglia di polinomi ortogonali. Il procedimento di ortogonalizzazione di Gram-Schmidt ci consente di generare due successioni di polinomi:

13. Ortogonalizzazione di Gram-Schmidt

14. Le due successioni La successione è costituita da polinomi ortogonali. Pertanto a norma di definizione: La successione è costituita da polinomi ortonormali. Quindi: Scegliendo tali successioni la matrice del sistema di equazioni è diagonale, quindi non più malcondizionata !

15. Esempio 4 : min_quad_continuo2.m Il polinomio di migliore approssimazione ai minimi quadrati nel caso continuo è calcolato come segue: 1) Si genera l’insieme utilizzando il procedimento di ortogonalizzazione di Gram-Schmidt 2) Si calcolano i coefficienti generalizzati di Fourier: j= 0…n 3) Si costruisce il polinomio di approssimazione:

16. Teoria dei polinomi ortogonali • -Esistono delle “famiglie” di polinomi ortogonali utilizzate in diversi ambiti della analisi numerica. • Ciascuna famiglia è caratterizzata da una proprietà principale: ciascun polinomio pn della famiglia è ortogonale a tutti i polinomi di grado minore o uguale a n. • Nel caso della approssimazione si osserva che, scegliendo w(x)=1 in [-1,1], si ottiene la famiglia dei polinomi di Legendre: -Scegliendo invece w(x)=1/sqrt(1-x2) in [-1,1], si ottiene la famiglia dei polinomi di Chebichev:

17. Esempio 5 : polinomiLegendre.m Genera i primi n polinomi di Legendre secondo la seguente ricorrenza: (Rif. Naldi-Pareschi-Russo pag. 262) L0 (x)=1 L1 (x)=x (n+1) Ln+1(x)=(2(n+1)-1)xLn(x)-nLn-1(x)

18. Esempio 6 : polinomiChebichev.m Genera i primi n polinomi di Chebichev secondo la seguente ricorrenza: T0 (x)=1 T1 (x)=x Tn+1(x)=2xTn(x)-Tn-1(x)