C05-Rezolvarea sistemelor de ecuatii_3.pdf

Cursul 5Cursul 5Metode de rezolvare a sistemelor liniare

bazate pe factorizarea ortogonal

R l i t l d ii li i Rezolvarea sistemelor de ecuaii liniare n sensul celor mai mici ptrate.

Sistemul de ecuaii liniareAx=b, ARmxn, b Rm, m>nnu admite n general soluie.Soluia n sensul celor mai mici ptrate (sauSoluia n sensul celor mai mici ptrate (sau pseudosoluia) se definete ca vectorul x*din Rn care asigur minimizarea normei geuclidiene a vectorului reziduu:

* xAbminxAb2Rx2

xAbminxAbn

=

C id ti l t i l V K (K RConsideram spatiul vectorial V = Kn(K=R sau K=C), nN*. Pe acest spatiu orice doua

t hi l tnorme sunt echivalente. Vom nota cu ||||

, ||||1, ||||2 urmatoarele

norme uzuale pe Kn:

Metode de rezolvare a sistemelor liniare bazate pe factorizare ortogonal.

Pseudosoluia este soluie exact pentru sistemul normal ATAx=ATb

Sistemul normal este ru condiionat astfel nct metodele obinuite de rezolvare (Gauss, Cholesky, etc) nu dau rezultate satisfctoare, fiind necesare metode mai t bil di t d d istabile din punct de vedere numeric,

bazate pe triunghiularizare ortogonal. P d l i t i d t i APseudosoluia este unic dac matricea A

are coloanele liniar independente, caz n care se exprim cacare se exprim ca

x*=(ATA)-1ATb=A#b


n care se defineteA#=(ATA)-1ATRnxm

ca invers generalizat sau pseudoinvers P M t i i d t hi lPenrose - Moore a matricei dreptunghiulare

A Rmxn.O t d t l t f i t lO metod ortogonal transform sistemul Ax=b

ntr-un sistem echivalent cu matrice superior triunghiulartriunghiular

HAx=Hbn care matricea de transformare HRmxn esten care matricea de transformare HRmxn este

ortogonal.


Matricea sistemului echivalent (n sensul c are aceeai pseudosoluie cuc are aceeai pseudosoluie cu sistemul iniial, deci i acelai numr de condiionare n norm 2.

Transformrile ortogonale conserv norma euclidian

22

TTTT22 xxxHxHx)Hx()Hx()Hx,Hx(Hx ====== 22 )()(),(

1xH

supH 2 =

= 1x

supH2

0x2 ==


22222AAHHARHAR ===

22

T

2

T

2

T RRHRHARHA ===

Metoda Householder22

AHA =

Metoda Householder

Reflectorii Householder sunt matrice de forma:

pTp

Tpp

mp vv

vv2IH

=

pp

2

2pp

Tp v

21

2

vv =

=

222

Reflectori Householder

Matricea Householder se folosete sub forma echivalent:H=I-2uuT, cu uTu=1

T

2

T

2T

T

uu2Ivv

vv

2Ivv

vv2IH ===

yx Teorem: x,yRn, xy, ||x||2=||y||2, 2yx

yxu

=

atunci (I-2uuT)x=y.

||x-y||22=(x-y)T(x-y)=xTx-xTy-yTx+yTy=2(xTx-yTx)

(I-2uuT)x=( )( ) ( )( )

( )xyxxyx2

xxyxyx

2x TTTT

2

T =

(I 2uu )x=

(I-2uuT)x=x-(x-y)=y

( )yyxx2xxyx2x TT22

Reflectori Householder

Consecin: =||x||2, v=x+e1, u=v/||v||2atunci (I-2uuT)x=-e1( ) 1Caz particular pentru x=-e1

uuT=vvT/||v||22

||v||22=(x+e1)(x+e1)=xTx+xTe1+e1Tx+2e1Te1 2 2 2=||x||22+2+2x1=22+2x1=2(+x1)

ntruct apare la numitor evitm anularea:=sign(x1)||x||2

Metoda Householder

Matricile reflectori elementari sunt simetrice

( ) TTTTT

ppT Hvv1

Ivv1

Ivv

IH

Matricile sunt ortogonale

( ) pppmppmmp HvvIvvIIH ===

=

( )m2

Tpp

Tpp

Tpp

Tpp

mTpp I

vvvv

vv

vvIHH =+=

Un reflector Householder Hp se aplic unui vector xRn pentru a i anula ultimele x i=p+1:m

xRn pentru a-i anula ultimele xi, i=p+1:mcomponente

( ) 0m:1pxx2

2

m

pi

2i

2 +== = 2pi=

Metoda HouseholderMetoda HouseholderFixm componentele vectorul vp care intervine n

expresia reflectorului Hp la valorile:p p

=i

1p:1i0

+==+=

m:1pix

pixv

i

pip

Acest vector poart numele de vector HouseholderHouseholder.

Pentru a evita ca vpp s se anuleze sau s aib valori mici vom alege ca , s aib acelaivalori mici vom alege ca , s aib acelai semn cu xp

Metoda Householder

( ) ( ) ( )2p

m

pi

2ip m:pxxsignxxsign ==

=

Aplicm transformarea Hp vectorului x:

pi=

pp

Tp

Tpp

mp vxv

xvxx

vvIxH =

=

=

xv

Tp =

+=== m

22m

22 xv1

v1

v1

+=== +== 1pi

ipppi

ip2pxv

2v

2v

2

Metoda Householder

( )[ ] ( ) ppp2p22p vxxx21

=+=++=2

( ) ( )xxvxv

xv

m

1pi

2ippp

m

piiipT

p

+=

=

=

+==

( ) ( )xx pp ++

( ) xxx 22++( )( ) 1x

xxx

p

ppp =+

++=

( )

===

== pivx1p:1ix

vxxH ppp

i

ipiip( )

+= m:1pi0

ppppipiip

Metoda Householder

n urma aplicrii transformrii Hp, vectorul x se modific astfel:primele p-1 componente rmn nemodificatecomponenta p devine mare n valoare absolutp prestul componentelor (p+1:m) se anuleaz.

n locul vectorului Householder vp am putea folosi p pvectorul normat

pv

2p

p

vu = .1u =

1uu,uu2IH TT ==

Metoda Householder

Particularizm pentru p=1

( ) 1i, =( )i1 m:2i,0

,xH

==

11211e)x(signxexH ==

=+

=m:2ix

xv 11i

= m:2i,xi

exv += 11 exv +=

Metoda Householder

Transformarea H1 este determinat numai de vectorul v i de .

Aadar

0

x

x

1H2

1

0x m

LL

vn care

= ,v

v

H 2

1

1 L

Cu excepia primei componente, vectorii v i x

vm

L

Cu excepia primei componente, vectorii v i xcoincid, ceea ce ne sugereaz s pstrm pe vn x, n poziiile 2:m, unde apar zerouri.

Metoda Householder

Cum vectorul v este determinat esenial numai ca direcie, l vom norma cu : vv = ,

+=+=

1x1x

v1

1

Prima component v1 va fi pstrat n ntr-

==

=m:2i,x

xv

ii

1i

Prima component v1 va fi pstrat n . ntradevr, nainte de scalarea cu , =.v1

x

12H2

1

v,v

x

x

1 =

1

mm vx

LL

Metoda Householder

ntruct transformarea Hp modific numai componentele din poziiile p:m, ea se poate realiza aplicnd H unui vector ( m)aplicnd H1 unui vector x(p:m)

Aplicm reflectorul H1 unui vector oarecare y

11

T1

T11

m1 vyvyv

yyvv

IyH =

=

=

yvm

i1iT

y

yv

1ii1iT

1==

=

( ) m:1i,vyyH 1iii1 ==

Teorema lui Householder

n locul rezolvrii sistemului A.x=b, cu ARmxn, bRm se rezolv sistemul echivalent H.A.x=H.b, cu HRmxm, H ortogonal sistem care are aceeai condiionareH-ortogonal, sistem care are aceeai condiionare cu sistemul iniial (cond(H.A)=cond(A))

Pentru orice matrice ARmxn exist o matricePentru orice matrice AR exist o matriceortogonal HRmxm astfel nct H.A=R n care R Rmxn este o matrice superior triunghiularR R este o matrice superior triunghiular

Matricea ortogonal se formeaz dintr-un produs

121nn HHHHH = L

AHHHHAH 121nn = K


cu p=1:n, pornind cu A1=A. De exemplu reflectorul H1 care anuleaz prima

pp1p AHA =+

coloan a matricei ARmxn este:[beta, A(:,1)]=HSH1x(A(:,1))for j=2:n

A(:,j)=HSH1y(A(:,1), beta, A(:,j));end

Matricea are structura

= 1pp V0

0IH

+ 1pnp V0


Spre deosebire de metodele gaussiene, metodele ortogonale asigur elemente diagonale mari n valoare absolut, ceea ce confer o ,stabilitate deosebit a metodei.

Vom determina matricea Hp, care transform matricea A utiliznd algoritmul HSHx()matricea Ap, utiliznd algoritmul HSHx(),punnd ca vector x coloana p din HpAp

[ ] [ ]aHaHaHaaaHAHColoanele situate n dreapta coloanei p, (j>p) se modific folosind algoritmul HSHy()

[ ] [ ]nppp1pnp1ppp aHaHaHaaaHAH == LLLL

(j>p) se modific folosind algoritmul HSHy()Transformarea Hp nu modific coloanele aj

pentru care =0 (coloanele din stnga coloanei )p, cu j


Matricea H poate fi pstrat n triunghiul inferior din A (care se anuleaz), iar ( )elementele diagonale n vectorul beta. n triunghiul superior din A vom avea matricea R

In rezumat, transformarea aplicat matricei :i las neschimbate primele coloane pmodific elementele din coloana astfel:

primele elemente rmn neschimbatprimele elemente rmn neschimbatelementul diagonal devine mareelementele subdiagonale se anuleazg


modific elementele din coloanele j>p astfel:primele p-1 elemente rmn neschimbateelementele din poziiile p:m se calculeaz cu

relaia Householder. Pornind de la sistemul iniial, prin aplicarea

transformrilor:

n:1p,AA,AHA 1pp1p ===+

bb,bHb 1pp1p ==+

se ajunge la sistemul echivalent superior triunghiular

1n1n bxA ++ =


Matricile factori i din factorizarea a matricei A:

AR 1nAR +=

( ) mT121nnT IHHHHHQ K==

adic transformrile se aplic matricei unitate i

( ) m121nnQ K

adic transformrile se aplic matricei unitate i se transpune rezultatul, sau se calculeaz direct:

n21Tn

T2

T1

T HHHHHHHQ KK ===

Rotatori Givens (matrice de rotaie elementare)

1k

I

sc

I

G

=

lm

1klkl

I

cs

IG

c2+s2=1, definete o rotaie de ordinul m n planul(k l) hi l d i i

lm

(k,l) cu unghiul , unde c=cos i s=sin

Matricea Gkl este ortogonal: Gkl*GTkl=InTransformarea Gkl aplicat unui vector xRm i

modific numai componentele xk i xlGklx=[x1cxk+sxl-sxk+cxlxm]

Rotatori Givens (matrice de rotaie elementare)

function x=rotvec(k, l, c, s, x)%aplica transformarea Gkl vectorului xt = zeros(2,1);t = [c s;-s c]*[x(k); x(l)];(k) t(l)x(k) = t(l);x(l) = t(2);

Determinm transformarea (c s) care anuleazDeterminm transformarea (c, s), care anuleaz componenta l :

-s x +c x =0-s.xk+c.xl=0

rx

xx

xc k

22

k =+

=rxx lk +

rx

xx

xs l

22

l =+

=rxx lk +

Rotatori Givens (matrice de rotaie elementare)22

rxxxx

xxxsxcx 2l

2k2

l2k

2l

2k

lkk =+=+

+=+=

Rotaia Gkl, care modific componentele xk=r i xl=0 este:

f nction [c s ] rot(k l )function [c, s, x]=rot(k, l, x)% determina rotatia (c,s) care% anuleaza x(l)% anuleaza x(l)r=sqrt(x(k)^2+x(l)^2);c=x(k)/r;s=x(l)/r;

n metoda Givens matricea ortogonal G se formeaz ca un produs de matriceformeaz ca un produs de matrice elementare de "rotaie" de forma

Rotatori Givens (matrice de rotaie elementare).GGGGGGG

12n1n G

121n,1n1

G

1n,2nn,2n

G

n,1n 444 3444 21L

444 3444 21321 =

Matricea Gkl afecteaz prin nmulire numai liniile k i l:

,1n:1kAGA )p(kl)1p( ==+

.n:1klAA

,1n:1kAGA)0(

kl

+==

Plecnd cu matricea A i folosind transformrile ortogonale G se obinetransformrile ortogonale Gkl se obine o matrice superior triunghiular:

Factorizarea QROrice matrice ARmxn poate fi descompus sub

forma A=Q1R1, cunoscut sub numele de factorizare QR n care Q = este o matrice cufactorizare QR n care Q1= este o matrice cu coloane ortogonale i R1 este o matrice ptrat superior triunghiular.

[ ] ,RQ0

RQQQ 11

121 =

=

( ) .RR,RQ,RQ nn1nnm

2nm

1 ,, 121

( ) n21TnT2T1T11nnT HHHHHHHHHHQ KKK ==== [ ]

=

=

=

0

IHHH

0

IQ

0

IQQQ nn21

nn211 K

000

Factorizarea QR

for i = n:-1:1 s = A(i,i+1:n) * x(i+1:n);x(i) = (b(i) - s) / A(i,i);

endQ Q'R = A; Q = Q';

n MATLAB apelul [Q,R] = qr(A) produce descompunerea matricei A ntr-o matrice superiordescompunerea matricei A ntr-o matrice superior triunghiular RRmxn (de aceeai dimensiune cu A) i o matrice QRmxm unitar astfel nct A=Q*R.)

[Q,R] = qr(A,0) produce o descompunere "economic", n care se calculeaz numai primele n coloane ale matricei Q.

Ortogonalizarea Gram-Schmidt

[ ] [ ]

n1j111

rr00

0

rrr

qqqaaa

LLLL

LL

[ ] [ ]

=

nn

jnjjnj1nj1

r000

rr00qqqaaa

L

LLLLL

LLLLL

Coloanele matricei Q sunt ortogonale (formeaz o baz ortonormat):

=

=ji0

ji,1qq j

Ti

ji,0j

1111 rqa =de unde i 1q1 = 111 ar = 1111 r/aq =

Ortogonalizarea Gram-Schmidt+ Ta2=q1r12+q2r22, q1Ta2=r12,

||q2||=1, r22=||a2-q1r12||, ( )/q2=(a2-q1r12)/r22

a3=q1r13+q2r23+q3r33, q2Ta3=r23,|| || 1 || ||||q3||=1, r33=||a3-q1r13-q2r23||, q3=(a3-q1r13-q2r23)/r33aj=q1r1j+q2r2j++qjrjjqkTaj=rkj, k=1:j-1

1j

||qj||=1, =

=1k

kjkjjj rqar

1j

jj1k

kjkjj r/rqaq

=

=

Algoritmul Gram-Schmidt clasicpentru j=1:n

qj = ajj jpentru k=1:j-1

rkj = qk*qjj j

=1j

kjkjj rqqq

rjj=||qj||

qj = qj/rjj

=1k

qj qj/ jjAlgoritmul Gram-Schmidt clasic este numeric

instabil. Rezultate mult mai bune ne furnizeaz o variant modificat

Algoritmul Gram-Schmidt modificatq1ak=r1kq1 ak r1kq1q1ak=q1r1k(I-q1q1)ak=q1r2k++qkrkkq2(I-q1q1)ak=r2kq2q2(I-q1q1)ak=q2r2k(I )(I )(I-q1q1)(I-q2q2)ak=q3r3k++qkrkk

Q = AQ = Apentru k=1:nrkk=||qk||kk ||q ||qk=qk/rkkpentru j=k+1:n

r =q qrkj=qkqjqj=qj-qkrkj

Algoritmul Gram-Schmidt clasic

[Q, R] = Gram_Schmidt(m, n, A)% Intrri: A matricea de factori at (ba a iniial)A = matricea de factorizat (baza iniial)% Ieiri:Q = factorul ortogonal(baza% ortonormat)% ortonormat)% R = factorul superior triunghiular

[m,n]=size(A);for i = 1 : n

R(1:i-1,i) = Q(1:m,1:i-1)'*A(1:m,i);A(1 i) Q(1 1 i 1)*R(1 i 1 i)y = A(1:m,i)-Q(1:m,1:i-1)*R(1:i-1,i);

R(i,i) = norm(y);Q(1:m,i) = y ./ R(i,i);Q(1:m,i) y ./ R(i,i);

end

Algoritmul Gram-Schmidt modificat[Q R] G S h idt M difi t( A)[Q, R] = Gram_Schmidt_Modificat(m, n, A) %Intrri: A = matricea de factorizat% (baza iniial)% (baza iniial)%Ieiri:Q = baza ortonormat % R = factorul superior triunghiular[m,n]=size(A);for i = 1 : n R(i i) =norm(A(1:m i));R(i,i) norm(A(1:m,i));Q(1:m,i) = A(1:m,i) / R(i,i);for j = i+1 : n R(i,j) = Q(1:m,i)' *A(1:m,j);A(1:m,j)= A(1:m,j)- Q(1:m,i)*R(i,j);

endendend

C05-Rezolvarea sistemelor de ecuatii_3.pdf

Documents

Transcript of C05-Rezolvarea sistemelor de ecuatii_3.pdf