Vectori Proprii Si Valori Proprii

Seminar 4

Calculul valorilor si vectorilorproprii

Valorile proprii ale unei matrice joaca un rol important ın multe domenii cum ar fi teoriastabilitatii. De asemenea, ele sunt legate ındeaproape de problema rezolvarii ecuatiilor alge-brice. Prin urmare, un calcul corect al valorilor si vectorilor proprii este de mare importanta.Din cauza faptului ca valorile proprii sunt radacinile ecuatiei polinomiale caracteristice, cal-culul lor pentru o matrice de ordin mai mare decat patru nu poate fi facut exact printr-osecventa finita de operatii aritmetice si este nevoie de un proces infinit care trebuie sa fietrunchiat astfel ıncat sa se obtina o aproximare acceptabila. Cel mai bun proces de acesttip este algoritmul QR care calculeaza o forma Schur aproximativa a matricei date. Acestseminar este dedicat diferitelor aspecte numerice ale calculului valorilor proprii si ale imple-mentarilor algoritmului QR.

4.1 Preliminarii

Fie A ∈ IRn×n. Un vector x ∈ ICn este un vector propriu al lui A si λ ∈ IC este valoareaproprie asociata daca sunt satisfacute urmatoarele conditii

a) x 6= 0b) Ax = λx.

Astfel, numarul λ ∈ C este o valoare proprie a lui A daca si numai daca matricea λI − Aeste singulara, i.e. valorile proprii sunt zerourile polinomului caracteristic

p(λ) = det(λIn −A),

si matricea A are exact n valori proprii, incluzand multiplicitatile. In mod evident, valorileproprii complexe ale unei matrice reale apar sub forma unor perechi conjugate. Multimeavalorilor proprii ale unei matrice A se noteaza cu λ(A) si se numeste spectrul (valorilorproprii) ale lui A.

Daca x este un vector propriu al lui A, vectorul y = αx este de asemenea un vectorpropriu al lui A, asociat aceleiasi valori proprii pentru orice scalar α diferit de zero. Decivectorii proprii sunt determinati doar prin directia lor.

1

2 SEMINAR 4. DESCOMPUNEREA VALORILOR PROPRII

Valorile proprii sunt conservate de asa numitele transformari de asemanare definite prin

B = TAT−1,

unde T este o matrice nesingulara. Daca T este o matrice ortogonala atunci A si B se numescortogonal asemenea. Daca A,B sunt asemenea atunci ele au acelasi spectru de valori propriiλ(A) = λ(B) iar vectorii proprii corespunzatori sunt legati prin relatia xB = TxA.

Cea mai simpla forma care poate fi obtinuta printr-o transformare de asemanare esteasa numita forma canonica Jordan. Daca o matrice n × n are n vectori proprii liniariindependenti, atunci forma sa canonica Jordan este diagonala si matricea se numeste simplasau diagonalizabila. Din pacate, forma canonica Jordan este foarte sensibila la perturbarinumerice. De aceea, ın calculul numeric este preferata o forma mai robusta si anume formaSchur (reala) (FSR) care poate fi obtinuta printr-o transformarea ortogonala de asemanare:

S = QT AQ =

S11 S12 · · · S1q

0 S22 · · · S2q

· · · · · · · · · · · ·0 0 . . . Sqq

,

unde Sii ∈ R1×1 or Sii ∈ R2×2 si toate blocurile diagonale 2× 2 au valori proprii complexe.Valorile reale ale unei matrice FSR pot fi determinate printr-o simpla inspectie a diagonaleiprincipale iar cele complexe pot fi calculate prin rezolvarea ecuatiilor algebrice de gradul doicorespunzatoare det(λI − Sii) = 0.

Fie o pereche (valoare proprie, vector propriu)=(λ, x), (un sistem propriu) al unei matriceA. Daca unul dintre cele doua elemente este cunoscut celalalt poate fi calculat cu usurinta(vezi problema rezolvata 4.1). Daca nici valoarea proprie nici vectorul propriu nu suntcunoscute, problema calculului lor devine una foarte dificila. Motivul sta ın imposibilitatearezolvarii ecuatiilor algebrice de grad mai mare decat patru prin secvente finite de operatiielementare.

Cele mai folosite metode pentru calculul iterativ al vectorilor proprii (fara a se cunoastevaloarea proprie asociata) sunt metoda puterii si metoda puterii inverse. Aceste metode suntimplementate ın mod profesional ıntr-o faimoasa procedura iterativa numita algoritmul QR.Algoritmul QR construieste iterativ o secventa de matrice ortogonal asemenea, convergentacatre o forma Schur a matricei originale. Proprietatile excelente de convergenta se datorezafolosirii unor deplasari optimale si unei serii de subtilitati tehnice cum ar fi utilizarea ex-clusiva a formei Hessenberg si monitorizarea structurala. Toate variantele algoritmului QRsunt organizate ın doua etape:

1. Reducerea la forma superior Hessenberg printr-o transformare ortogonala de asemanareH = UT AU calculata ıntr-un numar finit de flopi (O( 5

3n3)) (algoritmul HQ).

2. Reducerea iterativa la forma Schur, i.e. construirea iterativa a unui sir de matriceH1 = H, H2, . . . ,Hk, . . ., convergent catre o forma Schur a lui A. Pentru o matrice cuun spectru real, fiecare iteratie QR este definita prin:

1. [ Qk, Rk ] = qr(Hk − µkIn) % Factorizarea QR a lui Hk − µkIn

2. Hk+1 ← Rk ∗Qk + µkIn,

unde µk = Hk(n, n) este deplasarea Rayleigh. Daca matricea A are valori propriicomplexe, atunci se foloseste o varianta mai complicata cu pasi dubli si deplasariimplicite (vezi cursul).

4.2. PROBLEME REZOLVATE 3

4.2 Probleme rezolvate

4.2.1 Aspecte teoretice

Problema 4.1 Fie A, B ∈ IRn×n doua matrice asemenea, i.e. exista o matrice nesingularaT ∈ IRn×n astfel ıncat B = T−1AT . Demonstrati ca

a. det A = det B, b. trA = trB,

unde trA =∑n

i=1 aii este, prin definitie, urma lui A. Ca aplicatie, demonstrati ca dacaλ(A) = {λi, λi, . . . , λi} este spectrul valorilor proprii ale lui A, atunci

c. detA =n∏

i=1

λi, d. trA =n∑

i=1

λi.

Rezolvare. Este bine-cunoscut ca det(GH) = det G det H pentru toate matricelepatrate G si H. Evident, acest lucru implica faptul ca pentru toate matricele nesingulare G

avem det G−1 =1

det G. a. Din B = T−1AT rezulta ca det B =

1detT

detA det T = det A.

b. Sa notam X = T−1. Atunci XT = TX = In, i.e. X(i, :)T (:, j) = T (i, :)X(:, j) = δij ,unde scalarul δij = 0 pentru toti i 6= j si δij = 1 daca i = j. Prin urmare

trB =n∑

k=1

B(k, k) =n∑

k=1

X(k, :)AT (:, k) =n∑

k=1

n∑

i=1

n∑

j=1

X(k, i)A(i, j)T (j, k) =

=n∑

i=1

n∑

j=1

A(i, j)n∑

k=1

T (j, k)X(k, i) =n∑

i=1

n∑

i=1

A(i, j)T (j, :)X(:, i) =

=n∑

i=1

n∑

j=1

A(i, j)δij =n∑

i=1

A(i, i) = trA.

c. Fie S = QT AQ forma Schur reala pentru A. Atunci det A = det S =∏p

i=1 detSii.Daca Sii este un bloc scalar det Sii = Sii; daca Sii este un bloc 2 × 2, atunci det Sii =Sii(1, 1)Sii(2, 2) − Sii(1, 2)Sii(2, 1) = λ1(Sii)λ2(Sii). Astfel, detA =

∏pi=1 Πni

j=1λj(Sii) =∏ni=1 λi(A).

d. Fie din nou S = QT AQ forma Schur reala pentru A. Atunci trA = trS =∑p

i=1 trSii.Daca Sii este un bloc scalar trSii = Sii; daca Sii este un bloc 2×2, atunci trSii = Sii(1, 1)+Sii(2, 2) = λ1(Sii) + λ2(Sii). Astfel, trA =

∑pi=1

∑ni

j=1 λj(Sii) =∑n

i=1 λi(A).

Problema 4.2 Fie A ∈ IRn×n si λ(A) = {λ1, λ2, . . . , λn} spectrul valorilor sale proprii. a.Demonstrati ca valorile proprii sunt situate ın urmatorul domeniu al planului complex

D =n⋃

i=1

Di,


unde Di sunt discurile

Di = {z ∈ IC | |z − aii| ≤n∑

j=1j 6=i

|aij |}, i = 1 : n,

numite discurile lui Gershgorin.

b. Numarul ρ(A) def= maxi=1:n(|λi|) se numeste raza spectrala a lui A. Demonstratiinegalitatea

ρ(A) ≤ ‖A‖,unde ‖ · ‖ este o familie oarecare de norme matriceale consistente.

c. Demonstrati urmatoarele inegalitati

c1. |detA| ≤n∏

i=1

n∑

j=1

|aij | , c2. |detA| ≤

n∏

j=1

(n∑

i=1

|aij |)

.

Rezolvare. a. Fie x un vector propriu asociat valorii proprii λ ∈ λ(A). Atunci linia ia egalitatii vectoriale Ax = λx poate fi scrisa sub forma

(λ− aii)xi =n∑

j=1j 6=i

aijxj .

Prin urmare, este evident ca |λ− aii||xi| ≤∑n

j=1j 6=i

|aij ||xj |. Alegand acum linia i astfel ıncat

|xi| = maxk=1:n(|xk|) 6= 0, i.e.|xj ||xi| ≤ 1 pentru toti j, avem

|λ− aii| ≤n∑

j=1j 6=i

|aij | |xj ||xi| ≤

n∑j=1j 6=i

|aij |,

i.e. λ ∈ Di.b. Datorita consistentei familiei de norme matriceale, pentru toate λ ∈ λ(A) si toti

vectorii proprii asociati x cu ‖x‖ = 1, avem |λ| = ‖λx‖ = ‖Ax‖ ≤ ‖A‖‖x‖ = ‖A‖. Prinurmare inegalitatea este adevarata.

c. Sa demonstram mai ıntai inegalitatea c1. Daca A are (cel putin) o linie zero, atunciinegalitatea devine egalitate. Daca nu este cazul, sa definim δi =

∑nj=1 |aij | > 0, i = 1 : n,

si D = diag(δ1, δ2, . . . , δn). Matricea B = D−1A are ρ(B) ≤ ‖B‖∞ ≤ 1. Prin urmare,|detB| =

∏ni=1 |λi(B)| ≤ 1. Inegalitatea dorita se obtine din |detA| = |detD| · |detB| ≤

|detD|. Pentru a demonstra c2 aplicam inegalitatea din cazul c1 matricei AT .

Problema 4.3 Fie A ∈ IRm×n, B ∈ IRn×m cu m > n. Demonstrati ca λ(AB) = λ(BA) ∪{0, 0, . . . , 0}, i.e. ca, ın particular, matricele AB si BA au aceleasi valori proprii nenule.

Rezolvare. Matricele (m + n) × (m + n) C =[

AB 0B 0

]si D =

[0 0B BA

]sunt

asemenea deoarece D = TCT−1 cu T =[

Im A0 In

]avand ca inversa T−1 =

[Im −A0 In

]


(verificati!). Prin urmare λ(C) = λ(D). Dar C si D sunt matrice bloc inferior triunghiulare.Astfel, λ(C) = λ(AB) ∪ {0, 0, . . . , 0} si λ(D) = {0, 0, . . . , 0} ∪ λ(BA). Atunci, deoarecem > n, λ(AB) = λ(BA) ∪ {0, 0, . . . , 0}. Evident, daca m = n atunci λ(AB) = λ(BA).

Problema 4.4 a. Fie λ1, λ2 doua valori proprii diferite ale matricei A ∈ IRn×n si x1

un vector propriu al lui A asociat lui λ1 si y2 un vector propriu al lui AT asociat lui λ2.Demonstrati ca cei doi vectori sunt ortogonali i.e. yT

2 x1 = 0.b. Daca λ ∈ λ(A) este o valoare proprie distincta a lui A si x, respectiv y, sunt vectori

proprii la dreapta, respectiv la stanga, asociati lui λ, atunci yT x 6= 0. Dati un exemplu ıncare aceasta proprietate nu este satisfacuta daca λ nu este o valoare proprie distincta.

Rezolvare. a. Avem Ax1 = λ1x1 si AT y2 = λ2y2 sau yT2 A = λ2y

T2 . Prin urmare

λ1yT2 x1 = λ2y

T2 x1 sau (λ1− λ2)yT

2 x1 = 0 care implica yT2 x1 = 0 din cauza ca (λ1− λ2) 6= 0.

b. Fara sa pierdem caracterul general, putem presupune ca ‖x‖2 = 1. Conform culema de deflatie ortogonala, daca matricea [ x X ] este ortogonala, atunci B = XT AX =[

λ bT

0 C

]. Acum, daca vectorul y este un vector propriu la stanga al lui A (i.e. yT A = λyT ,

sau yT A = λyT ; astfel, un vector propriu la stanga al lui A este un vector propriu al luiAT ), atunci yT XB = yT XXT AX = yT AX = λyT X, deci zT B = λzT unde z = XT y. Prinurmare z este un vector propriu la stanga al lui B. Dar λ este o valoare proprie distincta,deci matricea λIn−1 − C este nesingulara. Rezulta ca z(2 : n) = (λIn−1 − CT )−1bz1, cuz1 = xT y. Daca z1 = 0, atunci z(2 : n) = 0 si, prin urmare, z = 0. Dar un vector propriu nupoate sa fie zero. Astfel z1 = xT y 6= 0. Drept exemplu pentru situatia cand valoarea proprie

nu este distincta si conditia nu este satisfacuta, fie matricea A =[

0 10 0

]ai carei vectori

proprii (la dreapta) au forma x =[

α0

]si cei la stanga au forma y =

[0β

], α, β ∈ IR,

α 6= 0, β 6= 0; astfel yT x = 0.

4.2.2 Aspecte calculatorii

Problema 4.5 Fie A ∈ IRm×n o matrice data.a. Presupunand ca vectorul propriu x ∈ IRn al lui A este cunoscut, scrieti un algoritm

pentru calculul valorii proprii asociate λ.b. Presupunand ca valoarea proprie λ ∈ λ(A) este cunoscuta, scrieti un algoritm pentru

calculul unui vector propriu asociat.

Rezolvare. a. Avem λx = Ax, deci λxT x = xT Ax si, din cauza ca x 6= 0, xT x =‖x‖2 6= 0. Prin urmare

λ =xT Ax

xT x=

∑ni=1

∑nj=1 aijxixj∑ni=1 x2

i

si algoritmul este

1. σ = 02. τ = 0


3. pentru i = 1 : n1. yi = 02. pentru j = 1 : n

1. yi = yi + aijxj

3. σ = σ + xiyi

4. τ = τ + x2i

2. λ = σ/τ

b. Daca o valoare proprie λ ∈ λ(A) este cunoscuta, atunci un vector propriu asociatpoate fi calculat prin calculul unei solutii nenule a sistemului liniar omogen (λI − A)x =0. Pentru aceasta vom aplica eliminarea gaussiana cu pivotare partiala (algoritmul GPP)matricei λI −A

Mn−1Pn−1 · · ·M2P2M1P1(λI −A) = U

unde U este superior triunghiulara. Evident, solutiile sistemului Ux = 0 sunt solutii alesistemului (λI−A)x = 0. Matricea superior triunghiulara U este singulara deoarece matriceaλI−A este singulara. Astfel, matricea U are cel putin un element diagonal nul. Fie ukk = 0primul element diagonal nul a lui U . Sistemul Ux = 0 poate fi scris sub forma

U11 U12 U13

0 0 U23

0 0 U33

x′

xk

x′′′

=

000

,

unde U11 = U(1 : k− 1, 1 : k− 1) este un bloc superior triunghiular nesingular, U12 = U(1 :k − 1, k), U23 = U(k, k + 1 : n), U33 = U(k + 1 : n, k + 1 : n). O solutie nenula poate fiobtinuta alegand

x′′′ = 0, xk = α 6= 0, x′ = −U−111 U12α

unde α este un scalar arbitrar nenul.

Prin urmare, un algoritm pentru calculul unui vector propriu, atunci cand valoareaproprie asociata este cunoscuta, este

1. [M, p, U ] =GPP(λI −A)2. k = 13. cat timp ukk 6= 0

1. k ← k + 14. pentru i = k + 1 : n

1. xi = 05. xk = 16. daca k > 1

1. x(1 : k − 1) =UTRIS(U(1 : k − 1, 1 : k − 1),−U(1 : k − 1, k))

In algoritmul de mai sus scalarul α a fost considerat egal cu 1.

Problema 4.6 a. Calculati forma Schur reala a matricei A =[

1 −22 −3

].

b. Calculati forma Schur reala a matricei A =[

1 1−1 1

].


c. Presupunand ca o matrice data A are un spectru real (i.e. λ(A) ⊂ IR) si ca dispunemde o procedura cu sintaxa

x = eigenvector(A)

pentru calculul unui vector propriu al matricei A (fara a se cunoasste valoarea proprieasociata), scrieti un algoritm pentru calculul unei forme Schur a lui A.

Rezolvare. a. Polinomul caracteristic al lui A este p(λ) = det(λI2 −A) = (λ− 1)(λ +3)+4 = (λ− 1)2, deci valorile proprii ale lui A sunt ambele egale cu −1. Sistemul Ax = −x

conduce la x1 = x2, astfel ca toti vectorii proprii ai lui A sunt de forma x = α

[11

],

α 6= 0. Pentru a aplica lema de deflatie ortogonala, fie vectorul propriu de norma euclidiana

unitate x = (1/√

2)[

11

]. Considerand matricea ortogonala Q = [x y], unde e.g. y =

(1/√

2)[

1−1

](avand, evident, ‖y‖ = 1), rezulta S = QT AQ =

[ −1 40 −1

], care este o

forma Schur a lui A. Observati ca forma Schur poate fi obtinuta pe cale indicata chiar dacaA este defectiva (i.e. nu este simpla).

b. Polinomul caracteristic al lui A este p(λ) = det(λI2−A) = (λ− 1)2 +1, deci valorileproprii ale lui A nu sunt reale. Astfel o forma Schur reala (FSR) a lui A este chiar A.

c. Vom aplica ın mod sistematic lema de deflatie ortogonala. Dupa cum se demonstreazala curs, daca x este un vector propriu de norma euclidiana unitate al lui A si Q = [ x Y ]este o matrice ortogonala atunci matricea

B = QT AQ =[

λ bT

0 C

]

este ın forma Schur ın prima coloana. Pentru a continua, putem aplica aceeasi procedurablocului C, etc. Algoritmul dorit se poate baza pe urmatoarea schema de calcul:

1. pentru k = 1 : n− 11. xk = eigenvector(A(k : n, k : n))

2. xk =xk

‖xk‖3. Se calculeaza o matrice ortogonala Qk astfel ıncat Qke1 = xk

4. A(k : n, k : n) = QTk A(k : n, k : n)

5. A( : , k : n) = A( : , k : n)Qk.

In final, matricea A va fi suprascrisa cu forma sa Schur.Pentru a detalia schema de mai sus, remarcati faptul ca matricea Qk din instructiunea

1.3 poate fi un reflector U = I − uuT /β astfel ıncat (Uxk)(2 : n − k + 1) = 0. Folosindprocedurile cunoscute pentru calculul vectorului u si al scalarului β ce defineste reflectorul sipentru ınmultirea la stanga si la dreapta a unei matrice cu un reflector si folosind un vectorgeneric x pentru xk, algoritmul detaliat devine:

1. pentru k = 1 : n− 11. x = eigenvector(A(k : n, k : n))

2. α =√∑n−k+1

i=1 x2i


3. pentru i = 1 : n− k + 11. xi ← xi/α

4. u1 = x1 − 1, % σ = −15. pentru i = 2 : n− k + 1

1. ui = xi

6. β = −u1, % σ = −17. pentru j = k : n

1. τ = (∑n

i=k ui−k+1aij) /β2. pentru i = k : n

1. aij ← aij − τui−k+1

8. pentru i = 1 : n

1. τ =(∑n

j=k aijuj−k+1

)/β

2. pentru j = k : n1. aij ← aij − τuj−k+1.

Problema 4.7 a. Fie A =[

α γγ β

]∈ R2×2 o matrice simetrica. Dati un algoritmu

pentru calculul valorilor si vectorilor proprii ai lui A. Retineti faptul ca valorile proprii suntreale si vectorii proprii sunt ortogonali.

b. Aratati ca o matrice simetrica reala are spectrul real si vectori proprii ortogonalidoi cate doi. Ce simplificari ale calculului apar ın algoritmul QR simetric fata de cazulnesimetric?

Rezolvare. a. Polinomul caracteristic al lui A este det(λI2−A) = (λ−α)(λ−β)−γ2 =λ2 − (α + β)λ + αβ − γ2 si discriminantul sau este δ = (α + β)2 − 4(αβ − γ2) = (α− β)2 +

4γ2 ≥ 0. Prin urmare valorile proprii sunt ambele reale si egale cu λ1 =α + β +

√δ

2si

λ2 =α + β −

√δ

2. Vectorii proprii pot fi calculati prin rezolvarea sistemelor liniare singulare

omogene Ax1 = λ1x1 and Ax2 = λ1x2. Sistemul Ax1 = λ1x1 poate fi scris{

(λ1 − α)x11 − γx21 = 0−γx11 + (λ1 − β)x21 = 0

Daca γ = 0, atunci λ1 = α si o solutie nenula a sistemului de mai sus este x1 =[

µ0

]cu

µ 6= 0 arbitrar. Daca γ 6= 0, atunci considerand x11 = µ 6= 0, o solutie nenula a sistemului

de mai sus este x1 =

µ(λ1 − α)µ

γ

. Sistemul Ax2 = λ1x2 poate fi scris similar

{(λ2 − α)x21 − γx22 = 0−γx21 + (λ2 − β)x22 = 0

Daca γ = 0, atunci λ2 = β si o solutie nenula a sistemului de mai sus este x2 =[

0ν

]cu

ν 6= 0 arbitrar. Daca γ 6= 0, atunci considerand x22 = ν 6= 0, o solutie nenula a sistemului


de mai sus este x2 =

(λ2 − β)νγν

. Daca γ = 0, cei doi vectori proprii sunt evident

ortogonali. Daca γ 6= 0, atunci

xT2 x1 =

[(λ2 − β)ν

γν

]

µ(λ1 − α)µ

γ

=

(λ2 − β)νγ

µ + ν(λ1 − α)µ

γ=

=µν

γ(λ1 − α + λ2 − β) = 0.

Deci, cei doi vectori proprii sunt ortogonali.b. Daca A ∈ IRn×n este simetrica, i.e. AT = A, fie S = QT AQ forma sa Schur reala,

care este o structura cvasi-triunghiulara. Avem ST = QT AT Q = QT AQ = S, deci matriceaS este simetrica si, prin urmare, are o structura cvasi-diagonala S = diag(S11, S22, . . . , Spp)cu blocurile diagonale 2×2 avand valori proprii complex conjugate. Dar blocurile diagonalesunt de asemenea simetrice si am vazut mai sus ca o matrice 2 × 2 simetrica are valoriproprii reale. Astfel ca toate blocurile diagonale ale unei forme Schur reale a unei matricesimetrice sunt 1× 1, i.e. forma Schur reala a unei matrice simetrice este diagonala (ın acestcaz forma canonica Jordan si forma Schur sunt identice), toate valorile proprii sunt realesi din S = QT AQ = Λ = diag(λ1, λ2, . . . , λn) avem AQ = QΛ, i.e. AQ(:, j) = QΛ(:, j) =λjQej = λjQ(:, j), j = 1 : n. Prin urmare, xj = Q(:, j) sunt vectori proprii ai matriceisimetrice A, asociati valorilor proprii λj si, fiind coloane ale unei matrice ortogonale, sunttoti ortogonali doi cate doi. Aceste proprietati utile aduc simplificari importante algoritmuluiQR simetric. Astfel matricea superior Hessenberg obtinuta ın prima etapa a algoritmuluiQR este simetrica si prin urmare tridiagonala. De asemenea, toate matricele din secventaQR sunt tridiagonale si convergenta la forma diagonala este foarte rapida. Pentru mai multedetalii consultati cursul si bibliografia.

Problema 4.8 Elaborati un algoritm pentru reducerea unei matrice A ∈ Rn×n la forma su-perior Hessenberg H = TAT−1, unde T este un produs de transformari gaussiene elementarestabilizate MkPk (de tipul celor ıntalnite ın algoritmul GPP).

Rezolvare. Ideea este aceeasi ca si ın cazul algoritmului de reducere ortogonala HQ sise bazeaza pe faptul ca matricele de permutare Pk = Pkik

, ik ≥ k si matricele elementare

inferior triunghiulare Mk = In − mkeTk au structura

[Ik−1 0

0 X

]la fel cu a matricelor

MkPk, P−1k = PT

k , M−1k = In + mkeT

k si P−1k M−1

k . Prin urmare, forma Hessenberg poate fiobtinuta prin urmatoarele transformari de asemanare

A ← H = TAT−1 = Mn−1Pn−1 · · ·M3P3M2P2A(M2P2)−1(M3P3)−1 · · · (Mn−1Pn−1)−1 =

= Mn−1Pn−1 · · ·M3P3M2P2AP2M−12 P3M

−13 · · ·Pn−1M

−1n−1.

Luand ın calcul faptul ca o matrice de permutare Pkikeste definita complet de cei doi ıntregi

(k, ik) putem scrie urmatoarea schema de calcul

1. pentru k = 1 : n− 2


1. Se determina ik+1 astfel ıncat |aik+1k| = maxi=k+1:n(|aik|)2. A(ik+1, k : n) ↔ A(k + 1, k : n)3. Se determina matricea elementara inferior triunghiulara

Mk+1 astfel ıncat (Mk+1A)(k + 2 : n, k) = 04. A = Mk+1A5. A( : , k + 1) ↔ A( : , ik+1)6. A = AM−1

k+1.

Folosind cunostintele de utilizare a transformarilor gaussiene, schema de calcul de maisus este detaliata ın urmatorul algoritm.

1. pentru k = 1 : n− 21. ik+1 = k + 12. ν = |ak+1,k|3. pentru i = k + 2 : n

1. daca |aik| > ν1. ik+1 = i2. ν = |aik|

4. pentru j = k : n1. aik+1,j ↔ ak+1,j

5. pentru i = k + 2 : n1. µi,k+1 = aik/ak+1,k

2. aik = 06. pentru j = k + 1 : n

1. pentru i = k + 2 : n1. aij ← aij − µi,k+1ak+1,j

7. pentru i = 1 : n1. aik+1,k+1 ↔ ai,k+1

8. pentru j = k + 1 : n1. pentru i = k + 2 : n

1. aij ← aij − µi,k+1ak+1,j

Solutia este de doua ori mai eficienta decat algoritmul HQ, dar acesta din urma este preferatdatorita avantajelor transformarilor ortogonale.

Problema 4.9 Fie H ∈ Rn×n o matrice superior Hessenberg. Descrieti detaliile unui pasQR simplu cu deplasare explicita folosind valoarea optima µ = hnn a deplasarii.

Rezolvare. Fiind data matricea H curenta din sirul QR, calculul H ← H ′ al succesoareiei reprezinta un pas simplu QR. Versiunea cu deplasare explicita este definitaa prin

{[ Q,R ] = qr(H − µIn) %factorizareaQRH ← H ′ = RQ + µIn

Matricea H − µI este superior Hessenberg; factorizarea sa QR poate fi calculata folosindo secventa de rotatii Givens pentru triangularizarea ortogonala a matricei H − µI, i.e.Pn−1,n · · ·P23P12(H − µI) = R; atunci Q = PT

12PT23 · · ·PT

n−1,n. Inmultirea matricelor RQ sepoate face fara calculul explicit al lui Q, folosind expresia RQ = RPT

12PT23 . . . PT

n−1,n. Toatecalculele pot fi facute prin suprascrierea matricei H. Schema de calcul pentru pasul QRsimplu cu deplasare explicita este:


1. µ ← hnn

2. H ← H − µI3. pentru k = 1 : n− 1

1. Se calculeaza rotatia Givens Pk,k+1 astfel ıncat (Pk,k+1H)(k + 1, k) = 02. H ← Pk,k+1H

3. pentru k = 1 : n− 12. H ← HPT

k,k+1

5. H ← H + µI

Notand cu (ck, sk) scalarii definitorii pentru rotatia Pk,k+1, schema de mai sus poate fidetaliata ın urmatorul algoritm.

1. µ ← hnn

2. pentru i = 1 : n1. hii ← hii − µ

3. pentru k = 1 : n− 1

1. ρ =√

h2kk + h2

k+1,k

2. ck = hk,k/ρ3. sk = hk+1,k/ρ4. hk,k = ρ5. hk+1,k = 06. pentru j = k + 1 : n

1. τ = hkj

2. hkj ← ckτ + skhk+1,j

3. hk+1,j ← −skτ + ckhk+1,j

4. pentru k = 1 : n− 11. pentru i = 1 : k + 1

1. τ = hik

2. hik ← τck + hi,k+1sk

3. hi,k+1 ← −τsk + hi,k+1ck

5. pentru i = 1 : n1. hii ← hii + µ

Problema 4.10 Elaborati un algoritm pentru calculul valorilor proprii ale matricei A =In + uvT , unde u, v ∈ Rn.

Rezolvare. Daca Q = QT este reflectorul Householder pentru care Qu = ‖u‖e1, atuncimatricea

A′ = QT AQ = I + ‖u‖e1vT QT = I + e1w

T ,

unde w = ‖u‖Qv, este ortogonal asemenea cu A si superior triunghiulara. Prin urmarevalorile proprii ale lui A sunt elementele diagonale ale lui A′. Avem a′11 = 1 + w1 =1 + ‖u‖vT q1, unde q1 = Qe1 este prima coloana a lui Q si a′ii = 1 pentru i ≥ 2. Darq1 = u/‖u‖, prin urmare a′11 = 1 + vT u. Astfel, λ1 = 1 + vT u si λi = 1, i = 2 : n este ovaloare proprie cu ordinul de multiplicitate n− 1.


Problema 4.11 a. Fiind dat un vector x ∈ IRn nenul, scrieti un algoritm pentru calcululunui vector v ∈ IRn astfel ıncat vT x = 1.

b. Presupunem ca matricea A ∈ ICn×n are valorile proprii reale λi, i = 1 : n, si caxi, i = 1 : n, sunt vectori proprii asociati. Fie v ∈ IRn un vector astfel ıncat vT x1 = 1 siconsideram matricea B = (In − x1v

T )A. Aratati ca λ(B) = { 0, λ2, . . . , λn }, si ca vectoriixB

1 = x1, xBi = xi − (vT

1 xi)x1 sunt vectori proprii ai matricei B.

Rezolvare. a. Fie U1 = In − uuT /β un reflector Householder astfel ıncat U1x = ρe1,ρ 6= 0. Atunci

v =1ρU1e1 =

1ρ(e1 − uτ)

unde τ = u1/β, este vectorul pe care ıl cautam. Intr-adevar, vT x =1ρeT1 U1x = 1. Vectorul

v poate fi calculat dupa cum urmeaza:

1. ρ = −√∑n

i=1 x2i

2. u1 = ρ + x1

3. pentru i = 2 : n1. ui = xi

4. β = u1ρ5. τ = u1β6. v1 = (1− u1τ)/ρ6. pentru i = 2 : n

1. vi = −uiτ/ρ

b. Avem Bx1 = (In−x1vT )Ax1 = λ1(In−x1v

T )x1 = 0; deci x1 este un vector propriu al luiB asociat unei valori proprii nule a lui B. De asemenea BxB

i = (In−x1vT )A(xi−(vT

1 xi)x1) =(In − x1v

T )(λixi − λ1x1(vT1 xi)) = λix

Bi , i = 2 : n.

Problema 4.12 Scrieti un algoritm pentru calculul polinomului caracteristic al unei ma-trice patrate A ∈ IRn×n date.

Rezolvare. Cea mai buna metoda numerica actuala consta ın calculul initial al valorilorproprii λi, folosind algoritmul QR, iar apoi ın calculul p(λ) =

∏ni=1(λ − λi). Pentru a

determina coeficientii lui p(λ), folosim schema de calcul

1. p(λ) = 12. for i = 1 : n

1. p(λ) = p(λ)(λ− λi)

Daca, la ınceputul pasului i avem

p(λ) = p(i−1)(λ) = λi−1 + p(i−1)2 λi−2 + . . . + p

(i−1)i−1 λ + p

(i−1)i

procesarea la pasul i duce la

p(λ) = p(i)(λ) = λi + (p(i−1)2 − 1 · λi)λi−1 + . . . + (p(i)

i−1 − p(i−1)i λi)λ + (0− p

(i−1)i λi) =

4.3. PROBLEME PROPUSE 13

= λi−1 + p(i)2 λi−2 + . . . + p

(i)i λ + p

(i)i+1.

Prin urmare, p(i)1 = p

(i−1)1 = 1 si p

(i)j = p

(i−1)j −p

(i−1)j+1 λi pentru j ≥ 2. Calculand coeficientii

ın ordinea inversa putem suprascrie coeficientii de la pasul anterior precum ın urmatorulalgoritm. Astfel, cel mai bun algoritm pentru calculul coeficientilor polinomului caracteristiceste:

1. Folosind algoritmul QR, calculati λi, i = 1 : n2. p1 = 13. pentru i = 2 : n + 1

1. pi = 04. pentru i = 1 : n

1. pentru j = i + 1 : −1 : 21. pj ← pj − λipj−1.

Evident, efortul principal de calcul consta ın determinarea valorilor proprii.

4.3 Probleme propuse

Problema 4.13 Fie date matricele

A =

3 −3 2−1 5 −2−1 3 0

, B =

−4 0 8−8 3 9−4 −1 9

.

Folosind definitiile valorilor si vectorilor proprii calculati valorile si vectorii proprii pentruA si B. Dupa aceea calculati-le si formele Schur.

Problema 4.14 a. O matrice patrata se numeste nilpotenta daca exista un numar ıntregk ≥ 0 astfel ıncat Ak = 0. Demonstrati ca toate valorile proprii ale unei matrice nilpotentesunt nule. Dati un exemplu de matrice 2× 2 nilpotenta.

b. O matrice patrata se numeste idempotenta daca A2 = A. Demonstrati ca toatevalorile proprii ale unei matrice idempotente sunt 0 sau 1. Dati un exemplu de matrice 2×2idempotenta care sa nu fie I2.

Problema 4.15 a. Calculati valorile proprii si un set de vectori proprii pentru un reflectorHouseholder Uk. b. Calculati valorile proprii si un set de vectori proprii pentru o rotatieplana Pij .

Problema 4.16 Care va fi rezultatul aplicarii metodei puterii matricelor:

A =

5 8 10 1 20 0 2

, B =

[1 −11 3

],

cu aproximarile initiale de rigoare? Ce ınseamna ”de rigoare” ın aceste cazuri?


Problema 4.17 Presupunem ca avem date matricea A ∈ IRn×n si vectorul z ∈ IRn. Scrietiun algoritm pentru calculul unei matrice ortogonale Q astfel ıncat QT AQ sa fie superiorHessenberg si QT z coliniara cu vectorul e1.

Problema 4.18 Fie data o pereche (valoare proprie, vector propriu)= (λ, x) a unei matricesuperior Hessenberg H ∈ IRn×n. Scrieti un algoritm pentru calculul unei matrice ortogonale

Q astfel ıncat matricea QT HQ sa aiba structura QT HQ =[

λ fT

0 G

], unde matricea

G ∈ IR(n−1)×(n−1) este superior Hessenberg.

Problema 4.19 Fie H ∈ IRn×n o matrice superior Hessenberg data si urmatoarea proce-dura pentru calculul matricei succesoare H ← H ′ ıntr-un algoritm iterativ:

1. Se aplica algoritmul GPP Mn−1Pn−1 . . . M1P1H = R pentru a obtine matricea supe-rior triunghiulara R.

2. H ← H ′ = RP1M−11 . . . Pn−11M−1

n−1,

Demonstrati ca:a. matricea succesoare H ′ este superior Hessenberg si ca,b. matricea succesoare H ′ este asemenea cu H.

Problema 4.20 a. Consideram matricea H ∈ IR2×2 cu spectru real si fie Hknot=

[α γε β

]

matricea curenta din sirul QR pentru matricea H. Folosind deplasarea µk = β calculatimatricea succesoare Hk+1. Examinati elementul Hk+1(2, 1) si deduceti unele proprietati deconvergenta ale sirului QR.

b. Consideram matricea simetrica T ∈ IR2×2 si fie Tknot=

[α εε β

]matricea curenta

a sirului QR pentru matricea T . Folosind deplasarea µk = β calculati matricea succe-soare Tk+1. Examinati elementele extradiagonale Tk+1(1, 2) = Tk+1(2, 1) si deduceti uneleproprietati de convergenta ale sirului QR simetric .

4.4 Bibliografie

1. B. Jora, B. Dumitrescu, C. Oara, Numerical Methods, UPB, Bucharest, 1995.

2. G.W. Stewart, Introduction to Matrix Computations, Academic Press, 1973.

3. G. Golub, Ch. Van Loan, Matrix Computations, 3-rd edition, John Hopkins Uni-versity Press, 1998.

4. G.W. Stewart, Matrix Algorithms, vol.1: Basic Decompositions, SIAM, 1999.

5. G.W. Stewart, Matrix Algorithms, vol.2: Eigensystems, SIAM, 2001.

6. B. Dumitrescu, C. Popeea, B. Jora, Metode de calcul numeric matriceal. Al-goritmi fundamentali, ALL, Bucuresti, 1998.

4.5. PROGRAME MATLAB 15

4.5 Programe MATLAB

Problema 4.5 a1: function [ lambda]= p45a(A,x)2: %——————————–3: % Rezolvarea problemei 4.5 subpunctul a).4: % Algoritmul calculeaza catul Rayleigh definit de perechea (A,x).5: % Daca x este un vector propriu al matricei patratice A, atunci6: % lambda este valoarea proprie asociata vectorului propriu x.7: % Apel:8: % lambda=p45a(A,x);9: % Stamatescu Grigore, mai 2006.10: %——————————–11: [m,n ]=size(A);12: % verificam daca matricea A este patratica13: if m∼=n,14: error(’Matricea A nu este patratica’);15: end16: sigma=0;17: tau=0;18: for i=1:n19: y(i)=0;20: for j=1:n21: y(i)=y(i)+A(i,j)*x(j);22: end23: sigma=sigma+x(i)*y(i);24: tau=tau+x(i)*x(i);25: end26: lambda=sigma/tau;

1: function [x ]= eigenvector(A)2: %——————————–3: % Rezolvarea problemei 4.5 subpunctul a).4: % Algoritmul calculeaza un vector propriu al matricei patratice A, folosind5: % functia MATLAB eig, necesar pentru testarea programului precedent.6: % Stamatescu Grigore, mai 2006.7: %——————————–8: [m,n ]=size(A);9: % verificam daca A este patratica10: if m∼=n,11: error(’Matricea A nu este patratica’);12: end13: [V,D]=eig(A);14: x=V(1:n);


Problema 4.5 b1: function [x ]= p45b(A,lambda)2: %——————————-3: % Algoritmul calculeaza vectorul propriu x al unei matrice A atunci4: % cand valoarea proprie lambda corespunzatoare este cunoscuta.5: % Apelul:6: % x=p45b(A,lambda)7: % Stamatescu Grigore, mai 2006.8: %——————————-9: [m,n ]=size(A);10: % verificam daca matricea A este patratica11: if m∼=n,12: error(’Matricea A nu este patratica’);13: end14: T=lambda*eye(size(A))-A;15: [T,M,p ]=GPP(T);16: k=1;17: for i=1:n-118: if A(i,i)∼=0,19: k=k+1;20: end21: end22:23: for i=k+1:n24: x(i)=0;25: end26: x(k)=1;27:28: if k>1,29: for i=1:k-130: x(i)=UTRIS(T(i,i),-T(i,k));31: end32: end

Problema 4.61: function [A]= Schr(A)2: %——————————-3: % Algoritmul calculeaza o forma Schur a matricei patratice A cu spectru real4: % folosind o procedura de calcul a unui vector propriu.5: % Apel:6: % A=Schr(A);7: % Stamatescu Grigore, mai 2006.8: %——————————-9: [m,n ]=size(A);10: if m∼=n,


11: error(’Matricea A nu este patratica’);12: end13: [V,D]=eig(A); % ineficient!14: for t=1:n15: if isreal(D(t,t))∼=1,16: error(’Matricea A nu are spectrul valorilor proprii real!’);17: end18: end19:20: for k=1:n-121: x=eigenvector(A(k:n,k:n));22: suma=0;23: for i=1:n-k+124: suma=suma+x(i)*x(i);25: end26: alfa=sqrt(suma);27: for i=1:n-k+128: x(i)=x(i)/alfa;29: end30: u(1)=x(1)-1;31: for i=2:n-k+132: u(i)=x(i);33: end34: beta=-u(1);35: for j=k:n36: suma1=0;37: for i=k:n38: suma1=suma1+u(i-k+1)*A(i,j);39: end40: tau=suma1/beta;41: for i=k:n42: A(i,j)=A(i,j)-tau*u(i-k+1);43: end44: end45: for i=1:n46: suma2=0;47: for j=k:n48: suma2=suma2+A(i,j)*u(j-k+1);49: end50: tau=suma2/beta;51: for j=k:n52: A(i,j)=A(i,j)-tau*u(j-k+1);53: end54: end55: end


Problema 4.81: function [A]= Hessenberg(A)2: %———————————3: % Algoritmul reduce matricea A patratica la forma Hessenberg.4: % Apel:5: % A=Hessenberg(A)6: % Stamatescu Grigore mai 20067: %———————————8: [m,n ]=size(A);9: % verificam daca A este patratica10: if m∼=n,11: error(’Matricea A nu este patratica’);12: end13:14: for k=1:n-215: z(k+1)=k+1;16: niu=abs(A(k+1,k));17: for i=k+2:n18: if abs(A(i,k))>niu,19: z(k+1)=i;20: niu=abs(A(i,k));21: end22: end23: for j=k:n24: A(z(k+1),j)=A(k+1,j);25: end26: for i=k+2:n27: miu(i,k+1)=A(i,k)/A(k+1,k);28: A(i,k)=0;29: end30: for j=k+1:n31: for i=k+2:n32: A(i,j)=A(i,j)-miu(i,k+1)*A(k+1,j);33: end34: end35: for i=1:n36: A(z(k+1),k+1)=A(i,k+1);37: end38: for j=k+1:n39: for i=k+2:n40: A(i,j)=A(i,j)-miu(i,k+1)*A(k+1,j);41: end42: end43: end

Problema 4.9


1: function [H]= pasQR(H)2: %———————————3: % Fiind data o matrice H superior Hessenberg algoritmul implementeaza4: % un pas QR cu deplasare explicita folosind deplasarea miu=H(n,n)5: % Apel:6: % H=pasQR(H)7: % Stamatescu Grigore, mai 2006.8: %———————————9: [m,n ]=size(H);10: % verificam daca H este patratica11: if m∼=n,12: error(’Matricea H nu este patratica’);13: end14: % analizam daca matricea este superior Hessenberg15: for i=3:n16: for j=1:i-217: if H(i,j)∼=018: error(’Matricea H nu e superior Hessenberg’)19: end20: end21: end22: miu=H(n,n);23: for i=1:n24: H(i,i)=H(i,i)-miu;25: end26: for k=1:n-127: ro=sqrt(H(k,k)*H(k,k)+H(k+1,k)*H(k+1,k));28: c(k)=H(k,k)/ro;29: s(k)=H(k+1,k)/ro;30: H(k,k)=ro;31: H(k+1,k)=0;32: for j=k+1:n33: tau=H(k,j);34: H(k,j)=c(k)*tau+s(k)*H(k+1,j);35: H(k+1,j)=-s(k)*tau + c(k)*H(k+1,j);36: end37: end38: for k=1:n-139: for i=1:k+140: tau=H(i,k);41: H(i,k)=tau*c(k)+H(i,k+1)*s(k);42: H(i,k+1)=-tau*s(k)+H(i,k+1)*c(k);43: end44: end45: for i=1:n46: H(i,i)=H(i,i)+miu;47: end


Problema 4.111: function [v ]= vtx1(x)2: %——————————–3: % Fiind dat un vector x nenul de dimensiune n, algoritmul calculeaza4: % un vector v (dimensiune n) astfel incat v T*x=15: % Apel:6: % v=vtx1(x)7: % Stamatescu Grigore, mai 2006.8: %——————————–9: n=length(x);f=0;10: % verificam daca x este nenul11: for i=1:n12: if x(i)==0,13: f=f+1;14: end15: end16: if f==n,17: error(’Vectorul x este nul!’);18: end19: s=0;20: for i=1:n21: s=s+x(i)*x(i);22: end23: ro=-sqrt(s);24: u(1)=ro+x(1);25: for i=2:n26: u(i)=x(i);27: end28: beta=u(1)*ro;29: tau=u(1)*beta;30: v(1)=(1-u(1)*tau)/ro;31: for i=2:n32: v(i)=-u(i)*tau/ro;33: end

Problema 4.121: function [p ]= polinom(A)2: %——————————–3: % Algoritmul calculeaza polinomul caracteristic p al matricei patratice4: % date A (nxn).5: % Apel:6: % p=polinom(A)7: % Stamatescu Grigore, mai 2006.


8: %——————————–9: [n,m]=size(A);10: % analizam daca matricea e patratica11: if n∼=m12: error(’Matricea nu e patratica’);13: end14:15: [V,D]=eig(A);16: for i=1:n17: lambda(i)=D(i,i);18: end19: p(1)=1;20: for i=2:n+121: p(i)=0;22: end23: for i=1:n24: for j=i+1:-1:225: p(j)=p(j)-lambda(i)*p(j-1);26: end27: end

Date post:	19-Jan-2016
Category:	Documents
Upload:	julianmscl
View:	231 times
Download:	4 times

Vectori Proprii Si Valori Proprii

Documents