Approssimazione ai minimi quadrati

Corso di Analisi Numerica A.A. 2004-2005
Implementazione del problema della
approssimazione ai minimi quadrati
Camillo Bosco
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:
m
 w  p *(x )  y 
i 1
2
i
i
i
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).
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:
m
F ( a, b)    yi  ( axi  b) 
2
i 1
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.
Le soluzioni del sistema
a
n
n
n
i 1
i 1
i 1
n  xi yi  (  xi )(  yi )
n
n
n  xi  (  xi )
2
i 1
b
2
i 1
n
n
n
n
i 1
i 1
n
i 1
i 1
(  xi 2 )(  yi )  (  xi )(  xi yi )
n
n (  xi 2 )  (  xi ) 2
i 1
i 1
Esempio 1 : fitting_lineare.m
xi
yi
0
10
1
25
2
51
3
66
4
97
5
118
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 x  Rn :
Ax  b  Ax  b 2 x  Rn
2
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”
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:
( k 1)
xi
i 1
 (bi   aij x j 
j 1
(k )
n
a x
j i 1
ij
(k )
j
) / aii
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
Esempio 2 : min_quad.m
xi
yi
1
2.1
2
2.9
3
3.2
0
1
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:

b
a
w( x)  p *( x)  f ( x) dx
2
La funzione w(x) è detta funzione peso.
A scopo didattico si sceglie w(x)=1 cioè peso uniforme ed
unitario.
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:
b
n
a
j 0
F (a0 ,..., an )   ( f ( x)   a j x j )2dx
Come minimizzare tale funzione ? Calcolando le derivate prime
di F rispetto ai valori a0,…,an e ponendole uguali a zero.
Otteniamo il sistema:
n
b
j 0
a
i
a
x
dx

f
(
x
)
x
dx
 j

i j
b
a
Esempio 3 : min_quad_continuo1.m
In tal caso [a,b] = [0,1] il sistema diventa:
n
aj
 i  j 1
j 0
1
  f ( x ) x i dx
0
PROBLEMA:
La matrice dei coefficienti di tale sistema è:
 1 
H 

i

j

1


Si tratta della
malcondizionata !
matrice
di
Hilbert.
i,j= 0…n
E’
una
matrice
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:
qk ( x )k 0
n
 pk ( x )k 0
n
Ortogonalizzazione di Gram-Schmidt
q0 ( x)  1
q0 ( x )
p0 ( x ) 
q0
q0 
k 1
qk ( x )  x   x , p j p j ( x )
k
k
j 0

b
a
q0 2 ( x )dx
qk ( x )
pk ( x ) 
qk
Le due successioni
La successione
qk ( x )k 0 è costituita da polinomi ortogonali.
n
Pertanto a norma di definizione:
q j , qk   w( x)q j ( x)qk ( x)dx  00
b
a
per k  j
per k=j
La successione  pk ( x )k 0 è costituita da polinomi ortonormali.
Quindi:
n
p j , pk   w( x) p j ( x) pk ( x)dx  10
b
a
per k  j
per k=j
Scegliendo tali successioni la matrice del sistema di equazioni
è diagonale, quindi non più malcondizionata !
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
 pk ( x )k 0
n
utilizzando il procedimento
di ortogonalizzazione di Gram-Schmidt
2) Si calcolano i coefficienti generalizzati di Fourier:
f , pj
j= 0…n
3) Si costruisce il polinomio di approssimazione:
n
pn* ( x )   f , p j p j ( x )
j 0
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:

1
1
L j ( x) Lk ( x)dx  0 per j  k
-Scegliendo invece w(x)=1/sqrt(1-x2) in [-1,1], si ottiene la
famiglia dei polinomi di Chebichev:
1
T j ( x )Tk ( x )
1
1 x

2
dx  0
per j  k
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)
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)