Bogdan Dumitrescubine precizate de operat˘ii elementare (cum ar adunarea, ^ nmult˘irea, etc.) prin...

Bogdan Dumitrescu

Metode Numericecurs pentru anul 2

v1.0 August 2019

Facultatea de Automatica si CalculatoareUniversitatea Politehnica Bucuresti

Cuvant ınainte

Acest volum contine varianta ideala a cursului de Metode Numerice de laanul 2 Ingineria Sistemelor tinut de mine ın anul scolar 2018-2019. Din di-verse cauze, ıncepand de la pierderea unor ore din motive administrative saude sarbatori legale si terminand cu dorinta mea de a micsora cantitatea deinformatie transmisa, cursul efectiv este cam 85-90% din cel scris.

Cursul nu contine exercitii si probleme pregatitoare pentru examen, aces-tea fiind prezente ın alte lucrari la care studentii au acces prin moodle saupe www.schur.pub.ro.

Pentru ca volumul este doar de uz intern si pur didactic, am eliminatcomplet referintele bibliografice. Lucrari mai ample, unde pot fi gasite maimulte detalii, inclusiv referinte, sunt:

G.H. Golub, C.F. Van Loan, Matrix Computation, Johns Hopkins Univer-sity Press, ed.3, 1996, ed.4, 2013.

B. Dumitrescu, B. Jora, C. Popeea, Calcul numeric matriceal. Algoritmifundamentali, ALL, 1998.

Aceasta fiind prima versiune, studentii care gasesc greseli de orice fel vorfi recompensati cu puncte aditionale atunci cand mi le comunica prin email.Doar prima raportare a unei greseli este premiata. Valoarea premiului depindede importanta greselii.

Voi fi bucuros sa primesc si orice alt fel de sugestii care pot duce laımbunatatirea lucrarii.

Bogdan Dumitrescu august 2019

v

www.schur.pub.ro

Cuprins

1 Specificul calculului numeric si algoritmi matricealielementari . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Reprezentarea ın virgula mobila . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Aritmetica ın virgula mobila . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Conditionarea problemelor de calcul . . . . . . . . . . . . . . . . . . . . . . 10

1.4 Stabilitatea numerica a algoritmilor . . . . . . . . . . . . . . . . . . . . . . . 13

1.5 Calitatile unui algoritm numeric . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.6 Operatii elementare cu vectori si matrice . . . . . . . . . . . . . . . . . . 16

2 Sisteme de ecuatii liniare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.1 Rezolvarea sistemelor triunghiulare . . . . . . . . . . . . . . . . . . . . . . . 22

2.2 Transformari elementare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.3 Triangularizare prin eliminare gaussiana . . . . . . . . . . . . . . . . . . . 28

2.4 Strategii de pivotare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.4.1 Pivotare partiala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.4.2 Pivotare completa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.5 Factorizari LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

2.5.1 Factorizari LU rezultate din eliminarea gaussiana . . . . . 38

2.5.2 Algoritmul Crout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

2.6 Rezolvarea sistemelor liniare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

2.7 Factorizarea Cholesky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

2.8 Calculul determinantului . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

2.9 Calculul inversei (optional) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

2.10 Conditionarea sistemelor liniare . . . . . . . . . . . . . . . . . . . . . . . . . . 56

2.11 Metode iterative pentru rezolvarea sistemelor liniare . . . . . . . . 59

vii

viii Cuprins

3 Problema liniara a celor mai mici patrate . . . . . . . . . . . . . . . . . 65

3.1 Solutia normala a problemei CMMP . . . . . . . . . . . . . . . . . . . . . . 67

3.2 Reflectori Householder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3.3 Triangularizarea ortogonala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

3.4 Rezolvarea problemei CMMP . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.5 Factorizarea QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

3.6 Triangularizare ortogonala cu rotatii . . . . . . . . . . . . . . . . . . . . . . 75

3.7 Rezolvarea problemei CMMP pentru sisteme subdeterminate . 79

3.8 O aplicatie: clasificare binara cu CMMP . . . . . . . . . . . . . . . . . . . 82

4 Valori si vectori proprii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

4.1 Fundamente matematice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

4.2 Forma Schur complexa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.3 Forma Schur reala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

4.4 Metodele puterii si puterii inverse . . . . . . . . . . . . . . . . . . . . . . . . 92

4.4.1 Metoda puterii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

4.4.2 Metoda puterii inverse . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

4.5 Algoritmul PageRank . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

4.6 Algoritmul QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

4.6.1 Asemanare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

4.6.2 Legatura cu metoda puterii . . . . . . . . . . . . . . . . . . . . . . . . 99

4.6.3 Legatura cu metoda puterii inverse . . . . . . . . . . . . . . . . . 100

4.6.4 Alegerea deplasarii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

4.7 Implementarea algoritmului QR . . . . . . . . . . . . . . . . . . . . . . . . . . 101

4.7.1 Aducerea la forma Hessenberg . . . . . . . . . . . . . . . . . . . . . 101

4.7.2 Pas QR explicit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

4.7.3 Testul de oprire si deflatia . . . . . . . . . . . . . . . . . . . . . . . . . 104

4.7.4 Pasul dublu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

4.8 Matrice simetrice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

4.9 Calculul vectorilor proprii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

5 Descompunerea valorilor singulare . . . . . . . . . . . . . . . . . . . . . . . . 109

5.1 Existenta DVS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

5.2 Semnificatii si utilizari ale DVS . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.2.1 Rangul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.2.2 Legatura cu valorile proprii . . . . . . . . . . . . . . . . . . . . . . . . 113

5.2.3 Calculul unor norme matriceale . . . . . . . . . . . . . . . . . . . . 114

Cuprins ix

5.2.4 Numarul de conditionare . . . . . . . . . . . . . . . . . . . . . . . . . . 114

5.2.5 Baze ortogonale pentru subspatii asociate unei matrice 115

5.3 Rezolvarea problemei CMMP generale . . . . . . . . . . . . . . . . . . . . 116

5.4 Problema CMMP totale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

5.5 Calculul DVS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

Capitolul 1

Specificul calculului numeric sialgoritmi matriceali elementari

Rezolvarea numerica—cu ajutorul calculatorului—a unei probleme de calculse face ıntotdeauna pe baza unui algoritm, adica a unei succesiuni finite sibine precizate de operatii elementare (cum ar fi adunarea, ınmultirea, etc.)prin care solutia problemei este calculata ın functie de datele initiale. Desimatematica furnizeaza deseori, ın demonstratiile numite constructive, algo-ritmi de rezolvare a problemelor pe care le trateaza, s-a dovedit ca de multeori implementarile numerice ale acestor algoritmi au dezavantaje majore careıi fac de nefolosit.

Discrepanta ıntre corectitudinea matematica a unui algoritm si adecvarealui la specificul calculului numeric practic poate proveni din doua cauze im-portante. Una, evidenta, consta ın (eventuala) complexitate (ın ce privestenumarul operatiilor de executat si memoria ocupata de date), care poate ducefie la timpi de executie foarte mari, deci nepractici, fie, pur si simplu, la im-posibilitatea executiei programului care implementeaza algoritmul. A douacauza, mai subtila, o constituie erorile care ınsotesc calculul numeric; prinacumulare, acestea pot periclita acuratetea rezultatului.

O serie de erori sunt inerente, cum ar fi cele care apar la determinareadatelor din masuratori experimentale. Altele sunt legate de algoritmul decalcul utilizat, care se poate baza pe simplificari ın modelarea unor fenomenefizice sau poate proveni dintr-o metoda matematica implicand o infinitate deoperatii, din care, prin trunchiere, se retine doar un numar finit de operatii;ın ambele cazuri, prin natura lucrurilor, algoritmul va furniza o solutie aprox-imativa.

In fine, exista erori avand cauze specifice calculului numeric, care opereazacu numere reale (si nu doar cu numere ıntregi); de acestea ne vom ocupa maipe larg. Desigur, numerele reale sunt reprezentate ın memoria calculatoru-lui ıntr-un format bine precizat si ocupand un numar relativ mic de locatiide memorie, de exemplu 4 sau 8 octeti; astfel, o infinitate de numere estereprezentata printr-un numar finit de combinatii de biti. De aici vine o primasursa de erori, cele de reprezentare. Apoi, pot aparea erori si la efectuarea

1

2 1 Specificul calculului numeric si algoritmi matriceali elementari

operatiilor cu numere reale aflate ın memoria calculatorului, erori numitede rotunjire; aceasta nu se ıntampla pentru ca unitatea centrala (micropro-cesorul, de exemplu) ar ”gresi” la calcule, ci, din nou, datorita faptului canumerele reale sunt reprezentate ıntr-un format finit, atat ın memorie, cat siın unitatea centrala.

Prin aceasta, calculul cu numere reale difera fundamental de cel cu ıntregi.Vom obtine 1 + 2 = 3 pe orice calculator, ın schimb 0.1 + 0.2 6= 0.3 pe oricecalculator (doar daca nu se folosesc programe speciale!); de exemplu — pe uncalculator si ıntr-un limbaj de programare pe care nu le precizam — adunand0.1 + 0.2 si scazand din rezultat 0.3 nu obtinem 0, asa cum ne-am astepta, ciaproximativ 5.5 · 10−17; e drept, eroarea este infima. Pericolul apare ın cazulacumularii unor astfel de mici erori, acumulare care poate duce la degradarea—uneori fatala—a rezultatului produs.

1.1 Reprezentarea ın virgula mobila

Fie x si x numere reale, x fiind interpretat ca o aproximare a lui x. Calitateaaproximarii poate fi cuantificata prin urmatoarele doua masuri.

Eroarea absoluta (cu care x aproximeaza x) se defineste prin

∆ = |x− x|.

Daca x 6= 0, atunci eroarea relativa se defineste prin

ε =

∣∣∣∣x− xx∣∣∣∣ =

∆

|x|.

Daca x ∈ Rn, se ınlocuieste ın relatiile de mai sus valoarea absoluta | · | cuo norma vectoriala ‖ · ‖.

Exemplul 1.1 Fie x = 1.0, si x = 0.999 o aproximare a sa. Atunci ∆ = 10−3

si ε = 10−3. Daca y = 0.009 este o aproximatie a lui y = 0.01, atunci eroareaabsoluta este aceeasi ca ın cazul precedent, ∆ = 10−3, dar eroarea relativaeste de o suta de ori mai mare: ε = 10−1. Raportandu-se la valoarea lui x,eroarea relativa este mult mai adecvata pentru estimarea calitatii aproximarii.

Erorile de reprezentare apar datorita memorarii ın calculator a numerelorreale printr-o secventa finita de simboluri (cifre binare). Pentru a prezenta oestimare a acestor erori, reamintim bine cunoscuta reprezentare pozitionalaa numerelor reale. Fie

• β ∈ N, β ≥ 2, baza de numeratie;

1.1 Reprezentarea ın virgula mobila 3

• C = {0, 1, . . . , β− 1}, multimea cifrelor ın baza β, adica primele β numerenaturale.

Se stie ca orice numar x ∈ R poate fi scris unic sub forma unei secventeinfinite

x = san−1an−2 . . . a0.a−1a−2a−3 . . . , (1.1)

care nu se termina printr-o secventa infinita de cifre egale cu β − 1 si ın careai ∈ C, iar s este semnul, s ∈ {+,−}. Valoarea lui x este

x = s

(n−1∑i=0

aiβi +

∞∑i=1

a−iβ−i

). (1.2)

Convenim sa eliminam din scriere secventa infinita de zerouri finale, atuncicand este cazul.

Exemplul 1.2 Numarul ın baza 10

3.25 = 3 · 100 + 2 · 10−1 + 5 · 10−2

se reprezinta ın baza 2 ın modul urmator (verificati egalitatea):

11.01 = 1 · 21 + 1 · 20 + 0 · 2−1 + 1 · 2−2.

In ambele baze, reprezentarile au un numar finit de cifre. In schimb, numarulexprimat simplu ın baza 10:

0.1 = 1 · 10−1

are o reprezentare de lungime infinita ın baza 2 (verificati din nou egalitatea):

0.0001100110011 . . . = (2−4 + 2−5)

∞∑k=0

2−4k.

Asadar, numarul 0.1 nu se poate reprezenta exact ıntr-un format finit.

Evident, reprezentarea numerelor reale pe calculator poate avea doar unnumar finit de cifre si deci, prin natura ei, este aproximativa. O ”buna”aproximare printr-o secventa finita de lungime fixata trebuie sa asigure:

• un domeniu suficient de mare de numere reprezentate;

• o eroare relativa de reprezentare suficient de mica;

• o distributie uniforma a erorii relative de reprezentare.

Un bun compromis ın satisfacerea acestor cerinte este dat de FormatulVirgula Mobila (FVM).


Definitia 1.3 Un Format Virgula Mobila (FVM) este definit de trei ıntregi(β, t, p), cu urmatoarea semnificatie:

β – baza de numeratie (β ≥ 2);

t – precizia, adica numarul de cifre semnificative (”lungimea” mantisei);

p – numarul de cifre ale exponentului.

Un numar ın virgula mobila este o pereche (f, e), unde

f = ±0.f1f2 . . . ft, f1 6= 0, fi ∈ C

este un numar fractionar cu semn (plus sau minus), normalizat (f1 6= 0), cut cifre, denumit mantisa, iar

e = ±e1e2 . . . ep

este un ıntreg cu semn, cu p cifre, numit exponent. Valoarea reala asociatanumarului ın virgula mobila (f, e) este

x = f · βe.

Fie L valoarea minima admisa a exponentului (de exemplu L = − ρ . . . ρ︸︷︷︸p

,

unde ρ = β− 1) si U cea maxima (de exemplu U = −L). Se observa ca toatenumerele ın virgula mobila x satisfac:

βL−1 = m ≤ |x| ≤M = βU (1− β−t)

adica domeniul reprezentarii ın virgula mobila este intervalul R = [−M,M ].Desigur, pentru majoritatea reprezentarilor pe calculator, β = 2.

Exemplul 1.4 Consideram FVM cu (β, t, p) = (2, 2, 1). Atunci L = −1,U = 1, m = 2−2 = 0.2510, M = 2(1 − 2−2) = 1.510. Numerele reale care auo reprezentare exacta ın acest FVM sunt:

x1 = (−0.11,+1) = −1.510 x7 = (+0.10,−1) = 0.2510x2 = (−0.10,+1) = −1.010 x8 = (+0.11,−1) = 0.37510x3 = (−0.11, 0) = −0.7510 x9 = (+0.10, 0) = 0.510x4 = (−0.10, 0) = −0.510 x10 = (+0.11, 0) = 0.7510x5 = (−0.11,−1) = −0.37510 x11 = (+0.10,+1) = 1.010x6 = (−0.10,−1) = −0.2510 x12 = (+0.11,+1) = 1.510

si vor fi folosite pentru reprezentarea (aproximativa) a tuturor numerelorreale din intervalul [−1.5, 1.5]. Reprezentarea geometrica a tuturor acestornumere pe axa reala este prezentata ın Figura 1.1. Numerele cu reprezentareexacta ın virgula mobila sunt relativ uniform distribuite; mai precis, |xi −xi−1| / |xi| are aproximativ aceleasi valori pentru orice i.

1.1 Reprezentarea ın virgula mobila 5

1

x11

-x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x12

1.50.5-0.5-1-1.50

Fig. 1.1 Reprezentarea numerelor ın virgula mobila din exemplul 1.4

0 1 1 0 1 0 0 1

-�exponent (7 biti)

0 1 0 0 1 1 0 0 1 0 1 1 1 0 0 0 0 0 0 0 1 0 1 1

mantisa (23+1 biti)6

semn numar (mantisa)

� -6

semn exponent

Fig. 1.2 Structura unui numar ın virgula mobila, pentru (β, t, p) = (2, 24, 7)

Exemplul 1.5 Un FVM pe 32 biti poate fi urmatorul: (β, t, p) = (2, 24, 7).Doi biti sunt atribuiti reprezentarii semnelor mantisei si exponentului; deoa-rece primul bit al mantisei are ıntotdeauna valoarea 1, semnul mantisei sememoreaza ın locul acestui bit, economisind astfel un bit al reprezentarii. Unexemplu de numar ın virgula mobila ın acest format este prezentat ın Figura1.2.

Pentru acest FVM avem: L = −127, U = 127, m = 2−127 = (2−10)1323 ≈10−38, M ≈ 1038 (putem aproxima 210 ≈ 103). Aceasta arata ca FVMacopera un domeniu larg de numere, utilizand un numar relativ mic de cifrebinare. Acest exemplu este foarte apropiat de formatul standard IEEE pentrunumere ın virgula mobila ın simpla precizie.

In formatul din acest exemplu, numarul real 0.1 are mantisa

0.110011001100110011001100

si exponentul 10000011 (−3 ın baza 10; prima cifra 1 a exponentului reprezintasemnul). Asadar numarul 0.1 este reprezentat aproximativ (vezi exemplul1.2).

Observatia 1.6 Definitia 1.3 nu permite reprezentarea exacta a numarului0, acesta neputand avea o mantisa normalizata. Practic, pentru a reprezenta0, precum si unele simboluri speciale despre care se va vorbi mai tarziu, sefolosesc valori dedicate ale exponentului, ın afara intervalului [L,U ]. De aceea,ın standardul IEEE pentru simpla precizie, L = −126.

Pentru a estima erorile efectuate la reprezentarea unui numar real ın vir-gula mobila, consideram un FVM, adica (β, t, p), si notam

F = {x ∈ R | x are reprezentare exacta in VM} ∪ {0}. (1.3)

Desigur, F este o multime finita de numere rationale.


Reprezentarea numarului real x ∈ [−M,M ] ın FVM ınseamna aproxi-marea lui cu un numar x ∈ F . Aceasta aproximare poate fi exprimata printr-ofunctie de rotunjire.

Definitia 1.7 Daca (β, t, p) este un FVM si F este multimea definita ın(1.3), atunci o functie

fl : [−M,M ]→ F

care asociaza fiecarui x ∈ [−M,M ] o unica reprezentare ın virgula mobila

x = fl(x)

este numita functie de rotunjire. Eroarea relativa de aproximare

|x− fl(x)||x|

definita pentru orice x ∈ [−M,M ] nenul este denumita eroare de reprezentare.

Deoarece intervalul [−M,M ] este o multime infinita de numere reale,fiecare x ∈ F constituie ”reprezentarea ın VM” a unei multimi infinite denumere din [−M,M ]; suntem interesati sa gasim o margine superioara aerorii de reprezentare pentru o functie de rotunjire data. Exista mai multefunctii de rotunjire. O vom prezenta doar pe cea mai simpla, denumita ro-tunjire prin taiere. In acest scop, scriem numarul x ∈ [−M,M ] ın forma (1.1)normalizata:

x = f · βe = ±0.f1f2 . . . ftft+1 . . . · βe =

= ±0.f1f2 . . . ft · βe ± 0.ft+1ft+2 . . . · βe−t =

= f · βe + g · βe−t,

unde fi ∈ C, f1 6= 0, f = ±0.f1f2 . . . ft, g = ±0.ft+1ft+2 . . .

In mod evident:

1/β ≤ |f | < 1, 1/β ≤ |f | < 1, 0 ≤ |g| < 1. (1.4)

Functia de rotunjire prin taiere

fl1 : [−M,M ]→ F

este definita prin

x = fl1(x) =

{f · βe, pentru x 6= 0,

0, pentru x = 0.

1.2 Aritmetica ın virgula mobila 7

Pe scurt, reprezentarea ın VM se obtine prin taierea cifrelor mantisei nor-malizate a numarului x care se afla ın afara formatului (de la pozitia t + 1ıncolo).

Utilizand inegalitatile (1.4), este usor de stabilit o margine superioaraa erorii de reprezentare introduse de fl1. Intr-adevar, pentru orice x ∈[−M,M ] \ {0} avem

ε =|x− fl1(x)||x|

=|fβe − fβe||f |βe

=|g|βe−t

|f |βe<β−t

β−1= β−t+1.

Aceasta formula arata ca, indiferent de valoarea numarului x, o marginepentru marimea erorii de reprezentare este determinata exclusiv de numarulde cifre ale mantisei (ın baza de numeratie a FVM) si acesta este motivulpentru care t este numit precizia reprezentarii ın virgula mobila. Numarulβ−t+1 este numit epsilon masina si reprezinta distanta dintre 1 si urmatorulnumar cu reprezentare ın FVM.

Standardul IEEE este implementat pe practic toate procesoarele din cal-culatoarele de uz general. Acesta prevede doua FVM, anume simpla si dublaprecizie. In C, tipurile de variabile corespunzatoare se numesc float, respectivdouble. In MATLAB, variabilele reale sunt implicit ın dubla precizie. Nu-merele ın dubla precizie sunt memorate pe 64 biti, cu t = 53 si p = 10 (unbit suplimentar este dedicat memorarii semnului exponentului).

Cele t cifre ale mantisei sunt numite cifre semnificative ale numaruluireal reprezentat. In baza de numeratie zecimala, numerele reprezentate ınsimpla precizie conform standardului IEEE (t = 24) au 7 cifre semnificative(numaram 3 cifre zecimale pentru 10 cifre binare), iar cele ın dubla pre-cizie (t = 53) au 16 cifre semnificative; ın consecinta, eroarea maxima dereprezentare este de aproximativ 10−7, respectiv 10−16 (anume chiar epsilonmasina).

1.2 Aritmetica ın virgula mobila

Procesoarele cu unitati dedicate calculelor ın FVM implementeaza operatiielementare cu numere FVM ca, de exemplu, adunarea, scaderea, ınmultireasau ımpartirea. Trebuie subliniat ca operatiile mentionate difera de celematematice corespunzatoare deoarece rezultatul lor apartine ıntotdeaunamultimii finite F a numerelor ın FVM. Asadar operatiile aritmetice nu potfi executate exact. Eroarea introdusa de operatiile aritmetice ın virgula mo-bila se numeste eroare de rotunjire. Consecintele pot fi foarte importante;exista posibilitatea ca erorile de rotunjire sa se acumuleze si valorile finalesa fie afectate, chiar pana la a le face inutilizabile. De aceea este de dorit ca


fiecare algoritm sa fie analizat pentru a vedea ın ce masura erorile de rotunjireafecteaza rezultatele; ın general, aceasta este o sarcina dificila.

Exemplul 1.8 Consideram FVM cu (β, t, p) = (10, 3, 1). Adunand exactnumerele 1.23 si 20.5, amandoua reprezentabile exact ın FVM ales, se obtine21.73. Rotunjit prin taiere, rezultatul ın VM este 21.7, deci diferit de celexact.

Analiza erorilor de rotunjire se bazeaza pe cunoasterea marginilor supe-rioare ale erorilor ce apar ın operatiile ın virgula mobila. Standardul IEEEasigura ca eroarea relativa maxima la efectuarea unei operatii aritmetice esteidentica cu eroarea de reprezentare a rezultatului exact al operatiei. De ex-emplu, presupunand ca trebuie calculata suma numerelor x, y ∈ F , rezultatulcalculat este fl(x+y). Observam ca acesta este de fapt cel mai precis rezultat(ın FVM) la care putem spera. (Notam ın plus ca modul de calcul e nebanal,deoarece procesoarele nu calculeaza efectiv suma exacta x+ y.)

Atunci cand o operatie ın FVM produce un numar cu un modul prea mare(> M), apare eroarea numita de depasire superioara (overflow). In standardulIEEE pentru virgula mobila, exista o valoare speciala, numita Inf (infinit),care este atribuita rezultatului ın cazul unei depasiri superioare; calculelecontinua cu aceasta valoare; rezultatele se obtin conform regulilor uzualede lucru cu infinitati. O alta valoare speciala—NaN (Not a Number)—esteatribuita rezultatelor nedefinite, ca 0 · ∞, 0/0, ∞/∞; o operatie implicandNaN are ıntotdeauna ca rezultat NaN.

Asadar, ın afara depasirilor, o operatie aritmetica introduce erori relativemici, de ordinul β−t (adica al lui epsilon masina). O problema fundamentalaa calculului numeric este evaluarea marimii erorii ce afecteaza rezultatul ıncazul unei secvente de operatii.

Pentru a aborda problema, sa consideram doua exemple ın care apar erorinumerice mari, dar din cauze esential diferite. Lucram ın FVM cu (β, t, p) =(10, 3, 1). Pentru a evita confuziile, vom nota cu ⊕, , ⊗ si � adunarea,scaderea, ınmultirea, respectiv ımpartirea ın VM; deci x+y este suma exacta,iar x⊕ y = fl(x+ y) este suma calculata ın VM.

Exemplul 1.9 In calculul radacinilor polinomului de gradul al doilea ax2 +bx+c, cu a 6= 0, este necesar calculul expresiei b2−4ac. Considerand b = 3.34,a = 1.22, c = 2.28, avem b2− 4ac = 0.0292, ın timp ce, rotunjind prin taiere,b ⊗ b = 4 ⊗ a ⊗ c = 11.1, deci b ⊗ b − 4 ⊗ a ⊗ c = 0. Rezultatul calculatare toate cifrele semnificative eronate, iar eroarea relativa aferenta este egalacu 1; totusi, pentru fiecare operatie ın parte, eroarea relativa este mai micadecat 10−2.

Exemplul 1.10 Daca polinomul de gradul II are radacini reale, acestea secalculeaza de obicei utilizand formulele

1.2 Aritmetica ın virgula mobila 9

x1 =−b−

√b2 − 4ac

2a, x2 =

−b+√b2 − 4ac

2a. (1.5)

Luand b = 10.1, a = 0.0123, c = 32.4, valoarea exacta rotunjita la 3cifre semnificative a lui x2 este −3.22. Efectuand calculele ın formatul ales,obtinem ∆ = b ⊗ b 4 ⊗ a ⊗ c = 100 (ın loc de 100.41, dar eroarea rel-ativa e ınca de ordinul 10−3); dar −b ⊕ SQRT (∆) = −0.1, ın timp ce−b+

√b2 − 4ac = −0.0792, eroarea relativa fiind acum 0.26. Valoarea calcu-

lata a radacinii x2 este −4.06, ın timp ce valoarea exacta este −3.22; din noutoate cifrele semnificative ale rezultatului sunt eronate.

Se observa ca, ın ambele exemple, erorile mari apar dupa scaderea a douanumere de valori apropiate. Se produce o anulare catastrofala, adica o pierderecompleta a cifrelor semnificative ale rezultatului. In primul caz anularea tre-buie pusa integral pe seama operatiei matematice de scadere; problema cal-culului b2−4ac este cauza erorilor, si nu algoritmul de calcul. In schimb, ın aldoilea caz, algoritmul de calcul al radacinilor este de vina, deoarece nu evitascaderea a doua numere apropiate, desi aceasta este posibil, dupa cum vomvedea putin mai tarziu. Sa studiem ıntai mai amanuntit ce se ıntampla lascaderea a doua numere apropiate.

Exemplul 1.11 Consideram doua numere reale, x si y, de acelasi ordin demarime (cu x > y, pentru simplitate), ale caror prime k cifre coincid, cuurmatoarele reprezentari (exacte) ın baza β (d1 6= 0):

x = 0.d1d2 . . . dkdk+1 . . . dt1 · βe,y = 0.d1d2 . . . dkck+1 . . . ct0 · βe.

Rotunjind prin taiere, reprezentarile lor ıntr-un FVM cu (β, t, p) sunt

x = fl(x) = 0.d1d2 . . . dkdk+1 . . . dt · βe,y = fl(y) = 0.d1d2 . . . dkck+1 . . . ct · βe.

Diferenta lor, calculata ın VM, este

x y = 0.f1 . . . ft−k0 . . . 0 · β−k+e.

Observam ca x− y = x y, deci scaderea se efectueaza fara erori. In schimb,rezultatul scaderii cu valorile exacte ale operanzilor este

x− y = 0.f1 . . . ft−k10 . . . 0 · β−k+e,

iar eroarea relativa asociata lui este∣∣∣∣ (x− y)− (x− y)

x− y

∣∣∣∣ ≥ β−t+1+e

β−k+e= β−t−1+k.


In concluzie, observam ca operanzii scaderii sunt afectati de erori relative dereprezentare de ordinul β−t−1, ın timp ce rezultatul ei are o eroare relativade ordinul β−t−1+k; atunci cand x si y au acelasi ordin de marime (deci multecifre semnificative identice, deoarece k este aproape de t), eroarea asociatarezultatului scaderii poate deveni uriasa.

Exemplul 1.11 arata ca adunarea (scaderea) nu introduce erori mari prinea ınsasi, ci potenteaza—ın anumite cazuri—erorile operanzilor.

De aceea, scaderea numerelor apropiate trebuie evitata. Pentru aceasta,daca este posibil, algoritmul de calcul este cel care trebuie modificat.

Exemplul 1.12 Am constatat, ın exemplul 1.10, ca formulele (1.5) nu suntcele mai bune din punct de vedere numeric, mai ales cand b2 � 4ac, sideci |b| ≈

√∆. Pentru a evita scaderea unor numere apropiate, se pot folosi

formulele

x1 =−b− sgn(b)

√b2 − 4ac

2a, x2 =

−2c

b+ sgn(b)√b2 − 4ac

. (1.6)

Utilizand ın (1.6) valorile numerice din exemplul 1.10, se obtine, ın FVMde lucru, un rezultat cu toate cifrele semnificative exacte.

1.3 Conditionarea problemelor de calcul

Orice problema de calcul numeric impune obtinerea unor rezultate numerice,pornind de la date de intrare numerice, astfel ıncat rezultatele sa poata fiinterpretate ca o aproximatie a solutiei unei probleme matematice pentrudatele de intrare respective. Pentru ca o astfel de problema sa fie bine definita,trebuie asigurate existenta si unicitatea solutiei. In acest caz, o problema decalcul numeric poate fi ıntotdeauna descrisa ca evaluare a unei functii

f : D ⊂ Rn −→ Rm (1.7)

ıntr-un punct x ∈ D dat. Cele n componente ale argumentului constituiedatele de intrare, iar cele m componente ale vectorului f(x) sunt rezultatelesau datele de iesire.

Desigur, functia f este de cele mai multe ori definita implicit iar de-scoperirea unei cai prin care calculul valorii f(x) devine posibil echivaleazacu rezolvarea, printr-un algoritm specific, a problemei de calcul considerate.

Exemplul 1.13 a) Fie D o multime de trei numere reale si

1.3 Conditionarea problemelor de calcul 11

x =

abc

∈ D = {x ∈ R3 | a 6= 0, b2 − 4ac ≥ 0} ⊂ R3.

Rezolvarea ecuatiei de gradul II

ay2 + by + c = 0

este o problema numerica bine definita, deoarece pentru orice x ∈ D existaun unic vector

y = f(x) =

[y1y2

]=

−b− sgn(b)

√b2 − 4ac

2a

−2c

b+ sgn(b)√b2 − 4ac

∈ R2,

care este solutia problemei.

b) Calculul integralei ∫ b

a

e−t2

dt,

unde a, b ∈ R sunt date (a 6= b), este de asemenea o problema numerica bine

definita. Intr-adevar, pentru datele de intrare x =

[ab

]∈ D = R2, exista si

este unica valoarea integralei y ∈ R, desi integrala nedefinita∫

e−t2

dt nu esteexprimabila prin functii elementare.

Abaterea solutiei calculate fata de solutia exacta a unei probleme numericepoate avea doua surse esential diferite, legate de problema ın sine si, re-spectiv, de algoritmul de calcul utilizat. De prima sursa de erori—numitaconditionarea problemei—ne ocupam ın continuare, de cea de-a doua—stabilitatea algoritmului—ın sectiunea urmatoare.

Conditionarea problemei caracterizeaza sensibilitatea solutiei ın raport cuperturbatiile datelor de intrare. In practica este cunoscuta doar o aproximatiex a datelor de intrare x si, ın cel mai bun caz, presupunand ca toate calculelese efectueaza exact, putem calcula f(x) ın loc de f(x). Daca functia f definindproblema noastra este prea sensibila la variatii ale argumentului, f(x) si f(x)pot diferi semnificativ chiar ın cazul unei erori relative mici ‖x − x‖ / ‖x‖ adatelor de intrare. O astfel de problema este numita rau conditionata.

Daca ıncercam rezolvarea unei probleme rau conditionate pornind de ladate inexacte, atunci solutia obtinuta va fi ınsotita de erori mari indiferentde felul ın care este calculata!

Definitia 1.14 Conditionarea numerica sau sensibilitatea locala a unei prob-leme f : D ⊂ Rn → Rm ıntr-un punct x ∈ D este exprimata de amplificareaerorii relative


κ(x) =‖f(x)− f(x)‖‖f(x)‖

‖x‖‖x− x‖

, (1.8)

pentru x 6= 0 si f(x) 6= 0. Cand κ(x) este mic (de ordinul unitatii, de exem-plu) problema este numita bine conditionata ın punctul x. Altfel, problemaeste rau conditionata.

Subliniem faptul ca conditionarea este o proprietate locala. De obicei,totusi, o problema este numita bine conditionata daca este bine conditionataın orice punct x ∈ D; de asemenea, este numita rau conditionata, daca existapuncte x ∈ D ın care este rau conditionata.

Exemplul 1.15 Problema calculului sumei a doua numere x, y ∈ R poatefi bine sau rau conditionata functie de valorile lui x si y (datele de intrare).(Aici f(x, y) = x + y, f : R2 −→ R; convenim sa scriem si f(x, y) ın loc de

f(

[xy

]), cum ar cere (1.7).)

Presupunem ca x, y sunt aproximari ale lui x, respectiv y, afectate de erorirelative de acelasi ordin de marime, γ; mai precis,

x = x+ εxx, cu |εx| ≤ γ � 1,y = y + εyy, cu |εy| ≤ γ.

Eroarea relativa asociata sumei este (calculul sumei se face exact)

ε =|x+ εxx+ y + εyy − (x+ y)|

|x+ y|=|εxx+ εyy||x+ y|

. (1.9)

Daca x si y au acelasi semn, atunci

ε ≤ |εxx|+ |εyy||x|+ |y|

≤ γ

si suma este afectata de o eroare de acelasi ordin de marime ca si operanzii;ın acest caz, adunarea este bine conditionata.

Daca x si y au semne diferite, este posibil (atunci cand εx si εy au si elesemne diferite) ca ε si

θ =γ ·max(|x|, |y|)|x+ y|

sa aiba acelasi ordin de marime. Dar daca |x+ y| � max(|x|, |y|), ceea ce seıntampla atunci cand x si y au valori apropiate, atunci θ � γ. In acest caz,problema calculului sumei este rau conditionata. Exemplele 1.9, 1.10, 1.11ilustreaza aceasta situatie.

1.4 Stabilitatea numerica a algoritmilor 13

1.4 Stabilitatea numerica a algoritmilor

Pentru a rezolva o problema de calcul numeric definita de functia (1.7), uncalculator executa o secventa (finita) bine precizata de operatii, denumitaalgoritm. In general exista mai multi algoritmi pentru rezolvarea unei prob-leme. Un algoritm f pentru problema f produce un raspuns unic determinatpentru date de intrare fixate si deci poate fi exprimat matematic printr-ofunctie

f : D ∩ F ⊂ Rn −→ Rm.

Deoarece executia unui algoritm se face folosind numere ın FVM si calculelesunt ınsotite de erori, este clar ca f(x) este de regula diferit de f(x). De

asemenea, nu ne putem astepta ca algoritmul f sa calculeze solutia uneiprobleme rau conditionate mai precis decat o garanteaza datele de intrare.Totusi, e de dorit ca f sa nu introduca erori exagerat de mari. Notiuneade stabilitate numerica exprima marimea erorilor numerice introduse de unalgoritm, ın ipoteza ca datele initiale sunt exacte.

Definitia 1.16 Un algoritm f , destinat rezolvarii problemei f , este numericstabil daca una dintre urmatoarele doua conditii alternative este ındeplinita:

i) f(x) este aproape de f(x) pentru orice intrare x, adica solutia calculataaproximeaza bine solutia exacta.

ii) pentru orice intrare x, exista un x ∈ D apropiat de x astfel ıncat

f(x) sa fie egal cu f(x), adica solutia calculata de algoritmul f cu date deintrare exacte este egala cu solutia exacta pentru niste date de intrare usorperturbate.

Altfel, algoritmul este numeric instabil.

Cele doua conditii sunt similare, doar ca prima se formuleaza ın spatiuldatelor de iesire, iar a doua ın spatiul datelor de intrare, dupa cum se poateremarca ın Figura 1.3; ın mod plastic, eroarea (absoluta) f(x) − f(x) estenumita eroare ”ınainte”, deoarece corespunde sensului de calcul al algorit-mului, iar x− x eroare ”ınapoi”, deoarece gasirea unui x care sa satisfaca ii)ınseamna ıntoarcerea ın spatiul datelor de intrare, ulterior aplicarii algorit-mului. Analiza stabilitatii numerice a unui algoritm ınseamna cuantificareaerorilor din definitia 1.16; evaluarea ‖f(x)−f(x)‖ este caracteristica analizeiınainte (directe); evaluarea ‖x− x‖ defineste analiza ınapoi (inversa).

Notiunea vaga ”aproape” se poate reexprima prin cea de eroare relativamica (notiune destul de vaga si ea); ın calculul numeric se folosesc deseoriastfel de exprimari. In mod sigur, o eroare relativa de ordinul lui epsilonmasina caracterizeaza un algoritm numeric stabil; solutia este calculata cuprecizia intrinseca a FVM utilizat.


x

x

Eroare ”ınapoi”

f

f

f

f(x)

Eroare ”ınainte”

f(x) = f(x)

Fig. 1.3 Erori ınainte si ınapoi

Este usor de aratat ca utilizarea unui algoritm stabil la rezolvarea uneiprobleme bine conditionate conduce la un rezultat de ınalta acuratete1.

In mod similar, se poate observa ca, ın alte situatii (de exemplu pentruprobleme rau conditionate sau pentru algoritmi instabili), marimea erorii aso-ciate rezultatului nu poate fi apreciata. Cand un algoritm stabil este folositpentru a rezolva o problema rau conditionata, nu exista nici o garantie caf(x) si f(x) au valori apropiate si, desi f(x) este (aproximativ) egal cu f(x),

diferenta dintre f(x) si f(x) poate fi semnificativa. Cand un algoritm instabileste folosit pentru rezolvarea unei probleme bine conditionate, solutia calcu-lata f(x) pentru date exacte de intrare este solutie exacta a problemei cu xdeparte de x, si tot asa este f(x) de f(x). In sfarsit, nu e nici o sperantade a obtine erori mici atunci cand un algoritm instabil este utilizat pentru oproblema rau conditionata.

In aceasta lucrare vom fi interesati mai mult de prezentarea celor mai bunialgoritmi care rezolva o problema, si mai putin de demonstrarea—aproapeıntotdeauna extrem de ”tehnica”—a stabilitatii lor numerice.

1.5 Calitatile unui algoritm numeric

Pentru a ıncheia discutia generala despre algoritmi, precizam criteriile dupacare se apreciaza calitatea unui algoritm numeric.

• Numarul de operatii. Timpul de executie al unui program—pe un cal-culator precizat—este o caracteristica de prim ordin a acestuia. De regula,un program se executa de mai multe ori, pentru a rezolva o anume prob-lema pentru diverse date de intrare; timpul de executie permite apreciereanumarului de solutii care se pot obtine ıntr-un anumit timp, de obicei impus.

1 Atragem atentia asupra sensurilor diferite pe care le au ın calculul numeric termenii

precizie si acuratete, aproape sinonimi ın limbajul uzual. Precizia se refera la numarul decifre semnificative din FVM, iar acuratetea la marimea erorii dintre rezultatul calculat

numeric si cel exact.

1.5 Calitatile unui algoritm numeric 15

Desigur, este de dorit a estima timpul de executie al unui program ınaintede crearea sa efectiva, numai pe baza algoritmului implementat de program.O buna masura este numarul de operatii aritmetice efectuate; ın algoritmiinumerici ele constituie de obicei partea majoritara a operatiilor, cele cu ıntregisau logice avand o pondere mica, atat ca numar, cat si ca durata de executie.

Numim flop o operatie aritmetica elementara (adunare, scadere, ınmultire,ımpartire). Numarul de operatii al unui algoritm este numarul de flopi nece-sari obtinerii rezultatului. Trebuie precizat ca, desi numarul de operatii poateindica doar cu aproximatie timpul de executie pe un calculator dat, ın schimbeste o foarte buna masura pentru compararea a doi algoritmi.

De regula, numarul de operatii Nop al unui algoritm e functie de dimen-siunea problemei care, de multe ori, poate fi apreciata prin numarul datelorde intrare, Nop = f(n). In calculul numeric, de obicei f este un polinom(de grad mic), f(n) = a0n

k + a1nk−1 + . . .+ ak. Pentru ca f poate avea o

expresie complicata, se pastreaza doar termenul cel mai semnificativ si sescrie Nop ≈ a0n

k. O astfel de aproximare este cu atat mai buna cu cat neste mai mare. Pentru aprecieri pur calitative, vom omite a0 si vom scrie2

Nop = O(nk).

• Memoria ocupata. Pentru executia unui algoritm (program) este nece-sara memorarea datelor de intrare, rezultatelor, precum si a altor valori in-termediare. Numarul de elemente ın FVM necesare ın acest scop se numestememoria ocupata de algoritm.

In majoritatea cazurilor, datele de intrare nu mai sunt necesare la ter-minarea algoritmului. Spatiul de memorie ocupat de ele poate fi utilizat peparcursul executiei algoritmului pentru a memora rezultate partiale si/saufinale. In acest caz se spune ca efectuam calculele pe loc.

• Stabilitatea numerica, discutata ın paragraful anterior, arata acuratetearezultatelor obtinute de algoritmul numeric (ın ipoteza ca datele initiale suntexacte).

Aceste trei criterii trebuie corelate cu viteza calculatorului, memoriadisponibila, numarul de cifre semnificative al FVM utilizat.

Desigur, nu trebuie uitata nici un moment conditionarea problemei, careinfluenteaza si ea acuratetea rezultatului. Intereseaza ın special performantelealgoritmului pentru probleme rau conditionate, adica atunci cand eroriledatelor de intrare sunt mult amplificate. In orice caz, calculul solutiei e binea fi ınsotit de determinarea (sau estimarea) conditionarii problemei, adicade calculul unei cantitati de tip (1.8). Putem adauga astfel un al patruleacriteriu de apreciere a unui algoritm:

• Siguranta ın functionare – capacitatea de a semnala situatiile ın carerezultatul poate fi afectat de erori importante datorate relei conditionari a

2 Se scrie f(n) = O(g(n)) daca exista doua constante c1, c2 astfel ıncat c1g(n) ≤ f(n) ≤c2g(n).


problemei. Decizia de a utiliza sau nu un astfel de rezultat revine utilizatoruluisau unui program expert.

1.6 Operatii elementare cu vectori si matrice

Discutam aici algoritmii la nivel de element pentru cateva operatii simple cuvectori si matrice.

Produsul scalar. Dandu-se doi vectori x, y ∈ Rn, produsul lor scalar este

yTx =

n∑i=1

xiyi. (1.10)

Un vector este implicit coloana (adica este organizat pe verticala). Notam yT

vectorul linie obtinut prin transpunerea lui y.

Algoritmul 1.1 (DOT – Calculul produsului scalar) (Se dau vectoriix, y ∈ Rn. Se calculeaza α = yTx.)

1. α← 02. Pentru i = 1 : n

1. α← α+ xiyi

Comentarii. Algoritmul necesita 2n operatii. Il vom apela prin α =DOT(x, y).

Produsul matrice-vector. Dandu-se o matrice A ∈ Rm×n si un vectorx ∈ Rn, produsul y = Ax ∈ Rm este definit prin

yi = aTi x =

n∑j=1

aijxj , (1.11)

unde aTi este linia i a matricei A. Produsul matrice-vector este ilustrat ınFigura 1.4. Algoritmul care calculeaza produsul este urmatorul.

Algoritmul 1.2 (PMV – Produs matrice-vector) (Se dau A ∈ Rm×n,x ∈ Rn. Se calculeaza y = Ax ∈ Rm utilizand operatii elementare.)

1. Pentru i = 1 : m1. yi ← 02. Pentru j = 1 : n

1. yi ← yi + aijxj

Se observa imediat ca bucla interioara reprezinta o operatie DOT.

Comentarii. Algoritmul necesita∑mi=1

∑nj=1 2 = 2mn flopi.

1.6 Operatii elementare cu vectori si matrice 17

A x

·

Fig. 1.4 Produsul matrice vector: produse scalare ıntre liniile matricei si vector.

xT y

· =

c x yT

· =

C

Fig. 1.5 Produs scalar (stanga) si produs exterior (dreapta).

Produsul matriceal. Fie A ∈ Rm×p si B ∈ Rp×n, doua matrice; produsullor este matricea C = AB ∈ Rm×n, definita prin

cij =

l∑k=1

aikbkj , i = 1 : m, j = 1 : n.

Discutam ıntai cazurile particulare ın care cel putin una dintre dimensiunieste egala cu 1.

Daca m = n = 1, atunci Adef= xT este un vector linie, B

def= y este

un vector coloana, ambii ın Rp, iar produsul lor coincide cu simplul produsscalar AB = xT y.

Daca p = 1, atunci Adef= x este un vector coloana ın Rm, B

def= yT este un

vector linie ın Rn, iar produsul lor este matricea C = xyT ∈ Rm×n, definitaprin cij = xiyj ; aceasta operatie cu doi vectori poarta numele de produsexterior. Produsul scalar si cel exterior sunt ilustrate ın Figura 1.5.

Daca n = 1, atunci Bdef= y este un vector coloana si operatia AB este o

ınmultire matrice-vector.

Daca m = 1, atunci Adef= xT este un vector linie si AB = xTB este un

vector linie (ınmultire vector linie - matrice).

Aplicand direct definitia, se obtine algoritmul urmator.

Algoritmul 1.3 (MM – Inmultire matriceala) (Se dau A ∈ Rm×p,B ∈ Rp×n. Se calculeaza C = AB ∈ Rm×n.)


1. C ← 02. Pentru i = 1 : m

1. Pentru j = 1 : n1. Pentru k = 1 : p

1. cij ← cij + aikbkj

Comentarii. Algoritmul necesita∑mi=1

∑nj=1

∑pk=1 2 = 2mnp flopi. Ob-

servam ca ordinea celor trei bucle poate fi schimbata ın orice mod, fara aschimba rezultatul sau numarul de operatii. In forma ijk de mai sus, ma-tricea C este parcursa pe linii, bucla ın k calculand elementul cij . Similar, ın

forma jik, matricea C este parcursa pe coloane. In celelalte ordini, elementelelui C nu sunt calculate succesiv, ci prin calcule ıntretesute, dar suma finalapentru fiecare este aceeasi.

Forme ale produsului matriceal. Vom evidentia acum mai multe formede prezentare a ınmultirii a doua matrice, folosind ca operatii elementarecazurile particulare prezentate mai sus. Pentru aceasta, este util sa privimo matrice ca o concatenare de vectori. De exemplu, evidentiind coloanele, omatrice poate fi vazuta ca

A = [a1 a2 . . . an], (1.12)

cu a1, a2, . . . , an ∈ Rm notandu-se coloanele matricei A. Evidentiind liniile,o matrice poate fi scrisa

A =

aT1aT2...aTm

. (1.13)

Atentie: acum aTi e o notatie pentru vectorul format de linia i a matricei A;ai nu este coloana i din (1.12).

1. Privim A ca alaturare de vectori linie, ca ın (1.13), si B ca alaturare devectori coloana, ca ın (1.12). Se observa imediat din definitie ca cij = aTi bj ,adica fiecare element al rezultatului poate fi exprimat printr-un produs scalar,asa cum se vede si din Figura 1.6.

2. Privim acum A partitionata pe coloane si B pe linii. Produsul lor sepoate exprima ca o suma de produse exterioare:

C = AB = [a1 . . . ap]

bT1...bTp

=

p∑k=1

akbTk . (1.14)

Demonstratie: cij =∑pk=1(akb

Tk )ij =

∑pk=1(ak)i(b

Tk )j =

∑pk=1 aikbkj .

Ilustratie grafica: Figura 1.7.

3. Punem acum ın evidenta numai coloanele matricei B. Atunci

1.6 Operatii elementare cu vectori si matrice 19

A B

· =

C

Fig. 1.6 Produsul matriceal exprimat prin produse scalare.

A B

· = · + · + ·

Fig. 1.7 Produs matriceal exprimat prin produse exterioare.

C = AB = A[b1 . . . bn] = [Ab1 . . . Abn], (1.15)

deci fiecare coloana a produsului este obtinuta prin ınmultirea matrice-vectordintre matricea A si coloana respectiva a lui B.

4. Fie acum A partitionata pe linii. Atunci

C = AB =

aT1...aTm

B =

aT1 B...

aTmB

, (1.16)

deci fiecare linie a produsului este obtinuta prin ınmultirea vector linie -matrice dintre linia respectiva a lui A si matricea B.

5. Exprimam acum coloanele produsului C ın functie de coloanele matriceiA. Cu ajutorul relatiei (1.15) se obtine

cj = Abj =

p∑k=1

bkjak, (1.17)

deci orice coloana din C este combinatie liniara a coloanelor matricei A.

6. In sfarsit, o ultima forma, ın care exprimam liniile produsului C ınfunctie de liniile matricei B. Din (1.16) rezulta

cTi = aTi B =

p∑k=1

aikbTk . (1.18)


Deci, o linie din C este combinatie liniara a liniilor matricei B.

Capitolul 2

Sisteme de ecuatii liniare

Un sistem de m ecuatii liniare cu n necunoscute are formaa11x1 + a12x2 + . . .+ a1nxn = b1a21x1 + a22x2 + . . .+ a2nxn = b2. . .am1x1 + am2x2 + . . .+ amnxn = bm

, (2.1)

unde aij ∈ R, bi ∈ R, pentru i = 1 : m, j = 1 : n, sunt date si xj , j = 1 : n,

sunt necunoscutele. In forma matriceala, sistemul este

Ax = b, (2.2)

A ∈ Rm×n fiind numita matricea coeficientilor, b ∈ Rm termenul liber (saupartea dreapta a sistemului), iar x ∈ Rn vectorul necunoscutelor.

O caracteristica importanta a unui sistem liniar este relatia dintre m si n:

• ın cazul unei matrice A patrate (m = n), sistemul (2.2) este determinat;

• cand numarul de ecuatii este mai mare decat numarul de necunoscute(m > n), sistemul (2.2) se numeste supradeterminat;

• cand, dimpotriva, numarul de necunoscute este mai mare (n > m), sis-temul este numit subdeterminat.

In acest capitol ne vom ocupa exclusiv de sisteme determinate. Celelaltedoua tipuri de sisteme, supra- sau subdeterminate, vor fi discutate ın Cap. 3.In aceste capitole va fi considerat doar cazul ın care matricea A are rangmaxim, anume min(m,n). Cazul ın care rangul nu e maxim va fi tratat ınCap. 5.

Presupunem deci, ın restul acestui capitol, ca m = n si rangA = n. Inacest caz, sistemul (2.2) are solutie unica, anume

x = A−1b. (2.3)

21

22 2 Sisteme de ecuatii liniare

Subliniem ca (2.3) nu este o formula adecvata calculului numeric al solutieisistemului Ax = b; calculul inversei este ineficient si are potential mai mare deinstabilitate decat alte metode. Vom vedea ın acest capitol metode numericeeficiente si precise pentru a face acest lucru.

In practica numerica actuala se folosesc doua categorii principale demetode numerice pentru calcularea solutiei unui sistem determinat:

• Metode directe, bazate pe reducerea sistemului (2.2), printr-o secventafinita de transformari elementare, la unul sau doua sisteme triunghiu-lare, care se rezolva simplu utilizand proceduri de substitutie descrise ınsectiunea 2.1. Din aceasta categorie fac parte metodele de eliminare gaus-siana si de factorizare LU , care sunt recomandate pentru sisteme de di-mensiuni mici si medii; orientativ, pentru n = 1000 rezolvarea unui sistemdureaza cam o secunda pe un laptop obisnuit, deci dimensiuni de ordinulmiilor pot cu siguranta fi considerate ca nefiind mari. Metodele directesunt subiectul principal al acestui capitol.

• Metode iterative, bazate pe construirea recursiva a unui sir de vectori careconverge catre solutia sistemului (2.2). Aceste metode sunt recomandatepentru sisteme de dimensiuni foarte mari, sau/si atunci cand matricea Aare o structura speciala, eventual rara. Metodele iterative vor fi ilustrateın sectiunea 2.11.

2.1 Rezolvarea sistemelor triunghiulare

Definitia 2.1 O matrice T se numeste inferior triunghiulara daca tij = 0,pentru i < j, si superior triunghiulara daca tij = 0, pentru i > j. Figura 2.1prezinta astfel de matrice.

O matrice se numeste triunghiulara unitate daca este triunghiulara si, ınplus, toate elementele sale diagonale sunt egale cu 1.

O matrice se numeste strict triunghiulara daca este triunghiulara si toateelementele sale diagonale sunt nule. De exemplu, o matrice A este strictinferior triunghiulara daca aij = 0, pentru i ≤ j.

Sistemul Ax = b este numit inferior (superior) triunghiular daca matriceaA este inferior (superior) triunghiulara.

O matrice triunghiulara este inversabila daca si numai daca toate ele-mentele sale diagonale sunt nenule.

Algoritmii pentru rezolvarea sistemelor triunghiulare sunt foarte sim-pli, deoarece necunoscutele pot fi calculate, ıntr-o ordine precizata, prinsubstitutie numerica.

Sa consideram ıntai un sistem inferior triunghiular, Lx = b, unde L ∈Rn×n, cu lij = 0 pentru i < j, lii 6= 0, i ∈ 1 : n, si b ∈ Rn. Prima ecuatie este

2.1 Rezolvarea sistemelor triunghiulare 23

Fig. 2.1 Matrice inferior (stanga) si superior (dreapta) triunghiulare

l11x1 = b1, din carex1 = b1 / l11. (2.4)

In general, daca se cunosc x1, x2, . . . , xi−1, se poate rezolva a i-a ecuatie,

i−1∑j=1

lijxj + liixi = bi,

pentru a obtine

xi =

bi − i−1∑j=1

lijxj

/ lii. (2.5)

Formulele (2.4) si (2.5) definesc un algoritm pentru calculul solutiei sis-temului Lx = b, cunoscut ca metoda substitutiei ınainte.

Algoritmul 2.1 (LTRIS – Calculul solutiei unui sistem inferior tri-unghiular) (Se dau L ∈ Rn×n, inferior triunghiulara, inversabila, sib ∈ Rn. Se calculeaza solutia x a sistemului Lx = b.)

1. x← b2. Pentru i = 1 : n

1. Pentru j = 1 : i− 11. xi ← xi − lijxj

2. xi ← xi / lii

Comentarii. Fiecare executie a buclei Pentru i necesita 2(i − 1) flopi ın2.1.1, pentru calculul sumei din (2.5), si 1 flop pentru ımpartirea din 2.2. Intotal, sunt necesari NLTRIS =

∑ni=1 (2(i− 1) + 1) = n(n− 1) + n = n2 flopi

si MLTRIS = n(n+1)/2+2n+1 ≈ n2/2 locatii de memorie (pentru elementeın virgula mobila). Daca vectorul b nu este necesar ulterior, solutia x poate fimemorata ın aceeasi zona de memorie cu b; se spune ca se efectueaza calcululpe loc ın b. Ulterior, vom apela acest algoritm prin x = LTRIS(L, b).

Un algoritm cu aceleasi proprietati se obtine daca schimbam ordinea celordoua bucle, ceea ce revine la a parcurge matricea L pe coloane, nu pe linii.Substitutia decurge acum astfel: dupa calculul lui x1 cu (2.4), se actualizeaza


valorile celorlalte necunoscute xi, i = 2 : n (initializate cu bi), cu termenulli1x1 care contribuie la sumele din (2.5); se poate calcula acum x2, dupa carese actualizeaza din nou restul necunoscutelor, etc.

Algoritmul 2.2 (LTRIS – versiunea pe coloane) (Se dau L ∈ Rn×n,inferior triunghiulara, inversabila, si b ∈ Rn. Se calculeaza solutia x asistemului Lx = b.)

1. x← b2. Pentru j = 1 : n

1. xj ← xj / ljj2. Pentru i = j + 1 : n

1. xi ← xi − lijxj

Desigur, numarul operatiilor este acelasi ın cele doua versiuni ale algorit-mului.

Pentru a rezolva sistemul superior triunghiular (nesingular) Ux = b, undeU ∈ Rn×n este o matrice superior triunghiulara, cu uii 6= 0, i ∈ 1 : n, sib ∈ Rn, observam ca ultima ecuatie are forma unnxn = bn, iar de aici

xn = bn/unn. (2.6)

Mai mult, daca xn, xn−1, . . . , xi+1, sunt deja calculate, atunci, din a i-aecuatie,

uiixi +

n∑j=i+1

uijxj = bi,

se obtine

xi =

bi − n∑j=i+1

uijxj

/uii. (2.7)

Din formulele (2.6) si (2.7), care definesc asa numita substitutie ınapoi,obtinem direct algoritmul urmator.

Algoritmul 2.3 (UTRIS – Calculul solutiei unui sistem superior tri-unghiular) (Se dau U ∈ Rn×n, superior triunghiulara, inversabila, sib ∈ Rn. Se calculeaza solutia x a sistemului Ux = b.)

1. x← b2. Pentru i = n : −1 : 1

1. Pentru j = i+ 1 : n1. xi ← xi − uijxj

2. xi ← xi/uii

Perfect analog cu algoritmul 2.2 avem versiunea pe coloane a UTRIS.

2.2 Transformari elementare 25

Algoritmul 2.4 (UTRIS – versiunea pe coloane) (Se dauU ∈ Rn×n, superior triunghiulara, inversabila, si b ∈ Rn. Se calculeazasolutia x a sistemului Ux = b.)

1. x← b2. Pentru j = n : −1 : 1

1. xj ← xj/ujj1. Pentru i = 1 : j − 1

1. xi ← xi − uijxj

Comentarii. Ca si ın cazul inferior triunghiular, efortul de calcul esteNUTRIS = n2 flopi, iar memoria ocupata MUTRIS ≈ n2/2 locatii ın virgulamobila. Vom apela acest algoritm prin x = UTRIS(L, b).

2.2 Transformari elementare

Transformarile matricei A necesare pentru rezolvarea sistemului (2.2) prinreducere la forma superior triunghiulara se pot descrie concis si elegant ıntermeni de operatii matriceale, mai precis ca produs ıntre A si anumite ma-trice elementare. (In loc de produs sau ınmultire, vom vorbi despre operatieca fiind aplicarea asupra lui A a matricei de transformare corespunzatoare.)

Definitia 2.2 O matrice inferior triunghiulara elementara (ITE) de ordin nsi indice k este o matrice de forma

Mk = In −mkeTk , (2.8)

undemk = [0 0 . . . 0 µk+1,k . . . µnk ]T (2.9)

este un vector avand primele k componente nule, iar ek este vectorul unitatede indice k.

Din motive ce vor fi evidente ulterior, Mk se mai numeste transformareelementara sau gaussiana, mk se numeste vector Gauss, iar scalarii µik senumesc multiplicatori (gaussieni).

Structura unei matrice ITE este

Mk =

1 0 . . . 0 . . . 00 1 . . . 0 . . . 0

. . . . . .0 0 . . . 1 . . . 00 0 . . . −µk+1,k . . . 0

. . . . . .. . . 0

0 0 . . . −µnk . . . 1

.


Notam de acum faptul ca o matrice ITE este complet definita de vectorulmk din (2.9); ın orice algoritm, matricea Mk nu se memoreaza integral, ci doarprin vectorul mk, adica prin multiplicatorii µik, i = k + 1 : n. Principaleleproprietati ale unei astfel de matrice sunt rezumate ın urmatoarea propozitie.

Propozitia 2.3 a) O matrice ITE este inversabila si

M−1k = In +mkeTk . (2.10)

b) Fie x ∈ Rn un vector dat, iar Mk o matrice ITE precizata. Atuncivectorul transformat y = Mkx are elementele

(Mkx)i =

{xi, pentru i = 1 : k,

xi − µikxk, pentru i = k + 1 : n.(2.11)

c) Fie x ∈ Rn. Daca xk 6= 0, atunci exista o matrice ITE Mk astfel ıncatvectorul transformat y = Mkx are ultimele n−k componente nule, mai precis

(Mkx)i =

{xi, pentru i = 1 : k,

0, pentru i = k + 1 : n.(2.12)

Daca xk = 0, atunci pentru orice Mk

Mkx = x. (2.13)

Demonstratie. a) Evident, det(Mk) = 1, adica Mk este inversabila. Maimult,

Mk(In +mkeTk ) = (In −mke

Tk )(In +mke

Tk ) = In −mk(eTkmk)eTk = In

deoarece, din (2.9), eTkmk = 0. Deci (2.10) este adevarata.

b) Utilizand (2.8), obtinem

Mkx = (In −mkeTk )x = x−mke

Tk x = xi − xkmk.

Scriind relatia de mai sus pe componente si tinand seama de (2.9), obtinem(2.11).

c) Tinand seama de (2.11) si alegand

µik = xi/xk, i = k + 1 : n, (2.14)

rezulta (Mkx)i = 0, pentru i = k+ 1 : n. Daca xk = 0, atunci (2.11) conducedirect la (2.13). De observat ca, daca xk 6= 0, atunci matricea ITE definitade vectorul Gauss cu elementele din (2.14) este unica matrice ITE de indicek care introduce zerouri ın ultimele n− k componente ale vectorului x.

2.2 Transformari elementare 27

Dupa cum vom vedea ın sectiunea urmatoare, proprietatile (2.12) si (2.13)sunt cruciale ın reducerea unei matrice la forma triunghiulara.

Vom introduce acum notiunea de matrice de permutare, care permite de-scrierea concisa a interschimbarilor de linii sau coloane ın termeni de operatiimatriceale.

Definitia 2.4 O matrice Pij ∈ Rn×n obtinuta din matricea unitate In prininterschimbarea a doua coloane (sau linii) i si j, adica o matrice de forma(aici i < j):

Pij = [e1 e2 . . . ei−1 ej ei+1 . . . ej−1 ei ej+1 . . . en]

=

Ii−1

0 1Ij−i−1

1 0In−j

ın care toate elementele nefigurate sunt nule, se numeste (matrice de) per-mutare elementara (PE).

De exemplu, permutarea P24 ∈ R4×4 are forma

P24 =

1 0 0 00 0 0 10 0 1 00 1 0 0

.Daca (i1, i2, . . . , in) este o permutare a multimii 1 : n, atunci matricea

P = [ei1 ei2 . . . ein ] este o matrice de permutare. Desigur, un produs dematrice elementare de permutare este o matrice de permutare si, reciproc,orice matrice de permutare se poate exprima ca produs de matrice PE.

Proprietatile matricelor elementare de permutare sunt expuse ın urmatoareapropozitie, a carei demonstratie e evidenta.

Propozitia 2.5 a) O PE este ortogonala si simetrica, deci P−1ij = Pij, iardetPij = −1.

b) Premultiplicarea unei matrice A cu PE Pij interschimba linia i cu liniaj din A, adica:

eTk (PijA) =

eTkA, pentru k 6= i, j,

eTj A, pentru k = i,

eTi A, pentru k = j.

c) Postmultiplicarea unei matrice A cu PE Pij interschimba coloana i cucoloana j din A, adica:


i

j

Pij ·

i j

·Pij

Fig. 2.2 Inmultirea cu o matrice de permutare elementara, la stanga (figura de sus) si la

dreapta (jos).

(APij)ek =

Aek, pentru k 6= i, j,

Aej , pentru k = i,

Aei, pentru k = j.

Aceste proprietati sunt ilustrate ın Figura 2.2.

2.3 Triangularizare prin eliminare gaussiana

Fie A ∈ Rn×n (cazul A ∈ Cn×n este identic) o matrice nu neaparat in-versabila. Eliminarea gaussiana este o metoda de reducere a matricei A laforma superior triunghiulara prin aplicarea la stanga lui A a unei secventeMk, k = 1 : n − 1, de matrice ITE, fiecare aleasa astfel ıncat sa anulezeelementele subdiagonale ın coloana corespunzatoare ak a matricei A (fara aaltera zerourile introduse anterior).

Pentru prezentarea metodei avem nevoie de doua rezultate cu caractertehnic.

Propozitia 2.6 Fie A ∈ Rn×n. Notam A[k] def= A(1 : k, 1 : k) submatricea

lider principala (sau blocul stanga sus) de ordin k a matricei A.

a) Daca L ∈ Rn×n este inferior triunghiulara, atunci

(LA)[k] = L[k]A[k].

b) Daca Li ∈ Rn×n, i = 1 : p, sunt matrice inferior triunghiulare, atunci

(L1L2 . . . LpA)[k] = L[k]1 L

[k]2 . . . L[k]

p A[k].

2.3 Triangularizare prin eliminare gaussiana 29

Demonstratie. a) Evidentiem blocurile din fiecare matrice (cel din stangasus avand dimensiune k × k iar blocurile care nu ne intereseaza direct fiindnotate cu ∗) si obtinem

LA =

[L[k] 0∗ ∗

] [A[k] ∗∗ ∗

]=

[L[k]A[k] ∗∗ ∗

].

b) Prin aplicarea succesiva a punctului a).

Urmatoarea teorema stabileste conditiile ın care este posibila reducereaunei matrice la forma triunghiulara, utilizand transformari ITE.

Teorema 2.7 Daca matricea A ∈ Rn×n satisface conditia

(i) submatricele lider principale A[k], k = 1 : n− 1, sunt nesingulare,

atunci exista o matricele ITE Mk, k = 1 : n− 1, astfel ıncat matricea

Mn−1 . . .M2M1A = U

este superior triunghiulara. Daca A este nesingulara, atunci U este nesingu-lara.

Demonstratia este constructiva, adica descrie o procedura pentru reducereaefectiva a matricei A la forma superior triunghiulara U . Procedura consta ınn− 1 pasi.

Pasul 1. Fie A1 = A si a1 = Ae1 prima sa coloana. Prin ipoteza, A[1]1 =

A[1] = a11 6= 0. Din Propozitia 2.3, exista o matrice ITE M1 astfel ıncat

(M1a1)i =

{a11, pentru i = 1,

0, pentru i = 2 : n,

adica matricea A2 = M1A1 are toate elementele subdiagonale ale primeicoloane egale cu zero:

A2 =

a(2)11 a

(2)12 . . . a

(2)1n

0 a(2)22 . . . a

(2)2n

0 a(2)32 . . . a

(2)3n

. . . . . . . . . . . .

0 a(2)n2 . . . a

(2)nn

.

Pasul k. Presupunem ca matricea curenta

Ak = Mk−1 . . .M2M1A

este superior triunghiulara ın primele k − 1 coloane, adica are toate ele-mentele subdiagonale ale acestor coloane egale cu zero. Aplicand Propozitia


2.6, obtinem

A[k]k = M

[k]k−1 . . .M

[k]1 A[k],

unde matriceleM[k]i sunt toate inferior triunghiulare unitate, deci det(M

[k]i ) =

1, i = 1 : k−1. Deoarece A[k]k este superior triunghiulara, ın virtutea conditiei

(i) din enuntul teoremei avem

det(A[k]k ) =

k∏i=1

a(k)ii = det(A[k]) 6= 0,

adica elementul a(k)kk , numit pivot, este nenul. Consideram partitionarea pe

coloane a matricei Ak:

Ak = [ a(k)1 a

(k)2 . . . a

(k)k . . . a(k)n ].

Putem utiliza acum Propozitia 2.3 pentru a conchide ca exista o matrice ITE

Mk astfel ıncat (Mka(k)k )i = 0, pentru i = k+1 : n. Mai mult, premultiplicarea

cu Mk a matricei Ak nu altereaza primele k − 1 coloane si ın particular

zerourile deja introduse la pasii anteriori, deoarece a(k)kj = 0, pentru j = 1 :

k − 1 (vezi Propozitia 2.3, punctul c). Prin urmare matricea

Ak+1 = MkAk (2.15)

este superior triunghiulara ın primele k coloane.

Astfel, procedura de introducere a zerourilor subdiagonale, ınceputa lapasul 1 si continuata pana la pasul n− 1 (inclusiv), are ca rezultat matricea

Udef= An = Mn−1Mn−2 . . .M1A, (2.16)

evident superior triunghiulara. Mai mult, deoarece produsul matriceal con-serva structura inferior triunghiulara unitate (demonstrati!), matricea

M = Mn−1Mn−2 . . .M1

este inferior triunghiulara unitate iar, daca A este nesingulara, atunci si Urezulta nesingulara, ca produs a doua matrice nesingulare.

Demonstratia de mai sus furnizeaza o procedura de triangularizare a uneimatrice, cunoscuta sub numele de eliminare gaussiana. Procedura este bazatape urmatoarea schema, ın care calculele se desfasoara pe loc ın tabloul A

ElG 1. Pentru k = 1 : n− 11. Calculeaza matricea ITE Mk (adica multiplicatorii µik, i =k + 1 : n), astfel ıncat (MkA)i = 0, pentru i = k + 1 : n

2. Calculeaza A←MkA

2.3 Triangularizare prin eliminare gaussiana 31

u11 u12 . . . u1k u1,k+1 . . . u1nµ21 u22 . . . u2k u2,k+1 . . . u2n

. . . . . .

µk1 µk2 . . . ukk uk,k+1 . . . ukn

µk+1,1 µk+1,2 . . . µk+1,k a(k+1)k+1,k+1 . . . a

(k+1)k+1,n

. . . . . .

µn1 µn2 . . . µnk a(k+1)n,k+1 . . . a

(k+1)nn

u11 u12 . . . u1k . . . u1nµ21 u22 . . . u2k . . . u2n

. . . . . .

µk1 µk2 . . . ukk . . . ukn. . . . . .. . . . . .

µn1 µn2 . . . µnk . . . unn

Dupa pasul k In final

Fig. 2.3 Continutul matricei A dupa pasul k al eliminarii gaussiene si la terminareaexecutiei algoritmului

La pasul 1.1, multiplicatorii gaussieni µik sunt calculati conform (2.14)prin

µik = a(k)ik /a

(k)kk , i = k + 1 : n, (2.17)

si pot fi memorati pe pozitiile elementelor anulate. Prin instructiunea 1.2 aschemei de mai sus, toate rezultatele intermediare sunt depuse ın spatiul dememorie ocupat de matricea A; vezi si Figura 2.3.

La pasul 1.2, transformarea A ← MkA este echivalenta cu calcululaj ← Mkaj , pentru toate coloanele aj = Aej ale matricei A. Transformareapoate fi efectuata eficient tinand cont ca primele k−1 coloane ale matricei Anu sunt afectate. Coloana k este modificata prin anularea elementelor subdi-agonale. Pentru restul submatricei, tinand seama de Propozitia 2.3b, coloanatransformata are elementele

aij ← (Mkaj)i = ((In −mkeTk )aj)i = aij − µikakj , i = k + 1 : n, (2.18)

unde j = k + 1 : n. Relatiile (2.17) si (2.18) ımpreuna cu schema de

calcul ElG definesc algoritmul de eliminare gaussiana detaliat ın contin-

uare. Conditia de terminare corecta a algoritmului este ca numerele a(k)kk ,

k = 1 : n − 1, numite elemente pivot, sa fie nenule, adica submatricele liderprincipale A[k], k = 1 : n− 1, ale matricei initiale, sa fie nesingulare.

Algoritmul 2.5 (G – Eliminare gaussiana) (Se da A ∈ Rn×n, cu sub-matricele lider principale A[k], k = 1 : n− 1, nesingulare. Se calculeazamatricea superior triunghiulara U si matricele ITE M1, . . . , Mn−1,astfel ıncat U = Mn−1Mn−2 . . .M1A. Matricea U se memoreaza ın tri-unghiul superior al lui A, inclusiv diagonala, iar multiplicatorii gaussieniµik care definesc transformarile Mk se memoreaza ın triunghiul inferioral matricei A, pe pozitiile elementelor anulate.)

1. Pentru k = 1 : n− 11. Pentru i = k + 1 : n


1. aik ← µik = aik/akk2. Pentru j = k + 1 : n

1. Pentru i = k + 1 : n1. aij ← aij − µikakj

Comentarii. Algoritmul G necesita

NG =

n−1∑k=1

(n− k + 2(n− k)2) =n(n− 1)

2+ 2

n(n− 1)(2n− 1)

6≈ 2n3

3,

operatii ın virgula mobila (flopi), iar memoria ocupata este de MG = n2

locatii (elemente ın virgula mobila).

Algoritmul G poate fi utilizat pentru rezolvarea sistemului liniar Ax = batunci cand A are toate submatricele lider principale nesingulare (inclusivA[n] = A). Intr-adevar, sistemul Ax = b are aceeasi solutie ca sistemul supe-rior triunghiular

Ux = Mn−1Mn−2 . . .M1b,

ın care matricea U este furnizata de algoritmul G. Acest sistem triunghiularpoate fi rezolvat cu algoritmul UTRIS.

Semnificatia operatiilor din algoritmul G este imediata; introducerea dezerouri subdiagonale ın coloana k ınseamna eliminarea necunoscutei xk dinecuatiile k + 1 : n; buclele 1.2 si 1.2.1 reprezinta scaderea din ecuatiile i =k + 1 : n a ecuatiei k ınmultite, pe rand, cu multiplicatorii µik, cu scopul dea elimina xk din fiecare ecuatie i.

Totusi, nesingularitatea submatricelor lider principale ale matricei A nueste o conditie necesara pentru existenta si unicitatea solutiei sistemului (2.8).De aceea, algoritmul G trebuie modificat astfel ıncat sa produca rezultatuldorit si atunci cand unele submatrice lider principale sunt singulare, dar Aeste nesingulara.

2.4 Strategii de pivotare

Daca o submatrice lider principala A[k] este singulara, atunci, la pasul k

al procedurii de eliminare gaussiana, pivotul a(k)kk este nul; deci, conform

Propozitiei 2.3c, nu exista nici o matrice ITE Mk care sa anuleze, prin pre-multiplicare, elementele subdiagonale ale coloanei k din Ak. Modificarea al-goritmului G recomandata pentru a evita aceasta dificultate consta ıntr-ointerschimbare de linii (sau/si coloane), care aduce ın pozitia pivotului unelement nenul. Un alt motiv pentru permutarea de linii si coloane este asig-urarea stabilitatii numerice a procesului de calcul. O explicatie intuitiva esteca, atunci cand un proces matematic nu poate fi definit pentru o anumita

2.4 Strategii de pivotare 33

Ak =

u11 . . . u1k . . . u1n

0. . . . . .

a(k)kk . . . a

(k)kn

0 . . .

a(k)ikk

. . . a(k)ikn

0 . . .

a(k)nk . . . a

(k)nn

PkikAk =

u11 . . . u1k . . . u1n

0. . . . . .

a(k)ikk

. . . a(k)ikn

0 . . .

a(k)kk . . . a

(k)kn

0 . . .

a(k)nk . . . a

(k)nn

Fig. 2.4 Pasul k al eliminarii gaussiene cu interschimbare de linii; a(k)kk este nul sau prea

aproape de zero; ın PkikAk, noul pivot este a(k)ikk

si are o valoare ”buna” (suficient de mare,

ın modul)

valoare ”critica” a unui parametru (de exemplu A[k] singulara), exista sansemari ca procesul numeric corespunzator sa fie instabil atunci cand parametrulefectiv este ın vecinatatea valorii critice. In cazul nostru, trebuie evitati pivotiimici (ın modul), corespunzatori cazului ın care submatricele lider principalesunt ”aproape” singulare; astfel se ımbunatateste stabilitatea numerica a al-goritmului.

2.4.1 Pivotare partiala

Vom modifica algoritmul G dupa cum urmeaza. La pasul k al proceduriide triangularizare a matricei A ∈ Rn×n (vezi demonstratia Teoremei 2.7),se alege ca pivot elementul maxim ın modul din coloana k, de la diagonala(inclusiv) ın jos; daca sunt mai multe elemente maxime egale, se poate alegeoricare dintre ele, de preferat primul pentru ca astfel el se poate afla chiarpe diagonala. Presupunand ca acest element se gaseste ın linia ik si este

nenul (a(k)ikk6= 0), interschimbarea liniilor k si ik, obtinuta matriceal prin

ınmultirea la stanga a lui Ak cu PE Pkik , muta a(k)ikk

ın pozitia (k, k), astfelıncat acest element devine noul pivot. Pentru a evita alterarea structuriimatricei obtinute pana la pasul k (superior triunghiulara ın primele k − 1coloane), trebuie ıntr-adevar sa alegem ik ≥ k. Vezi Figura 2.4.

Pentru simplitate, vom nota matricea Pikk cu Pk. Pasul k al algoritmuluiG devine

1. Se determina ik ≥ k astfel ıncat |aikk| = maxi=k:n |aik|.2. Se interschimba liniile ik si k, adica A← PkA.

3. Se determina matricea ITE Mk astfel ıncat (MkA)ik = 0, i = k + 1 : n.

4. Se aplica transformarea A←MkA.


Astfel la pasul k are loc transformarea

A← Ak+1 = MkPkAk. (2.19)

Procesul global, cunoscut sub numele de eliminare gaussiana cu pivotarepartiala, conduce ın cele din urma la matricea superior triunghiulara

Udef= An = Mn−1Pn−1Mn−2Pn−2 . . .M1P1A. (2.20)

Algoritmul corespunzator este prezentat ın continuare.

Algoritmul 2.6 (GPP – Eliminare gaussiana cu pivotare partiala)(Se da A ∈ Rn×n, nesingulara. Se calculeaza matricea superior triun-ghiulara U , matricele ITE M1, . . . , Mn−1 si PE P1, . . . , Pn−1 astfelıncat (2.20) sa fie satisfacuta. Matricea U se memoreaza ın triunghiulsuperior al lui A, inclusiv diagonala, iar multiplicatorii gaussieni µikcare definesc matricele Mk se memoreaza ın triunghiul inferior al luiA. Intregii ik care definesc PE Pk, sunt memorati ıntr-un vector p dedimensiune n− 1, cu elementele p(k) = ik, pentru k = 1 : n− 1.)

1. Pentru k = 1 : n− 11. Se determina ik ∈ k : n astfel ıncat |aikk| =

maxi=k:n |aik|.2. p(k)← ik3. Pentru j = k : n % interschimba liniile k si ik

1. akj ↔ aikj4. Pentru i = k + 1 : n



Comentarii. Operatiile suplimentare necesare ın algoritmul GPP, fata dealgoritmul G, apar numai ın cautarea pivotului si constau ın comparatiilenecesare gasirii maximului din instructiunea 1.1; numarul de comparatii este

n−1∑k=1

(n− k + 1) ≈ n2/2� 2n3/3,

deci pivotarea partiala nu este costisitoare (relativ la numarul de operatiiaritmetice). Necesarul de memorie este MGPP = n2.

Un rezultat important referitor la algoritmul GPP este dat de

Teorema 2.8 Daca A ∈ Rn×n este nesingulara, atunci algoritmul GPP seva executa complet iar matricea superior triunghiulara U rezultata ın final,vezi (2.20), este nesingulara.

2.4 Strategii de pivotare 35

Ak =

u11 . . . u1k . . . u1jk . . . u1n

0. . . . . . . . .

a(k)kk . . . a

(k)kjk

. . . a(k)kn

0 . . .

a(k)ikk

. . . a(k)ikjk

. . . a(k)ikn

0 . . .

a(k)nk . . . a

(k)njk

. . . a(k)nn

PkikAkQkjk =

u11 . . . u1jk . . . u1k . . . u1n

0. . . . . . . . .

a(k)ikjk

. . . a(k)ikk

. . . a(k)ikn

0 . . .

a(k)kjk

. . . a(k)kk . . . a

(k)kn

0 . . .

a(k)njk

. . . a(k)nk . . . a

(k)nn

Fig. 2.5 Pasul k al eliminarii gaussiene, cu interschimbari de linii si coloane

Demonstratie. Algoritmul GPP nu se executa complet daca, la un pask, elementul pivot aikk este nul. Intr-un astfel de caz, deoarece pivotul estemaxim ın modul pe coloana (sub diagonala), toate elementele aik, i = k : n,sunt nule, deci matricea Ak este singulara (vezi structura matricei Ak, cu

A[k]k superior triunghiulara).

Matricele ITE si PE fiind nesingulare, singularitatea matricei

Ak = Mk−1Pk−1 . . .M1P1A

implica A singulara; contradictie. Deci, toti pivotii aikk sunt nenuli si algo-ritmul se executa complet. Nesingularitatea matricei U rezulta din nesingu-laritatea tuturor matricelor din termenul drept din (2.20).

2.4.2 Pivotare completa

Proprietati numerice superioare se pot obtine daca, la pasul k al eliminarii

gaussiene, pivotul maxim ın modul este cautat printre elementele a(k)ij , i =

k : n, j = k : n, adica ın ıntregul colt dreapta jos al matricei Ak. Unelement cu modul maxim (de pilda primul ın ordinea explorarii pe coloane),sa zicem aikjk , este mutat ın pozitia pivot (k, k), prin intermediul a douainterschimbari, una de linii, cealalta de coloane, exprimate matriceal prinınmultirea lui Ak la stanga, respectiv la dreapta, cu PE Pkik , respectiv Qkjk ;vezi Figura 2.5.

Notand Pkik cu Pk si Qkjk cu Qk, pasul k al algoritmului G devine:

1. Se determina cei mai mici ik si jk astfel ıncat |aikjk | = maxi=k:n,j=k:n |aij |.2. Se interschimba liniile ik si k, adica A← PkA.

3. Se interschimba coloanele jk si k, adica A← AQk.

4. Se determina matricea ITE Mk astfel ıncat (MkA)ik = 0, i = k + 1 : n.


5. Se aplica transformarea A←MkA.

Astfel, la pasul k are loc transformarea

A← Ak+1 = MkPkAkQk, (2.21)

iar procesul global de calcul, cunoscut sub numele de eliminare gaussiena cupivotare completa, conduce la matricea superior triunghiulara

A← Udef= Mn−1Pn−1 . . .M1P1AQ1 . . . Qn−1. (2.22)

Algoritmul rezultat este prezentat ın continuare.

Algoritmul 2.7 (GPC – Eliminare gaussiana cu pivotare completa)(Se da A ∈ Rn×n, nesingulara. Se calculeaza matricea superior triun-ghiulara U , matricele ITE M1, . . . , Mn−1 si PE P1, Q1, . . . , Pn−1,Qn−1, astfel ıncat (2.22) sa fie satisfacuta. Matricea U se memoreazaın triunghiul superior al lui A, inclusiv diagonala, iar multiplicatoriigaussieni µik care definesc matricele Mk se memoreaza ın triunghiulinferior al matricei A. Intregii ik si jk, definind PE Pk, respectiv Qk,sunt memorati ın vectorii p, respectiv q, de dimensiune n − 1, astfelıncat p(k) = ik si q(k) = jk, pentru k = 1 : n− 1.)

1. Pentru k = 1 : n− 11. Se determina ik ∈ k : n si jk ∈ k : n astfel ıncat|aikjk | = maxi=k:n,j=k:n |aij |.

2. p(k)← ik3. q(k)← jk4. Pentru j = k : n % interschimba liniile k si ik

1. akj ↔ aikj5. Pentru i = 1 : n % interschimba coloanele k si jk

1. aik ↔ aijk6. Pentru i = k + 1 : n



Comentarii. Operatiile suplimentare ce apar ın algoritmul GPC, fata dealgoritmul G, sunt doar cele necesare cautarii pivotului, ca si la pivotareapartiala; acum ınsa, numarul de comparatii este

n−1∑k=1

(n− k + 1)2 =

n∑k=2

k2 ≈ n3/3;

spre deosebire de pivotarea partiala, pivotarea completa introduce un numarde comparatii de acelasi ordin de marime cu numarul operatiilor aritmetice;

2.5 Factorizari LU 37

aceasta poate conduce la cresterea semnificativa a timpului de executie, deaceea pivotarea completa trebuie folosita doar atunci cand se doreste o maresiguranta ın calculul rezultatului. Vom prezenta ulterior argumentele careconduc la alegerea strategiei de pivotare. Necesarul de memorie este si aiciMGPC = n2.

Teorema urmatoare corespunde Teoremei 2.8, demonstratia fiind asemana-toare.

Teorema 2.9 Daca A ∈ Rn×n este nesingulara, atunci algoritmul GPC seexecuta complet, iar matricea superior triunghiulara rezultata U din (2.22)este nesingulara.

2.5 Factorizari LU

In multe situatii, este convenabil a exprima o matrice data ca produs a douamatrice triunghiulare. Aceasta forma poate fi utilizata pentru a simplificaanumite calcule implicand matricea respectiva.

Definitia 2.10 Fie A ∈ Rn×n. Daca exista o matrice inferior triunghiularaL ∈ Rn×n si o matrice superior triunghiulara U ∈ Rn×n astfel ıncat

A = LU (2.23)

atunci se spune ca A admite o factorizare (sau descompunere) LU. L si Use numesc factori triunghiulari ai matricei A.

Daca A ∈ Rn×n este o matrice care are o factorizare LU, iar D esteo matrice diagonala nesingulara, atunci, D−1 fiind de asemenea diagonala,(2.23) poate fi scrisa ın forma

A = LU = LDD−1U = L′U ′, (2.24)

unde L′ = LD este inferior triunghiulara si U ′ = D−1U este superior tri-unghiulara. Daca D 6= I, atunci factorizarile (2.23) si (2.24) difera si decifactorizarea LU, daca exista, nu este unica. De aceea este util a se introducedefinitia urmatoare.

Definitia 2.11 Fie A ∈ Rn×n. Daca exista matricele L ∈ Rn×n inferiortriunghiulara unitate, U ∈ Rn×n superior triunghiulara unitate si D ∈ Rn×ndiagonala astfel ıncat

A = LDU,

atunci expresia de mai sus este numita factorizare LDU a matricei A.


Factorizarea LDU exista si este unica daca este satisfacuta conditia defunctionare a algoritmului G de eliminare gausiana.

Teorema 2.12 O matrice A ∈ Rn×n are o unica factorizare LDU daca sinumai daca submatricele lider principale A[k], k = 1 : n−1, sunt nesingulare.

Factorizarile LU utilizate ın practica numerica asociaza matricea diago-nala D din cu una dintre matricele U , L ale factorizarii LDU. Astfel se obtin,respectiv:

• factorizarea Doolittle A = LU , cu L inferior triunghiulara unitate si Usuperior triunghiulara,

• factorizarea Crout A = LU , cu L inferior triunghiulara si U superior triun-ghiulara unitate.

Din Teorema 2.12 rezulta ca ambele factorizari, Doolittle si Crout, existasi sunt unice daca si numai daca submatricele lider principale A[k], k = 1 :n−1, sunt nesingulare. Vom pastra ın continuare aceasta ipoteza, specificandcazurile ın care ea nu este necesara.

2.5.1 Factorizari LU rezultate din eliminareagaussiana

Factorizarea Doolittle este calculata chiar de algoritmul G de triangularizaregaussiana! Este una din ”surprizele” placute ale calculului numeric.

Teorema 2.13 Algoritmul G, aplicat matricei A ∈ Rn×n, cu detA[k] 6= 0,k = 1 : n−1, calculeaza factorizarea Doolittle A = LU , cu L inferior triunghi-ulara unitate si U superior triunghiulara. Factorii L (mai putin elementelediagonale) si U sunt memorati peste A.

Demonstratie. Algoritmul G produce matricele ITE Mk, k = 1 : n−1, astfelıncat din (2.16) rezulta A = LU , cu L = (Mn−1 . . .M2M1)−1 si U memoratapeste A. Asadar

L = M−11 M−12 . . .M−1n−1.

Dar (vezi Propozitia 2.3), M−1k = In +mkeTk , deci

L = (In +m1eT1 )(In +m2e

T2 ) . . . (In +mn−1e

Tn−1) = In +

n−1∑k=1

mkeTk + S,

unde S este o suma de matrice de forma

T = . . .mpeTpmqe

Tq . . .


cu p < q. Deoarece eTpmq = 0, pentru orice p < q, obtinem T = 0 si S = 0.Astfel,

L = In +

n−1∑k=1

mkeTk =

1 0 . . . 0 0µ21 1 . . . 0 0µ31 µ32 . . . 0 0

. . .. . . . . .

µn−1,1 µn−1,2 . . . 1 0µn1 µn2 . . . µn,n−1 1

,

deci matricea L este formata chiar din multiplicatorii gaussieni memorati ınA, ın algoritmul G.

Vom studia acum legatura ıntre algoritmul GPP de eliminare gaussianacu pivotare partiala si factorizarile LU.

Teorema 2.14 Triangularizarea (2.20) (Mn−1Pn−1 . . .M1P1A = U) a ma-tricei A ∈ Rn×n, obtinuta prin algoritmul GPP, este echivalenta cu factor-izarea Doolittle a matricei A cu liniile permutate, adica

PA = LU, (2.25)

unde P = Pn−1 . . . P1. In partea subdiagonala, coloana k ∈ 1 : n−1 a matriceiL este o permutare a vectorului Gauss corespunzator mk; mai precis, dacahk = Pn−1 . . . Pk+1mk, atunci L(k + 1 : n, k) = hk(k + 1 : n).

Demonstratie. Se poate verifica usor echivalenta dintre (2.20) si

Mn−1 . . . M1PA = U, (2.26)

unde Mn−1 = Mn−1 si, pentru k < n− 1,

Mk = Pn−1 . . . Pk+1MkPk+1 . . . Pn−1 =

= I − Pn−1 . . . Pk+1mkeTk Pk+1 . . . Pn−1 =

= I − Pn−1 . . . Pk+1mkeTk = I − hkeTk .

In relatiile de mai sus am tinut seama ca Mk = I−mkeTk si P−1s = Ps; ın plus,

eTk Ps = eTk , pentru s > k, deoarece Ps este o PE de indici s si t ≥ s. Deci Mk

este o matrice ITE definita de vectorul Gauss hk. Deoarece (2.26) si (2.16)sunt analoage, rezulta, ca ın Teorema 2.13, ca matricea L = M−11 . . . M−1n−1este inferior triunghiulara unitate cu coloanele egale ın partea subdiagonalacu vectorii hk, k ∈ 1 : n− 1. Prin urmare relatia (2.25) este satisfacuta.

Expresia vectorului hk arata cum trebuie modificat algoritmul GPP pen-tru a produce factorizarea LU cu pivotare (2.25); se observa ca asuprafiecarui vector Gauss mk trebuie aplicate toate permutarile ulterioare Ps,


cu s > k. Pentru a implementa aceasta, modificam instructiunea 1.3 a GPPdin ”Pentru j = k : n” ın ”Pentru j = 1 : n”, adica, ceea ce e echivalentcu afirmatia anterioara, fiecare permutare Pk se aplica si tuturor vectorilorGauss anteriori mj , cu j ∈ 1 : k− 1. Dupa executia algoritmului GPP astfelmodificat, matricele triunghiulare L si U din (2.25) se vor afla ın A.

Cazul eliminarii gaussiane cu pivotare completa este similar. Modificandalgoritmul GPC ın acelasi mod ca pe GPP, se va obtine pe loc ın matriceaA o factorizare LU a unei matrice obtinute din A prin permutari de linii sicoloane:

PAQ = LU, (2.27)

unde P , L si U au aceeasi semnificatie (dar nu aceleasi elemente) ca ın Teo-rema 2.14, iar Q = Q1 . . . Qn−1 este construita din matricele de permutareelementara din relatia (2.22).

2.5.2 Algoritmul Crout

Factorizarea Crout poate fi calculata luand ın relatia (2.24) expresia LU cafiind factorizarea Doolittle si D = diag(u11, u12, . . . , unn). Atunci, L′U ′ estedescompunerea Crout a matricei A.

Totusi, exista si alte posibilitati de a organiza calculele, care de obicei senumesc scheme de factorizare compacta. Vom deduce o astfel de factorizareCrout pe baza ideii de a utiliza direct identitatea A = LU , calculand succesivcate o coloana din L si cate o linie din U . Procedura poate fi initializataconsiderand prima coloana a identitatii A = LU ,

Ae1 = LUe1 = Le1,

de undeli1 = ai1, i = 1 : n, (2.28)

respectiv prima linie a aceleiasi identitati,

eT1 A = eT1 LU = l11eT1 U,

adicaa1j = l11u1j , j = 2 : n.

Deoarece l11 = a11 = A[1] 6= 0, din relatia anterioara rezulta

u1j = a1j / l11, j = 2 : n. (2.29)

Presupunem acum ca primele k − 1 coloane din L si primele k − 1 liniidin U sunt cunoscute (calculate); vezi Figura 2.6. Din identitatea A = LU


A

akk

aik

akj

utilizat neutilizat

=

L

lkk

liki

calculat necalculat

0

U

1 ukj

calculat

necalculat

0

Fig. 2.6 Stadiu intermediar al factorizarii Crout.

obtinem

aik =

k−1∑s=1

lisusk + lik, i = k : n,

ın care singura necunoscuta este scalarul lik. Deci,

lik = aik −k−1∑s=1

lisusk, i = k : n. (2.30)

De asemenea,

akj =

k−1∑s=1

lksusj + lkkukj , j = k + 1 : n,

ın care acum singura necunoscuta este scalarul ukj . Din Propozitia 2.6 rezultaA[k] = L[k]U [k]; avand ın vedere ca U [k] este nesingulara prin constructie, dacaA[k] este nesingulara, atunci si L[k] este nesingulara, adica lkk 6= 0. Asadar,

ukj =

(akj −

k−1∑s=1

lksusj

)/ lkk, j = k + 1 : n. (2.31)

Astfel, procedura initiata prin (2.28) si (2.29) poate fi continuata prin(2.30) si (2.31); daca submatricele A[k], k = 1 : n − 1, sunt nesingulare,ceea ce garanteaza ca ımpartirea din (2.31) poate fi executata, procedura vacalcula complet factorii L si U . Observand ca (2.28) si (2.29) sunt cazuriparticulare ale (2.30), respectiv (2.31), putem include primul pas ın buclagenerala, rezultand algoritmul urmator.

Algoritmul 2.8 (CROUT – Factorizare Crout) (Se da A ∈ Rn×n,cu submatricele lider principale A[k], k = 1 : n − 1, nesingulare. Se


calculeaza factorizarea Crout A = LU , cu L inferior triunghiulara si Usuperior triunghiulara unitate. Matricele L si U sunt memorate pesteelementele corespunzatoare din A, mai putin elementele diagonale alelui U , egale cu 1.)

1. Pentru k = 1 : n1. Pentru i = k : n

1. aik ← lik = aik −∑k−1s=1 lisusk

2. Pentru j = k + 1 : n

1. akj ← ukj =(akj −

∑k−1s=1 lksusj

)/ lkk

Comentarii. Numarul de operatii ın virgula mobila esteNCROUT ≈ 2n3/3 lafel ca ın algoritmii G, GPP sau GPC, iar memoria ocupata este MCROUT ≈n2.

Ca si algoritmul G de eliminare gaussiana, algoritmul CROUT esueaza peparcurs daca se obtine un lkk = 0 (ceea ce este echivalent cu A[k] singulara);mai mult, valori apropiate de zero ale unui lkk conduc la valori absolute foartemari ale elementelor matricei U calculate ın instructiunea 1.2.1 precum si laerori numerice inadmisibile ın rezultat.

Pentru a preveni aceasta situatie se introduce pivotarea partiala (pe linii).Spre deosebire de algoritmul G unde, la pasul k, permutarea liniilor se efec-tua ınaintea celorlaltor operatii, aici permutarea poate avea loc doar dupacalcularea elementelor lik, cu i ∈ k : n, pentru a putea aduce cea maiconvenabila valoare ın pozitia (k, k). (Sa observam ca pivotarea completaeste acum imposibila, deoarece ar implica anticiparea calculelor pentru toatamatricea.) Aceasta revine la ınmultirea la stanga a matricei A cu o PEPk ≡ Pkik , cu ik ≥ k, la fiecare pas k si, global, la o factorizare PA = LU ,cu P = Pn−1 . . . P1. Algoritmul corespunzator este urmatorul.

Algoritmul 2.9 (CROUTP – Factorizare Crout cu pivotare partia-la) (Se da A ∈ Rn×n nesingulara. Se calculeaza factorizarea Crout PA =LU , cu L inferior triunghiulara, U superior triunghiulara unitate siP = Pn−1 . . . P1 matrice de permutare. Matricele L si U sunt memoratepeste elementele corespunzatoare din A, mai putin elementele diagonaleale lui U , egale cu 1. Intregii ik definind PE Pk sunt memorati ıntr-unvector p astfel ıncat p(k) = ik, k = 1 : n− 1.)

1. Pentru k = 1 : n1. Pentru i = k : n

1. aik ← lik = aik −∑k−1s=1 lisusk

2. Se determina ik ∈ k : n astfel ıncat |likk| = maxi=k:n |lik|3. Pentru j = 1 : n {interschimba liniile k si ik}

1. akj ↔ aikj4. Pentru j = k + 1 : n

1. akj ← ukj =(akj −

∑k−1s=1 lksusj

)/ lkk

2.6 Rezolvarea sistemelor liniare 43

Comentarii. Desigur, numarul comparatiilor efectuate ın instructiunea 1.2este de ordinul O(n2), adica neglijabil fata de cel al operatiilor aritmetice.

2.6 Rezolvarea sistemelor liniare

Teoremele 2.8 si 2.9 arata ca algoritmii GPP si GPC constituie un buninstrument pentru rezolvarea sistemului liniar nesingular Ax = b.

Sa consideram urmatoarea sintaxa neformala pentru utilizarea algoritmilorGPP si respectiv GPC:

[M,U, p] = GPP(A)[M,U, p, q] = GPC(A)

unde M denota multiplicatorii gaussieni µik, k = 1 : n−1, i = k+1 : n, U estematricea superior triunghiulara produsa de algoritm, iar p si q sunt vectoride ıntregi definind permutarile de linii si, respectiv, coloane. Desi matriceleM si U se memoreaza peste A (asa cum am si procedat ın GPP si GPC),vom utiliza sintaxa de mai sus pentru claritate.

Rezolvarea sistemelor utilizand GPP. Consideram sistemul liniarAx = b, cu A nesingulara si utilizam GPP pentru triangularizarea matriceiA conform relatiei (2.20). Doarece toate matricele Mk si Pk sunt nesingulare,sistemul Ax = b este echivalent (are aceeasi solutie) cu sistemul

Mn−1Pn−1 . . .M1P1Ax = Mn−1Pn−1 . . .M1P1b. (2.32)

In consecinta, (2.32) se scrie ın forma

Ux = d,

unde matricea superior triunghiulara U se afla ın A, iar vectorul

d = Mn−1Pn−1 . . .M1P1b, (2.33)

se poate calcula pe loc ın b, utilizand urmatoarea schema de calcul:

b 1. Pentru k = 1 : n− 11. b← Pkb2. b←Mkb

In final, solutia x a sistemului Ax = b se obtine rezolvand (pe loc ın b)sistemul superior triunghiular Ux = b. Evident, ın acest scop se utilizeazaalgoritmul UTRIS, cu sintaxa

x = UTRIS(U, b).


Pentru detalierea schemei facem apel la Propozitiile 2.3 si 2.5. Se obtinealgoritmul de mai jos.

Algoritmul 2.10 (S GPP – Calculul solutiei unui sistem liniar uti-lizand GPP) (Se dau A ∈ Rn×n nesingulara si b ∈ Rn. Se calculeazasolutia x ∈ Rn a sistemului liniar Ax = b, utilizand procedura de elimi-nare gaussiana cu pivotare partiala.)

1. [M,U, p] = GPP(A)2. Pentru k = 1 : n− 1

1. bk ↔ bp(k)2. Pentru i = k + 1 : n

1. bi ← bi − µikbk3. x = UTRIS(U, b)

Comentarii. Numarul de operatii este:

NS GPP = NGPP+

n−1∑k=1

2(n− k)+NUTRIS ≈2n3

3+n2+n2 ≈ 2n3

3≈ NGPP,

si deci principalul efort consta ın triangularizarea matricei A. Evident, memo-ria necesara este MS GPP ≈ n2.

Rezolvarea sistemelor utilizand GPC. Alternativ, dar mai costisitor,putem utiliza GPC, vezi (2.22); ın acest caz, sistemul Ax = b este echivalentcu

Mn−1Pn−1 . . .M1P1AQ1 . . . Qn−1Qn−1 . . . Q1x = Mn−1Pn−1 . . .M1P1b.(2.34)

NotandQn−1 . . . Q1x = y, (2.35)

din (2.22) rezulta ca vectorul y poate fi calculat prin rezolvarea sistemuluisuperior triunghiular

Uy = d, (2.36)

unde d este vectorul dat de (2.33). In final, din (2.35), si stiind ca Q−1k = Qk,obtinem

x = Q1Q2 . . . Qn−1y. (2.37)

Din relatiile (2.33)–(2.37), folosind aceeasi variabila pentru x si y, rezultaalgoritmul urmator.

Algoritmul 2.11 (S GPC – Calculul solutiei unui sistem liniar uti-lizand GPC) (Se dau A ∈ Rn×n nesingulara si b ∈ Rn. Se calculeazasolutia x ∈ Rn a sistemului liniar Ax = b, utilizand procedura de elim-inare gaussiana cu pivotare completa.)

2.6 Rezolvarea sistemelor liniare 45

1. [M,U, p, q] = GPC(A)2. Pentru k = 1 : n− 1

1. bk ↔ bp(k)2. Pentru i = k + 1 : n

1. bi ← bi − µikbk3. x = UTRIS(U, b)4. Pentru k = n− 1 : −1 : 1

1. xk ↔ xq(k)

Comentarii. Complexitatea acestui algoritm este asimptotic aceeasi cu ceaa algoritmului GPC.

Observatia 2.15 Eliminarea gaussiana, prezentata aici ıntr-un mod sistem-atic, este binecunoscuta metoda de calcul manual prin reducere si substitutie.Pivotarea pe linii corespunde unei reordonari a ecuatiilor ın (2.2), iar piv-otarea pe coloane corespunde unei reordonari a componentelor necunoscuteix. Subliniem ca aceasta metoda este cea mai eficienta procedura de rezolvarea sistemelor de ecuatii liniare nesingulare de forma generala.

Observatia 2.16 Algoritmii GPP si GPC au proprietati foarte bune ınprivinta stabilitatii numerice. Algoritmul GPP nu este neconditionat stabil;au fost construite artificial exemple de matrice pentru care este instabil. Cutoate acestea, pentru probleme practice, GPP poate fi considerat stabil, deaceea este recomandat ca prima optiune pentru rezolvarea sistemelor liniare.Algoritmul GPC este numeric stabil, deci poate fi utilizat pentru orice ma-trice; din motive de complexitate, el este recomandat doar ın cazul ın care sedoreste siguranta absoluta (iar matricea A nu se ıncadreaza ıntr-o categoriecu proprietati cunoscute).

Rezolvarea sistemelor utilizand factorizarea LU. Este clar ca sis-temul Ax = b se poate rezolva si recurgand explicit la factorizarea LU amatricei A. Presupunand, pentru simplitate, ca nu a fost utilizata pivotarea,factorizarea este A = LU . Avem LUx = b, sau Ly = b, unde am notatUx = y. Prin urmare, solutia sistemului poate fi gasita dupa urmatoareaschema:

Ax = b 1. A = LU (calculeaza factorizarea LU)2. Se rezolva sistemul inferior triunghiular Ly = b3. Se rezolva sistemul superior triunghiular Ux = y

Efortul de calcul este acelasi ca ın eliminarea gaussiana si partea sa ceamai importanta se consuma ın procesul de factorizare. Desigur, ın practicaeste recomandat a se folosi pivotarea, efectuand o factorizare LU de formaPA = LU .

Rezolvarea sistemelor cu parte dreapta multipla. Daca sistemulliniar este de forma AX = B cu B ∈ Rn×p, p > 1, atunci X rezulta de


aceeasi dimensiune ca B; sistemul AX = B poate fi redus la p sisteme cu partedreapta simpla, Axj = bj , j = 1 : p, unde vectorii bj si xj sunt coloanele j dinB, respectiv X. Utilizarea naiva de p ori a algoritmului S GPP ar conducela un efort de calcul de pn3/3 flopi, inacceptabil de mare. Dar, deoareceA este aceeasi pentru toate cele p sisteme, factorizarea LU sau, ın general,triangularizarea matricei A, se efectueaza o singura data, ramanand apoi derezolvat doar sisteme triunghiulare. Ideea algoritmului, ın cazul factorizariiLU, este urmatoarea:

AX = B 1. A = LU (calculeaza factorizarea LU)2. Pentru j = 1 : p

1. Se rezolva sistemul inferior triunghiular Ly = bj2. Se rezolva sistemul superior triunghiular Uxj = y

Numarul de operatii va de numai aproximativ 2n3/3 + 2pn2 flopi, menti-nandu-se complexitatea O(n3). Scrieti varianta schemei de mai sus ın carese utilizeaza GPP pentru triangularizarea matricei A, modificand adecvatalgoritmul S GPP (instructiunea 1 se executa o singura data, celelalte de pori, cate o data pentru fiecare coloana din B si X).

Rezolvarea sistemelor cu necunoscuta ın stanga. Un sistem liniarpoate avea si forma XA = B cu B ∈ Rp×n, solutia X avand din nou aceeasidimensiune ca X. Rezolvarea unui astfel de sistem nu necesita un algoritmspecial, deoarece prin transpunere se obtine ATXT = BT , adica un sistemcu necunoscuta aflata ın dreapta, ca de obicei. (Pozitionarea necunoscutei ındreapta este deci o alegere arbitrara, adica este o simpla conventie.)

2.7 Factorizarea Cholesky

Consideram acum sistemul Ax = b, unde matricea A ∈ Rn×n este simetricasi pozitiv definita.

Definitia 2.17 Matricea A este simetrica daca A = AT , adica aij = aji,∀i, j ∈ 1 : n.

Asadar simetria este fata de diagonala principala a matricei; se observa capentru o matrice nepatrata notiunea de simetrie nu are sens, A si AT neavandaceleasi dimensiuni.

Definitia 2.18 O matrice simetrica A este pozitiv definita daca xTAx > 0,pentru orice vector x ∈ Rn, x 6= 0; se noteaza A � 0. Matricea A este pozitivsemidefinita daca xTAx ≥ 0, ∀x ∈ Rn, ceea ce se noteaza A � 0.

Atentie, faptul ca A � 0 nu ınseamna ca toate elementele matricei sunt

pozitive; de exemplu,

[2 −1−1 2

]este pozitiv definita, deoarece 2x21− 2x1x2 +

2.7 Factorizarea Cholesky 47

2x22 > 0, ∀[x1 x2]T 6= 0. Pe de alta parte, matricea

[1 22 1

]nu este pozitiv

definita, deoarece x21 − 4x1x2 + x22 poate avea orice semn cand x1, x2 ∈ R.

Pentru matrice simetrice si pozitiv definite factorizarea LU capata formadescrisa mai jos.

Teorema 2.19 Pentru orice matrice simetrica si pozitiv definita A ∈ Rn×n,exista o unica matrice inferior triunghiulara L ∈ Rn×n, cu elemente diago-nale pozitive, astfel ıncat

A = LLT , (2.38)

expresie care este numita factorizare Cholesky (iar L este numit factorulCholesky al matricei A). Reciproc, daca factorizarea Cholesky exista, atunciA este pozitiv definita.

Evident, (2.38) se poate scrie sub forma A = RTR, unde R = LT estesuperior triunghiulara. De asemenea, factorizarile A = LTL, respectiv A =RRT , se numesc tot Cholesky.

Demonstratie. Deoarece A este pozitiv definita atunci submatricele liderprincipale A[k], k = 1 : n − 1, sunt pozitiv definite, deci nesingulare. AtunciA are o unica factorizare LDU (vezi Teorema 2.12), A = L′DU ′ si, din motivede simetrie, U ′ = (L′)T . Deci A si D sunt congruente, rezultand ca D estepozitiv definita, adica dii > 0, i = 1 : n. Fie F ∈ Rn×n matricea diagonaladefinita de fii =

√dii, i = 1 : n. Notand L = L′F , prima parte a teoremei

este demonstrata (lii =√dii, i = 1 : n).

Reciproca este evidenta; fie x ∈ Rn nenul, altfel arbitrar; atunci avemxTAx = xTLLTx = ‖LTx‖22 > 0 (deoarece L este nesingulara, LTx 6= 0).

Exista mai multi algoritmi pentru calculul factorizarii Cholesky a uneimatrice pozitiv definite, care difera doar prin ordinea de efectuare a calculelor.Vom deduce aici o procedura de factorizare compacta, similara cu cea folositala calculul factorizarii Crout, ın care elementele factorului Cholesky L suntcalculate ın ordinea (crescatoare a) coloanelor. Din identitatea A = LLT

(doar triunghiul inferior stang al lui A este reprezentat, deoarece A estesimetrica)

a11...

. . .

ak1 . . . akk...

.... . .

an1 . . . ank . . . ann

=

l11...

. . . 0lk1 . . . lkk...

.... . .

ln1 . . . lnk . . . lnn

l11 . . . lk1 . . . ln1. . .

......

lkk . . . lnk

0. . .

...lnn

obtinem pentru prima coloana:

a11 = l211 ⇒ l11 =√a11,

ai1 = li1l11 ⇒ li1 = ai1 / l11, i = 2 : n,(2.39)


si deci procesul de calcul poate fi initiat.

Acum, presupunand ca primele k − 1 coloane din L au fost calculate,identificarea elementului din pozitia (k, k) ın A = LLT furnizeaza

akk =

k−1∑j=1

l2kj + l2kk, (2.40)

ın care singura necunoscuta este lkk. Conform Teoremei 2.19, daca A estepozitiv definita, atunci exista un unic lkk pozitiv astfel ıncat (2.40) sa aibaloc. Atunci,

akk −k−1∑j=1

l2kj > 0 (2.41)

si

lkk =

√√√√akk −k−1∑j=1

l2kj . (2.42)

De asemenea,

aik =

k−1∑j=1

lij lkj + liklkk, i = k + 1 : n,

ın care, ın ordinea de calcul considerata, singura necunoscuta este lik. Asadar

lik =

aik − k−1∑j=1

lij lkj

/ lkk, i = k + 1 : n. (2.43)

Formulele (2.39), (2.42) si (2.43) definesc algoritmul de mai jos; pentruconcizie, am inclus si primul pas ın bucla de parcurgere a coloanelor. Toto-data, acest algoritm este cel mai bun mijloc de testare a pozitivitatii uneimatrice; testul se bazeaza pe faptul ca factorizarea Cholesky a unei matricepoate fi calculata daca si numai daca matricea este pozitiv definita (Teorema2.19); asadar, daca termenul stang al inegalitatii (2.41) nu este pozitiv, atuncimatricea A nu este pozitiv definita.

Algoritmul 2.12 (CHOL – factorizare Cholesky) (Se da A ∈ Rn×nsimetrica. Se stabileste daca A este pozitiv definita si, ın acest caz,se scrie peste triunghiul inferior al lui A matricea L din factorizareaCholesky A = LLT .)

1. Pentru k = 1 : n

1. α← akk −∑k−1j=1 l

2kj

2. Daca α ≤ 0 atunci

1. Tipareste ’A nu este pozitiv definita’2. Stop

2.8 Calculul determinantului 49

3. akk ← lkk =√α

4. Pentru i = k + 1 : n

1. aik ← lik =(aik −

∑k−1j=1 lij lkj

)/ lkk

Comentarii. Algoritmul CHOL necesita aproximativ NCHOL = n3/3 flopisi, ın plus, calculul a n radacini patrate (care e neglijabil). Memoria necesaraeste de aproximativ MCHOL = n2/2.

Algoritmul este numeric stabil.

Desigur, solutia sistemului liniar Ax = b, cu A pozitiv definita se calculeazacu schema

Ch 1. Se factorizeaza A = LLT cu algoritmul CHOL2. Se rezolva sistemul inferior triunghiular Ly = b3. Se rezolva sistemul superior triunghiular LTx = y

efectuandu-se de doua ori mai putine operatii decat daca s-ar aplica algorit-mul de eliminare gaussiana.

2.8 Calculul determinantului

Determinantul matricei A ∈ Rn×n este un numar real, notat det(A), care

poate fi definit recursiv. Daca Adef= a ∈ R, atunci det(A) = a. Pentru A ∈

Rn×n,

det(A) =

n∑j=1

(−1)j+1 a1j det(A1j),

unde A1j ∈ R(n−1)×(n−1) este matricea obtinuta din A prin eliminarea liniei1 si coloanei j.

Cateva proprietati importante ale determinantului sunt:

1. det(I) = 1;

2. det(αA) = αndet(A), ∀α ∈ R;

3. det(AT ) = det(A);

4. det(AB) = det(A) · det(B).

Precizam ca, din punct de vedere numeric, calculul determinantului nu esteutil pentru determinarea nesingularitatii unei matrice. Un algoritm fiabil vafi prezentat ın capitolul 5.

Definitia recursiva de mai sus nu conduce la algoritmi eficienti pentrucalculul determinantului. Algoritmii GPP si GPC sunt utili ın acest scop.

Amintim ca interschimbarea a doua linii sau coloane ale unei matriceschimba semnul determinantului acesteia, sau, echivalent, daca P este o


permutare elementara, atunci det(P ) = −1. In plus, matricele ITE Mk,k = 1 : n − 1, din (2.20), au elementele diagonale egale cu 1 si decidet(Mk) = 1. Asadar

det(A) = (−1)s det(U) = (−1)sn∏i=1

uii,

unde s ≤ n− 1 este numarul de PE proprii (adica Pk 6= In) ın (2.20).

Algoritmul corespunzator celor de mai sus este:

Algoritmul 2.13 (DET GPP – Calculul determinantului) (Se daA ∈ Rn×n. Se calculeaza detA, utilizand algoritmul GPP.)

1. [M,U, p] = GPP(A)2. det← 13. Pentru k = 1 : n

1. det← det · ukk4. Pentru k = 1 : n− 1

1. Daca p(k) 6= k atunci

1. det← −det

Principalul efort de calcul ın algoritmul de mai sus este efectuat ın executiaalgoritmului GPP. Desigur, algoritmul GPC poate fi utilizat ın loc de GPP,obtinandu-se, de regula, o acuratete superioara a rezultatului. Scrierea algo-ritmului corespunzator este propusa cititorului.

Daca este disponibila o factorizare LU a matricei A, atunci det(A) =det(LU) = = det(L) · det(U) = (

∏ni=1 lii) (

∏ni=1 uii), ceea ce permite evalu-

area imediata a determinantului.

2.9 Calculul inversei (optional)

Precizam de la ınceput ca, ın majoritatea cazurilor, calculul explicit al inverseiunei matrice nu este necesar, el putand de multe ori fi ınlocuit eficient curezolvarea unui sistem liniar.

Daca A ∈ Rn×n este nesingulara, atunci procedura de eliminare gaus-siana cu pivotare partiala se termina cu succes si produce o matrice superiortriunghiulara nesingulara U astfel ıncat

Mn−1Pn−1 . . .M1P1A = U, (2.44)

unde Mk sunt matrice ITE si Pk sunt matrice elementare de permutare. Din(2.44) rezulta

X ≡ A−1 = U−1Mn−1Pn−1 . . .M1P1,

deci inversa matricei A poate fi calculata astfel

2.9 Calculul inversei (optional) 51

A← A←M

U

M

U−1

=⇒

Fig. 2.7 Memorarea datelor ın inversarea matriceala

A−1 1. [M,U, p] = GPP(A)2. X = UINV(U)3. Pentru k = n− 1 : −1 : 1

1. X ← XMk

2. X ← XPk

Algoritmul UINV pentru calculul inversei unei matrice superior triunghi-ulare va fi discutat mai tarziu, tot an aceasta sectiune.

In algoritmul GPP, matricea U si multiplicatorii gaussieni µik, k = 1 :n− 1, i = k+ 1 : n, notati global cu M ın schema de mai sus, se memoreazapeste matriceaA ca ın partea stanga a Figurii 2.7. Vom vedea ca, ın algoritmulUINV, U−1 se poate memora peste U , astfel ca, dupa instructiunea 2 aschemei de mai sus, datele sunt memorate ca ın partea dreapta a Figurii 2.7.

Pornind de la schema precedenta, vom prezenta un algoritm care calculeazainversa matricei A, folosind minimum de memorie. Pentru aceasta, observamca primele doua instructiuni necesita memorie aditionala doar pentru vectorulp ∈ Rn−1, ın care se memoreaza permutarile de linii din GPP. Analizamacum bucla Pentru din instructiunea 3. Prima atribuire, X ← XMk, poatefi detaliata utilizand o partitionare pe linii a lui X:

xTi Mk = xTi (In −mkeTk ) = xTi − (xTi mk)eTk (2.45)

= [xi1 xi2 . . . xi,k−1 xik − xTi mk xi,k+1 . . . xin], i = 1 : n.

Se poate observa ca sunt modificate doar elementele de pe coloana k a lui X.

A doua instructiune a buclei Pentru, X ← XPk, interschimba coloanele ksi ik (ik ≥ k).

Tinand cont de toate aceste consideratii si de ordinea inversa de cal-cul a coloanelor, pasul k al buclei Pentru din schema de mai sus, adicaX ← XMkPk, poate fi descris ca ın Figura 2.8, unde U−1 si M reprezintaelementele din U−1 si, respectiv, multiplicatorii gaussieni care nu au fostutilizati (si nici modificati) ın executia buclei pana la pasul curent; elementelemodificate (parte a lui X) sunt notate prin ×.

Astfel, singura informatie care trebuie salvata consta ın multiplicatoriigaussieni din coloana curenta, µik, i = k+ 1 : n, necesari pentru calculul din(2.45).


u(−1)1k

...u(−1)kk

µk+1,k

...µnk

× . . . ×

× . . . ×

× . . . ×

. . .

. . .

M

U−1

k

@@@@@ X ← XMkPk

×...×

×...×

× . . . ×

× . . . ×

× . . . ×

. . .

. . .

M

U−1

k

@@@@@-

Fig. 2.8 Memorarea datelor ın pasul k al inversarii matriceale

In acest fel, inversarea matricei A poate fi executata utilizand doar doi vec-tori suplimentari de lungime (n−1), anume p pentru memorarea permutarilorsi, sa zicem, g pentru salvarea temporara a multiplicatorilor gaussieni de pecoloana curenta.

Algoritmul care implementeaza consideratiile de mai sus este:

Algoritmul 2.14 (INV GPP – Calculul inversei unei matrice) (Seda A ∈ Rn×n nesingulara. Se calculeaza matricea A−1 peste A. Algo-ritmul utilizeaza algoritmii GPP pentru triangularizarea gaussiana siUINV pentru inversarea unei matrice superior triunghiulare.)

1. [M,U, p] = GPP(A) % M si U sunt memorate ın A2. A← X = UINV(U) % doar triunghiul sup. e modificat3. Pentru k = n− 1 : −1 : 1

1. Pentru i = k + 1 : n1. gi ← µik % salveaza multiplicatorii

2. Pentru i = 1 : k1. aik ← aik −

∑nt=k+1 aitgt% relatia (2.45), primele k linii

3. Pentru i = k + 1 : n1. aik ← −

∑nt=k+1 aitgt % relatia (2.45), celelalte linii

4. Daca p(k) 6= k atunci

1. Pentru i = 1 : n % interschimba coloanele k si p(k)1. aik ↔ ai,p(k)

Comentarii. Numarul de operatii ın virgula mobila este

NINV GPP = NGPP +NUINV +

n−1∑k=1

2n(n− k) ≈ 2n3

3+n3

3+ n3 = 2n3

si memoria utilizata MINV GPP = n2.

Este remarcabil faptul ca inversarea matriceala nu este mai complexa (canumar de operatii) decat ınmultirea de matrice.


O acuratete mai buna a matricei inverse calculate se poate obtine utilizandalgoritmul GPC ın locul lui GPP. In acest caz,

Mn−1Pn−1 . . .M1P1AQ1Q2 . . . Qn−1 = U,

si deciX ≡ A−1 = Q1Q2 . . . Qn−1U

−1Mn−1Pn−1 . . .M1P1.

Scrierea algoritmului corespunzator este lasata ın sarcina cititorului.

In practica, acuratetea obtinuta cu INV GPP este suficient de buna pen-tru majoritatea aplicatiilor.

Observatia 2.20 Efortul de calcul necesar inversarii unei matrice este deaproximativ trei ori mai mare decat cel pentru rezolvarea unui sistem liniarde aceeasi dimensiune. De aceea, de cate ori e posibil, inversarea matricealatrebuie evitata si ınlocuita de rezolvarea de sisteme liniare. De exemplu, pen-tru calculul scalarului real

α = cTA−1b,

unde A ∈ Rn×n, b, c ∈ Rn, schema recomandata este

cTA−1b 1. Se rezolva sistemul Ax = b % rezulta x = A−1b2. α← cTx

sensibil mai eficienta decat varianta cu inversarea matricei A

Observatia 2.21 Este clar ca inversa matricei A se poate calcula si re-curgand explicit la factorizarea LU a matricei. Presupunand ca A = LU ,avem

A−1 = U−1L−1,

iar o procedura de inversare ar putea consta ın:

1. A = LU (se calculeaza factorizarea LU)2. U ← X = U−1

3. L← Y = L−1

3. A−1 ← XY

Numarul de operatii este tot 2n3. Desigur, ın practica este indicat a sefolosi o factorizare LU cu pivotare.

Inversarea matricelor triunghiulare. Fie L ∈ Rn×n o matrice inferiortriunghiulara; presupunem ca L este inversabila, adica lii 6= 0, i ∈ 1 : n.Pentru a calcula X = L−1, utilizam egalitatea evidenta

LX = In,

pe care o partitionam pe coloane


Lxj = ej , j ∈ 1 : n, (2.46)

xj = Xej fiind cea de-a j-a coloana a matricei X. Sistemul inferior triunghiu-lar (2.46) poate fi rezolvat ın mod eficient prin adaptarea algoritmului LTRISla forma particulara a termenului drept. Intr-adevar, sistemul (2.46) poate fipartitionat astfel: [

L(j)11 0

L(j)21 L

(j)22

] [x′jx′′j

]=

[0e′′j

], (2.47)

unde L(j)11 este submatricea de dimensiune (j − 1)× (j − 1) din coltul stanga

sus al matricei L sie′′j = [1 0 . . . 0]T ∈ Rn−j+1.

Din (2.47) rezulta {L(j)11 x

′j = 0,

L(j)21 x

′j + L

(j)22 x

′′j = e′′j .

(2.48)

Deoarece L este inversabila si, prin urmare, la fel sunt L(j)11 si L

(j)22 , relatia

(2.48) devine {x′j = 0,

L(j)22 x

′′j = e′′j .

(2.49)

Aceasta relatie stabileste un fapt important, dat de

Propozitia 2.22 Inversa unei matrice inferior triunghiulare nesingulareeste inferior triunghiulara.

In concluzie, matricele inferior triunghiulare inversabile formeaza grup ınraport cu ınmultirea matriceala.

Din (2.49) poate fi dedus un algoritm de calcul al inversei unei matricenesingulare, avand la baza urmatoarea schema de calcul:

L−1 1. Pentru j = 1 : n1. Daca j > 1 atunci x′j = 0

2. x′′j = LTRIS(L(j)22 , e

′′j )

Se poate observa ca inversa calculata poate fi memorata ın aceeasi zona dememorie ca matricea L. Explicitand rezolvarea sistemului inferior triunghiu-lar, obtinem

Algoritmul 2.15 (LINV – Calculul inversei unei matrice inferiortriunghiulare) (Se da L ∈ Rn×n, inferior triunghiulara si inversabila. Secalculeaza pe loc ın L inversa acesteia.)

1. Pentru j = 1 : n1. ljj ← xjj = 1 / ljj


2. Pentru i = j + 1 : n

1. lij ← xij = −(∑i−1

k=j likxkj

)/ lii

Comentarii. Numarul de operatii necesar este

NLINV = n+

n−1∑j=1

n∑i=j+1

2(i− j) ≈ n3

3.

Memoria necesara este MLINV ≈ n2/2 locatii ın virgula mobila. Se vedeacum ca rezolvarea sistemului liniar Lx = b folosind relatia x = L−1b (adicaLINV) este ineficienta ın comparatie cu algoritmul direct LTRIS.

Pentru a calcula inversa unei matrice superior triunghiulare nesingulareU ∈ Rn×n, trebuie rezolvata ecuatia matriceala UX = In, sau, echivalent,cele n ecuatii liniare

Uxj = ej , j = 1 : n, (2.50)

unde xj este coloana j a matricei X = U−1.

Utilizand partitia [U

(j)11 U

(j)12

0 U(j)22

] [x′jx′′j

]=

[e′j0

], (2.51)

unde U(j)11 este submatricea de dimensiune j×j din coltul stanga sus a lui U si

e′j = [0 0 . . . 0 1]T ∈ Rj , si urmand o cale asemanatoare celei parcurse pentruinversarea matricelor inferior triunghiulare, se obtin analoagele relatiei (2.49),Propozitiei 2.22 si algoritmului LINV. Avem:{

U(j)11 x

′j = e′j ,

x′′j = 0.(2.52)

Propozitia 2.23 Inversa unei matrice superior triunghiulare nesingulareeste superior triunghiulara.

Asadar matricele superior triunghiulare inversabile formeaza grup ın ra-port cu ınmultirea.

Daca sistemele (2.52) sunt rezolvate ın ordine inversa (pentru j = n, n −1, . . . , 1), se poate vedea cu usurinta ca elementele matricei inverse se potmemora peste cele ale matricei originale, pe masura ce sunt calculate. Vomobtine deci urmatorul algoritm.

Algoritmul 2.16 (UINV – Calculul inversei unei matrice superiortriunghiulare) (Se da U ∈ Rn×n, superior triunghiulara si inversabila.Se calculeaza pe loc ın U inversa acesteia.)


1. Pentru j = n : −1 : 11. ujj ← xjj = 1/ujj2. Pentru i = j − 1 : −1 : 1

1. uij ← xij = −(∑j

k=i+1 uikxkj

)/uii

Comentarii. La fel ca ın cazul inferior triunghiular, NUINV ≈ n3/3 siMUINV ≈ n2/2. Algoritmii de inversare a matricelor triunghiulare sunt nu-meric stabili.

2.10 Conditionarea sistemelor liniare

Introducem ıntai cateva notiuni privind normele matriceale. Acestea se de-finesc la fel ca normele vectoriale.

Definitia 2.24 O norma matriceala este o functie ‖ · ‖ : Rm×n → R+ caresatisface conditiile

1. ‖A‖ > 0, ∀A ∈ Rm×n, A 6= 0 (pozitivitate);

2. ‖αA‖ = |α| · ‖A‖, ∀A ∈ Rm×n, ∀α ∈ R (omogenitate);

3. ‖A+B‖ ≤ ‖A‖+ ‖B‖, ∀A,B ∈ Rm×n (inegalitatea triunghiului).

In plus, o norma matriceala este consistenta daca

‖AB‖ ≤ ‖A‖ · ‖B‖, (2.53)

atunci cand produsul AB este definit.

Norma Frobenius este definita prin

‖A‖F =

√√√√ m∑i=1

n∑j=1

a2ij (2.54)

si este consistenta.

Normele induse sunt definite pornind de la o norma vectoriala ‖ · ‖ prinrelatia

‖A‖ = supx6=0

‖Ax‖‖x‖

= max‖x‖=1

‖Ax‖, (2.55)

ın care norma matriceala ‖A‖ se defineste cu ajutorul normelor vectoriale ‖x‖si ‖Ax‖. (Trecerea de la sup la max este justificata de faptul ca hipersferade ecuatie ‖x‖ = 1 este un compact iar norma vectoriala ‖ · ‖ este o functiecontinua, deci marginita si ısi atinge marginile.)

2.10 Conditionarea sistemelor liniare 57

Astfel, folosind (2.55), se pot defini normele p matriceale subordonatenormelor p (Holder) vectoriale. Cele mai utilizate norme p sunt normele 1, 2si ∞. Pentru p = 1,∞, relatia (2.55) se reduce la expresiile simple

‖A‖1 = max‖x‖1=1

‖Ax‖1 = maxj=1:n

(m∑i=1

|aij |

),

‖A‖∞ = max‖x‖∞=1

‖Ax‖∞ = maxi=1:m

n∑j=1

|aij |

. (2.56)

Pentru norma 2, atat semnificatia cat si modul de calcul vor fi prezentateın capitolul 5.

Normele p matriceale, aplicate unor vectori priviti ca matrice, dau aceleasivalori ca normele p vectoriale. Astfel, nu e nici o posibilitate de confuzie ınutilizarea aceluiasi simbol ‖ · ‖p pentru normele matriceale si cele vectoriale.

In sfarsit, normele matriceale subordonate unor norme vectoriale sunt con-sistente (demonstrati, folosind definitia) si, ın plus, au proprietatea ‖I‖ = 1.

Conditionare. Dupa cum am vazut ın capitolul 1, conditionarea uneiprobleme de calcul, ın cazul nostru rezolvarea sistemului liniar determinatnesingular Ax = b, este cruciala ın aprecierea acuratetii rezultatelor obtinute.Sa consideram o perturbatie a datelor de intrare ale sistemului

(A, b) −→ (A+∆A, b+∆b) ≡ (A, b),

care va avea ca efect o perturbatie a solutiei x

x −→ x+∆x ≡ x, (2.57)

astfel ıncat(A+∆A)(x+∆x) = b+∆b. (2.58)

Consideram doar perturbatii mici, respectand ‖∆A‖ � ‖A‖, ‖∆b‖ � ‖b‖(unde ‖ · ‖ este o norma consistenta oarecare), adica efectuam o analiza asensibilitatii locale a solutiei ın raport cu variatia datelor.

Teorema 2.25 Presupunem ca x 6= 0. Eroarea relativa a solutiei (2.57) asistemului perturbat (2.58) este marginita conform relatiei

‖∆x‖‖x‖

≤ κ(A)

(‖∆A‖‖A‖

+‖∆b‖‖b‖

), (2.59)

undeκ(A) = ‖A−1‖ · ‖A‖ (2.60)

se numeste numar de conditionare al matricei A.


Cu alte cuvinte, erorile relative ‖∆A‖ / ‖A‖ si ‖∆b‖ / ‖b‖ ale datelor deintrare pot fi amplificate de κ(A) ori ın solutie, presupunand ca toate calculelese efectueaza exact (deci indiferent de algoritmul folosit).

Demonstratie. Deoarece Ax = b, din (2.58) rezulta

A ·∆x+∆A · x+∆A ·∆x = ∆b.

Neglijand produsul ∆A ·∆x (deoarece perturbatiile sunt mici), obtinem

∆x ≈ −A−1 ·∆A · x+A−1 ·∆b.

Utilizand o familie consistenta de norme matriceale ‖ · ‖, din relatia de maisus rezulta

‖∆x‖ ≤ ‖A−1‖ · ‖∆A‖ · ‖x‖+ ‖A−1‖ · ‖∆b‖.

Prin urmare, daca x 6= 0, o limita superioara a erorii relative este data de

‖∆x‖‖x‖

≤ ‖A−1‖ · ‖A‖ · ‖∆A‖‖A‖

+ ‖A−1‖ · ‖A‖ · ‖∆b‖‖A‖ · ‖x‖

≤

≤ ‖A−1‖ · ‖A‖(‖∆A‖‖A‖

+‖∆b‖‖b‖

),

ceea ce demonstreaza (2.59). (La ultima inegalitate, am tinut seama ca ‖b‖ =‖Ax‖ ≤ ‖A‖ · ‖x‖.)

Asadar numarul de conditionare κ(A) definit ın (2.60) furnizeaza o masuraa conditionarii problemei rezolvarii unui sistem liniar Ax = b. Un κ(A) miccorespunde unei bune conditionari, adica unei sensibilitati reduse a solutieiın raport cu perturbatiile datelor. Se observa ca doar matricea A determinaconditionarea unui sistem, termenul liber b neavand nici un rol; de aceeavorbim de ”conditionarea matricei A” ın loc de ”conditionarea sistemului”.Vom atasa un indice notatiei κ(A) atunci cand vom preciza norma folositaın (2.60); de exemplu, κ∞(A) = ‖A−1‖∞ · ‖A‖∞.

Este intuitiv faptul ca un sistem este cu atat mai rau conditionat cu catmatricea A este mai ”aproape” de o matrice singulara. Trebuie subliniatca masura acestei apropieri este data de numarul de conditionare si nu dedeterminantul matricei A.

Exemplul 2.26 Sa consideram matricele D,T ∈ Rn×n

D =

0.1

0.1. . .

0.1

, T =

1 −1 . . . −1

1 . . . −1. . .

...1

.

2.11 Metode iterative pentru rezolvarea sistemelor liniare 59

Fiind diagonala, este evident ca D este bine conditionata; se verifica imediatca κ(D) = 1, ın schimb det(D) = 10−n, adica, daca de exemplu n = 20, ovaloare foarte mica. Pentru matricea triunghiulara T situatia este contrara:desi det(A) = 1, matricea este rau conditionata pentru ca κ∞(T ) = n2n−1

(U = T−1 este inferior triunghiulara unitate cu elementele uij = 2j−i−1,pentru j > i, deci ‖T−1‖∞ = 1 + 2 + 22 + . . .+ 2n−2 = 2n−1 − 1. Pe de altaparte, ‖T‖∞ = n).

2.11 Metode iterative pentru rezolvarea sistemelorliniare

Metodele iterative pentru rezolvarea sistemului (2.2) construiesc un sir devectori x(k) ∈ Rn care tinde (cand k → ∞) la solutia x a sistemului. Prindefinitie, convergenta acestui sir de vectori este echivalenta cu convergenta(scalara) ‖x(k)−x‖ → 0. Convergenta unui sir de matrice se defineste similar,norma implicata fiind evident una matriceala.

Incepem cu cateva detalii tehnice privind matricele convergente.

Definitia 2.27 O matrice T ∈ Rn×n este convergenta daca

T k → 0. (2.61)

In sens larg, matricele convergente sunt generalizarea multidimensionalaa scalarilor subunitari. Ele pot fi caracterizate de proprietati ale valorilorproprii. Raza spectrala a unei matrice este egala cu cel mai mare modul alunei valori proprii.

Propozitia 2.28 O matrice este convergenta daca si numai daca raza eispectrala este subunitara.

Asa cum vom vedea ın capitolul 4, calculul valorilor proprii este mai com-plicat decat rezolvarea sistemelor liniare, deci conditia de convergenta bazatape raza spectrala este mai greu de folosit. De aceea pot fi utile conditii sufi-ciente de genul celei care urmeaza.

Propozitia 2.29 O matrice T ∈ Rn×n este convergenta daca ‖T‖ < 1, unde‖ · ‖ este o norma matriceala consistenta.

Demonstratie. Definitia (2.53) a consistentei spune ca ‖T k‖ ≤ ‖T‖k, deciT k → 0 atunci cand ‖T‖ < 1.

Descriem ın continuare principiul unei singure categorii (poate cea maisimpla dintre multele existente) de metode iterative pentru rezolvarea sis-temelor liniare. Metoda este descrisa de iteratia


x(k+1) = Tx(k) + c, k = 1, 2, . . . (2.62)

ın care fiecare vector din sir este construit afin din precedentul, iar vectorulinitial x(0) ∈ Rn este dat. Matricea T ∈ Rn×n si vectorul c ∈ Rn sunt obtinutidin A si b prin procedee pe care le vom ilustra ulterior. Deocamdata studiemconvergenta sirului construit de (2.62).

Teorema 2.30 Sirul (2.62) este convergent pentru orice initializare x(0)

daca si numai daca matricea T este convergenta.

Demonstratie. Presupunem ıntai ca x(k) → x. Trecand la limita ın ambeleparti ale (2.62), rezulta ca x = Tx+ c, deci x = (I − T )−1c. (Observam decica are loc egalitatea (I − T )−1c = A−1b, ceea ce limiteaza alegerea lui T sic.) Scazand x = Tx+ c din (2.62) obtinem

x(k+1) − x = T (x(k) − x),

ceea ce implicax(k) − x = T k(x(0) − x),

deciT k(x(0) − x)→ 0, ∀x(0) ∈ Rn.

Avand ın vedere ca vectorul care ınmulteste T k este arbitrar, singura posi-bilitate ca relatia de mai sus sa fie adevarata este ca T k → 0, deci T esteconvergenta.

Reciproc, presupunem acum ca T este convergenta. Substituind repetat ın(2.62) obtinem

x(1) = Tx(0) + cx(2) = Tx(1) + c = T 2x(0) + (T + I)c

...x(k) = T kx(0) + (T k−1 + . . .+ T + I)c

Deoarece T k → 0, observand ca

(I − T )(T k−1 + . . .+ T + I) = I − T k,

rezulta cax(k) → (I − T )−1c.

Deci sirul x(k) este convergent.

Prezentam mai departe doua modalitati simple de a obtine, pe baza rezul-tatului de mai sus, algoritmi iterativi de rezolvare a sistemelor liniare. Ele sebazeaza pe scrierea

A = L+D + U, (2.63)


unde L este strict inferior triunghiulara, D este diagonala si U este strictsuperior triunghiulara. (Nu este nevoie de nici un calcul pentru a obtine(2.63), ci se extrag pur si simplu elementele respective din A.)

Metoda iterativa Jacobi. Presupunem ca matricea diagonala D estenesingulara. Din simpla observatie ca

Ax = b ⇔ Dx = −(L+ U)x+ b ⇔ x = −D−1(L+ U)x+D−1b

putem construi iteratia (numita metoda Jacobi)

x(k+1) = −D−1(L+ U)x(k) +D−1b. (2.64)

Prin modul de constructie suntem siguri ca, daca metoda converge, atunci eaconverge la solutia sistemului Ax = b. La nivel de element, relatia (2.64) areforma

x(k+1)i =

bi − n∑j=1,j 6=i

aijx(k)j

/aii, i = 1 : n (2.65)

Observam ca relatia (2.64) este un caz particular al formei generale (2.62)cu

T = −D−1(L+ U), c = D−1b. (2.66)

Pentru implementare, observam ca sunt suficienti doi vectori pentru calcu-lul aproximatiilor succesive ale solutiei, corespunzand lui x(k+1) si x(k) (ıi vomnota x, respectiv y, ın implementarea de mai jos). O problema tipica pentrualgoritmii iterativi este stabilirea numarului de iteratii efective; convergentaeste un proces asimptotic care trebuie oprit dupa un numar finit (si nu preamare) de iteratii. Necunoscand limita sirului x(k), care este chiar solutia sis-temului pe care vrem sa ıl rezolvam, decizia de oprire poate fi luata canddiferenta (relativa) dintre doi termeni succesivi ai sirului,

‖x(k+1) − x(k)‖‖x(k+1)‖

, (2.67)

este suficient de mica. In plus, trebuie sa fixam si un numar maxim deiteratii, Nmax, pentru eventualitatea ın care sirul nu converge. Evitam astfelramanerea ıntr-o bucla infinita. Tinand seama de toate acestea, algoritmulare forma urmatoare.

Algoritmul 2.17 (JACOBI – Metoda iterativa Jacobi pentru cal-culul solutiei unui sistem liniar) (Se dau A ∈ Rn×n nesingulara sib ∈ Rn; aproximatia initiala x(0) ∈ Rn; toleranta ε pentru oprire;numarul maxim de iteratii Nmax. Se calculeaza o aproximatie a solutieix a sistemului Ax = b.)

1. y ← x(0)


2. Pentru k = 1 : Nmax

1. Pentru i = 1 : n

1. xi ←(bi −

∑nj=1,j 6=i aijyj

)/aii

2. Daca ‖x− y‖/‖x‖ < ε atunci Stop

3. y ← x3. Mesaj: ”a fost atins numarul maxim de iteratii”

Comentarii. Toleranta ε este mijlocul prin care utilizatorul decide acurateteasolutiei calculate; valorile tipice pot merge de la 10−4 la 10−8 sau chiar maimici. Desigur, o valoare mai mica ınseamna mai multe iteratii executate.Numarul Nmax poate fi de ordinul sutelor, miilor, sau chiar mai mare, ınfunctie de dimensiunea matricei. El trebuie sa fie suficient de mare pentru anu ımpiedica functionarea normala a algoritmului, dar nu atat de mare ıncatsa prelungeasca inutil executia ın caz de neconvergenta sau convergenta foartelenta.

Numarul de operatii este O(n2) pe iteratie, dominat de ınmultirea matrice-vector din (2.64). Pentru a fi comparabil cu GPP, algoritmul ar trebui sa setermine ın mai putin de n iteratii, ceea ce se ıntampla doar ın putine cazuri.Algoritmul Jacobi este indicat ın special pentru matrice rare; acestea auputine elemente nenule; memorandu-le doar pe acestea, ınmultirea matrice-vector se poate implementa foarte eficient.

In mod evident, algoritmul Jacobi nu converge pentru orice matrice A. Maimult, daca A are un element nul pe diagonala, atunci nu poate fi executatanici macar o iteratie. Exista totusi o categorie de matrice (care apare si ınunele aplicatii) pentru care convergenta este asigurata.

Definitia 2.31 O matrice A ∈ Rn×n este diagonal dominanta (pe linii) dacaelementele diagonale sunt mai mari (ın modul) decat suma modulelor celor-lalte elemente de pe linia respectiva:

|aii| >n∑

j=1,j 6=i

|aij |, i = 1 : n (2.68)

(O definitie similara poate fi data pe coloane.)

Propozitia 2.32 Daca matricea A este diagonal dominanta, atunci algorit-mul Jacobi converge.

Demonstratie. Aratam ca matricea T din (2.66) satisface ‖T‖∞ < 1. In acestcaz, conform Propozitiei 2.29, sirul produs de algoritmul Jacobi converge.

Din (2.68) rezulta

n∑j=1,j 6=i

|aij ||aii|

< 1, ∀i = 1 : n,


de unde obtinem

‖T‖∞ = ‖ −D−1(L+ U)‖∞ = maxi=1:n

n∑j=1,j 6=i

| − aij ||aii|

< 1,

ceea ce ıncheie demonstratia.

Desigur, algoritmul Jacobi converge si pentru alte matrice decat cele diag-onal dominante, dar pentru acestea nu exista un mod simplu de caracterizare.

Algoritmul Gauss-Seidel. In relatia (2.65) care defineste algoritmul Ja-cobi, vectorul ”urmator” x(k+1) este construit folosind exclusiv elementelevectorului curent x(k). Deoarece este de asteptat ca x(k+1) sa fie o aproximatiemai buna a solutiei, putem folosi elementele lui, imediat dupa ce au fost cal-culate, ın locul celor ale lui x(k). Presupunand ca parcurgem elementele dela 1 la n, relatia (2.65) se transforma ın

x(k+1)i =

bi − i−1∑j=1

aijx(k+1)j −

n∑j=i+1

aijx(k)j

/aii, i = 1 : n. (2.69)

Deci, atunci cand calculam x(k+1)i , folosim elementele x

(k+1)j pentru j < i,

acestea fiind deja calculate, si elementele x(k)j pentru j > i. Relatia (2.69)

defineste algoritmul Gauss-Seidel, care foloseste mereu cea mai noua versiunea unui element al aproximarii solutiei.

Algoritmul Gauss-Seidel este similar algoritmului Jacobi 2.17, singuradiferenta fiind ca instructiunea 2.1.1 este ınlocuita de

xi ←

bi − i−1∑j=1

aijxj −n∑

j=i+1

aijyj

/aii.

Se poate demonstra usor ca iteratia Gauss-Seidel (2.69) este un caz par-ticular al formei generale (2.62) cu

T = −(L+D)−1U, c = (L+D)−1b. (2.70)

Observatia 2.33 Fara demonstratie, precizam ca, la fel ca metoda Jacobi,metoda Gauss-Seidel converge pentru matrice diagonal dominante. De aseme-nea, ea converge pentru matrice simetrice si pozitiv definite.

In general, metoda Gauss-Seidel este considerata superioara metodei Ja-cobi, desi convergenta nu este neaparat mai rapida pentru orice matrice (pen-tru care ambele metode converg). Desigur, masuram convergenta ın numarulde iteratii necesar atingerii unei erori fixate; numarul de operatii pe iteratieeste acelasi pentru cele doua metode.

Capitolul 3

Problema liniara a celor mai micipatrate

In acest capitol discutam despre rezolvarea sistemului Ax = b, cu A ∈ Rm×navand rang maxim, anume rang(A) = min(m,n). Consideram doar sistemecu m 6= n, cele determinate fiind subiectul capitolului anterior. Sunt douacazuri importante.

Sisteme supradeterminate. Cand m > n (sistemul are mai multeecuatii decat necunoscute), sistemul nu are solutii, ın general. Reamintimca imaginea unei matrice se defineste prin

ImA = {y ∈ Rm | ∃x ∈ Rn astfel ıncat y = Ax} (3.1)

si reprezinta subspatiul generat de coloanele lui A. Dimensiunea lui ImA estedeci n, mai mica decat cea a spatiului ın care se afla b, anume Rm. Sistemulare solutie doar daca b ∈ ImA. In lipsa unei solutii exacte, este natural sa oınlocuim cu acel vector x care face ca Ax sa fie cel mai aproape de b.

Definitia 3.1 Solutia ın sensul celor mai mici patrate (CMMP) a sistemuluisupradeterminat Ax = b este acel vector x∗ pentru care

‖b−Ax∗‖2 = minx∈Rn

‖b−Ax‖2. (3.2)

Vectorul b − Ax este numit reziduu. Solutia CMMP este deci cea careminimizeaza norma euclidiana a reziduului. Numele CMMP vine de la simplulfapt ca norma euclidiana se defineste printr-o suma de patrate.

Exemplul 3.2 Problemele CMMP apar natural atunci cand se doreste cal-cularea parametrilor unei legi liniare. Consideram un exemplu simplu carepoate fi usor ilustrat grafic. Presupunem ca doua marimi fizice, ξ si η, suntlegate prin relatia

η = αξ + β, (3.3)

adica η este o functie afina de ξ. Coeficientii α si β sunt necunoscuti si trebuiedeterminati din masuratori. Dispunem de m masuratori (ξi, ηi), i = 1 : m.

65

66 3 Problema liniara a celor mai mici patrate

×ξ1

η1

×

ξ2

η2×

ξ3

η3

×

ξm−1

ηm−1

×

ξm

ηm

. . .

...

× ×

× ×

Fig. 3.1 Solutie CMMP pentru o problema de aproximare liniara.

Problema este ilustrata ın Figura 3.1. Relatia ideala (3.3) este reprezentata deo dreapta. Din cauza erorilor de masurare, nu obtinem puncte pe dreapta ciın apropierea ei (ın figura distantele sunt exagerat de mari, pentru vizualizaremai buna). Putem calcula coeficientii din (3.3) minimizand norma erorii

m∑i=1

(ηi − αξi − β)2 =

∥∥∥∥∥∥∥∥∥

η1η2...ηm

−ξ1 1ξ2 1...

...ξm 1

[αβ

]∥∥∥∥∥∥∥∥∥2

= ‖b−Ax‖2. (3.4)

(Constructia matriceiA si a vectorului b ce definesc problema CMMP din dateeste evidenta.) Aceasta reprezinta suma patratelor distantelor de la punctela dreapta masurate pe verticala (cu rosu ın figura). Cu alte cuvinte, se pre-supune ca valorile ξi sunt masurate perfect si doar cele ηi contin erori. Aceastaeste de multe ori valabil ın practica: ξi este o intrare cu valoare usor de fixatsi doar iesirea ηi se masoara cu posibile erori. Desigur, cu cat luam mai multepuncte, cu atat este mai probabil ca erorile sa se compenseze si sa obtinemvalori ale parametrilor necunoscuti α si β mai aproape de cele adevarate.

Calculul solutiei CMMP (3.2) este subiectul principal al acestui capitol.

Sisteme subdeterminate. Cand m > n (mai multe necunoscute decatecuatii) sistemul Ax = b are o infinitate de solutii. In general, mai bine decata lua o solutie particulara oarecare, este preferabil sa o alegem pe cea mai”mica”.

Definitia 3.3 Solutia CMMP a sistemului subdeterminat Ax = b este aceasolutie x∗ pentru care

‖x∗‖2 = minAx=b

‖x‖2. (3.5)

Solutia problemei (3.5) va fi discutata ın sectiunea 3.7. Problema CMMPsubdeterminata este mai putin importanta decat cea supradeterminata,deoarece apare mai rar ın practica, si se rezolva cu instrumente asemanatoare.

3.1 Solutia normala a problemei CMMP 67

ImA

Ax∗

b

b−Ax∗

Fig. 3.2 Geometria problemei CMMP.

3.1 Solutia normala a problemei CMMP

Pentru a rezolva problema CMMP supradeterminata, analizam ıntai geome-tria acesteia, ilustrata ın Figura 3.2. Asa cum am vazut mai sus, ImA esteun subspatiu ın spatiul Rm (reprezentat printr-un plan ın figura). Vectorulb nu se afla, ın general, ın acest subspatiu. Reziduul b − Ax este un vectorcare porneste din subspatiul ImA si are varful comun cu b. Rezulta imediatca minimizarea normei reziduului b − Ax, impusa ın (3.2), este echivalentacu ortogonalitatea reziduului pe ImA. Acest fapt duce imediat la o formulapentru calculul solutiei. Inainte de a o prezenta, reamintim definitia unui altsubspatiu interesant, anume nucleul matricei A:

KerA = {x ∈ Rn | Ax = 0}. (3.6)

Nucleul contine vectorii avand ca elemente coeficientii combinatiilor liniarenule ale coloanelor lui A. Reamintim de asemenea ca subspatiile ImA siKerAT sunt ortogonale (uT v = 0, ∀u ∈ ImA, ∀v ∈ KerAT ) si complementare(ın sensul ca ImA + KerAT = Rm: orice vector w ∈ Rm se poate scrie subforma w = u+ v, cu u ∈ ImA si v ∈ KerAT ).

Teorema 3.4 Atunci cand matricea A ∈ Rm×n are rang maxim (si m > n),solutia problemei CMMP supradeterminate (3.2) este

x∗ = (ATA)−1AT b. (3.7)

Demonstratie. Am vazut mai sus ca b−Ax∗ ⊥ ImA, ceea ce implica b−Ax∗ ∈KerAT . Conform (3.6), aceasta ınseamna

AT (b−Ax∗) = 0,

ceea ce duce imediat la (3.7).

Observatia 3.5 Sistemul liniar ATAx = AT b a carui solutie este (3.7) se mainumeste si sistemul ecuatiilor normale asociat problemei CMMP. Observam


ca matricea ATA este simetrica si pozitiv definita. Asadar, solutia CMMP sepoate calcula cu urmatorul algoritm

1. Calculeaza B = ATA2. Calculeaza factorizarea Cholesky B = LLT

3. Rezolva sistemul inferior triunghiular Ly = b4. Rezolva sistemul superior triunghiular LTx = y

Acest mod de calcul este recomandabil doar pentru matrice A cu dimen-siuni (foarte) mici. Desi am prezentat notiunea de numar de conditionaredoar pentru matrice patrate, aceasta se poate extinde si la matrice drep-tunghiulare, cu notiuni prezentate ın capitolul 5. Fara a discuta definitia,remarcam faptul ca κ2(ATA) = κ2(A)2. Asadar, sistemul ecuatiilor nor-male are un numar de conditionare mult mai mare decat κ(A) (de exempluκ2(A) = 105 ⇒ κ2(ATA) = 1010), ceea ce ıl poate face usor rau conditionat.Este deci preferabil sa lucram direct cu matricea A.

Ideea de rezolvare a problemei CMMP este similara celei de la elim-inarea gaussiana: transformarea matricei A ıntr-una superior triunghiularaprin transformari elementare. Pentru a vedea ce transformari sunt posibile,ne mai reamintim cateva notiuni de algebra liniara.

Definitia 3.6 O matrice Q ∈ Rm×m este ortogonala daca QTQ = QQT = I.

Asadar, inversa unei matrice ortogonale este transpusa ei. Coloanele uneimatrice ortogonale sunt vectori ortogonali si au norma egala cu 1.

Propozitia 3.7 (Norma euclidiana a unui vector este pastrata la ınmultireacu o matrice ortogonala.) Daca Q ∈ Rm×m este o matrice ortogonala siy ∈ Rm este un vector, atunci

‖Qy‖2 = ‖y‖2. (3.8)

Demonstratie. ‖Qy‖2 = (Qy)T (Qy) = yTQTQy = yT y = ‖y‖2.

Putem deci transforma reziduul b−Ax din (3.2) prin ınmultire cu o matriceortogonala, pastrand echivalenta cu problema initiala datorita conservariinormei reziduului. Sectiunea urmatoare prezinta transformari ortogonale el-ementare, care conduc la algoritmi cu complexitate redusa.

3.2 Reflectori Householder

Reflectorii elementari sunt omologul ortogonal al matricelor inferior triunghi-ular elementare din capitolul precedent, ın sensul ca au un scop similar –introducerea de zerouri ıntr-un vector.

3.2 Reflectori Householder 69

Definitia 3.8 Un reflector elementar (Householder) de ordin m si indice keste o matrice Uk ∈ Rm×m de forma

Uk = Im −1

βkuku

Tk (3.9)

undeuk = [0 . . . 0 ukk . . . umk]T (3.10)

este un vector avand primele k − 1 elemente nenule si

βk = ‖uk‖2/2.

Un reflector are structura

Uk =

[Ik−1 0

0 Uk

],

unde Uk ∈ R(m−k+1)×(m−k+1) este un reflector elementar de ordin m− k+ 1si indice 1. Principalele proprietati ale unui reflector sunt urmatoarele. (Incontinuare, renuntam la calificativul ”elementar”.)

Propozitia 3.9 a) Un reflector Uk este o matrice ortogonala si simetrica.

b) Fie x ∈ Rm un vector dat, iar Uk un reflector dat. Atunci vectorultransformat y = Ukx are elementele

(Ukx)i =

{xi, pentru i = 1 : k − 1,

xi − τuik, pentru i = k : m.(3.11)

unde

τ =1

βk

m∑i=k

uikxi. (3.12)

c) Fie x ∈ Rm un vector dat ale carui elemente xi, i = k+ 1 : m, nu sunttoate nule. Atunci exista un reflector Uk astfel ıncat vectorul transformaty = Ukx are ultimele m− k elemente nule (yi = 0, i = k + 1 : m). Notand

σ2 =

m∑i=k

x2i 6= 0, (3.13)

reflectorul Uk cautat este definit de

uik =

0, pentru i = 1 : k − 1,

xk + σ, pentru i = k,

xi, pentru i = k + 1 : m.

(3.14)

Mai mult, vectorul rezultat y are forma


(Ukx)i =

xi, pentru i = 1 : k − 1,

−σ, pentru i = k,

0, pentru i = k + 1 : m.

(3.15)

Demonstratie. a) Din (3.9) este clar ca UTk = Uk, deci reflectorul estesimetric. De asemenea,

UTk Uk =

(I − 1

βkuku

Tk

)(I − 1

βkuku

Tk

)= I− 2

βkuku

Tk +

1

β2k

uk(uTk uk)uTk = I,

deoarece uTk uk = ‖uk‖2 = 2βk. Asadar un reflector elementar este propria luiinversa.

b) Calculam

y = Ukx =

(I − 1

βkuku

Tk

)x = x− 1

βkuku

Tk x = x− uTs x

βkuk

(3.12)= x− τuk.

La nivel de element, aceasta relatie este identica cu (3.11).

c) Din (3.13) si (3.14) rezulta ca

uTk x =

m∑i=1

uikxi = xk(xk + σ) +

m∑i=k+1

x2i = σ2 + xkσ (3.16)

si

βk =1

2

m∑i=k

u2ik =1

2

((xk + σ)2 +

m∑i=k+1

x2i

)

=1

2

(x2k + 2xkσ + σ2 +

m∑i=k+1

x2i

)= σ2 + xkσ = σ(xk + σ).

(3.17)

Deci, din (3.14) si (3.16) se obtine

yi = xi − uis =

xi, pentru i = 1 : k − 1

−σ, pentru i = k

0, pentru i = s+ 1 : n

adica relatia (3.15) este demonstrata.

Propozitia 3.9 da un algoritm aproape complet pentru calculul unui reflec-tor care anuleaza elemente ale unui vector. Constanta σ definita ın (3.13) nuare semnul determinat de aceasta relatie. Pentru a-l alege, remarcam ca ıncalcularea valorilor ukk si βs, ın (3.17 si (3.14), apare calculul sumei xk + σ.Pentru a evita o posibila anulare catastrofala (revedeti sectiunea 1.2) alegemca σ sa aiba acelasi semn ca xs. Obtinem astfel urmatorul algoritm.

3.3 Triangularizarea ortogonala 71

Algoritmul 3.1 (Se dau un vector x ∈ Rm si un indice k ∈ 1 : m−1.Se calculeaza reflectorul elementar Uk astfel ıncat (Ukx)i = 0, i = k+1 :m. Se calculeaza vectorul Ukx pe loc ın x.)

1. σ ← sign(xk) ·√∑m

i=k x2i

2. ukk = xk + σ3. uik = xi, pentru i = k + 1 : m4. βk = σukk5. xk ← −σ

Comentarii. Pentru claritate, am memorat separat elementele definitorii alereflectorului, anume vectorul uk (elementele k : m) si scalarul βk. Se poateface o memorare mai economica, folosind elementele anulate ale vectorului x,dar aceasta ar complica scrierea algoritmilor urmatori; ın practica se folosestememorarea economica. Evident, nu este necesara memorarea explicita a val-orilor nule ale vectorului x transformat, ci doar cea a singurului elementmodificat, xk.

3.3 Triangularizarea ortogonala

Asa cum am vazut la sfarsitul sectiunii 3.1, rezolvarea problemei CMMPse poate face prin transformarea matricei A cu transformari ortogonale. Inaceasta sectiune vom arata cum se poate aduce o matrice A ∈ Rm×n, cum > n la forma superior triunghiulara cu ajutorul reflectorilor elementari.Vom vedea apoi ca rezolvarea problemei CMMP devine mult mai simpla candmatricea este superior triunghiulara. Urmatoarea teorema, similara ın spiritcu Teorema 2.7, descrie procesul de triangularizare.

Teorema 3.10 Fie A ∈ Rm×n, cu m > n. Atunci exista reflectorii Uk,k = 1 : n, astfel ıncat matricea

Un . . . U2U1A = R (3.18)

este superior triunghiulara.

Demonstratie. Planul demonstratiei este similar celui de la Teorema 2.7.Pentru simplitate, presupunem direct ca transformarile se efectueaza ın A,descriind astfel algoritmul de triangularizare. Lucram ın mod natural cupartitionarea pe coloane

A = [a1 a2 . . . an]


Pasul 1. Folosind Algoritmul 3.1, calculam reflectorul U1 care introducezerouri subdiagonale ın prima coloana a matricei A, adica (U1a1)i = 0, i =2 : m. Apoi calculam

A← U1A = [U1a1 U1a2 . . . U1an],

produsele U1aj , j = 2 : n, putand fi calculate cu formulele din Propozitia3.9b (ın care k = 1 si locul vectorului x este luat de coloana aj). (Desigur,daca prima coloana este deja nula sub diagonala, reflectorul este U1 = I.)

Pasul k. Presupunem ca matricea curenta este superior triunghiulara ınprimele k − 1 coloane, adica are forma

A← Uk−1 . . . U2U1A =

a11 a12 . . . a1k . . . a1na22 . . . a2k . . . a2n

. . . . . .akk . . . akn

0 ak+1,k . . . ak+1,n

. . .amk . . . amn

.

Folosim Algoritmul 3.1 pentru a calcula reflectorul Uk astfel ıncat

(Ukak)i = 0, i = k + 1 : m,

obtinand astfel zerouri subdiagonale pe coloana k. Calculand produsul

A← UkA = [Uka1 . . . Ukak . . . Ukan],

observam ca produsele Ukaj se ımpart ın doua categorii. Cand j < k(coloanele din stanga celei curente), deoarece aij = 0 pentru i = k : m,rezulta ca uTk aj = 0, deci Ukaj = aj (vezi Propozitia 3.9b, unde pentrux = aj se obtine τ = 0); asadar zerourile subdiagonale introduse ın pasii1 : k − 1 nu sunt alterate. Pentru coloanele j = k + 1 : n, calculul Ukaj seface cu formulele din Propozitia 3.9b. Pe scurt, matricea A ← UkA este su-perior triunghiulara ın primele k coloane si avem un algoritm pentru calcululmatricei.

Deci procedura de calcul poate continua pana la pasul n (de care e nevoiepentru ca si ultima coloana are zerouri subdiagonale) si astfel relatia (3.18)este satisfacuta.

Observatia 3.11 Demonstratia de mai sus se poate utiliza si pentru cazulın care m = n, ın care sunt ınsa n− 1 pasi (si reflectori) ın loc de n.

Demonstratia de mai sus furnizeaza o procedura de triangularizare ortog-onala a unei matrice, organizata conform schemei

3.4 Rezolvarea problemei CMMP 73

1. Pentru k = 1 : n1. Calculeaza reflectorul Uk (adica vectorul uk si scalarul βk),

astfel ıncat (UkA)ik = 0, pentru i = k + 1 : m2. Calculeaza A← UkA

Renuntand din motive didactice la optimizarea memoriei si folosind rezul-tatele din Propozitia 3.9, obtinem algoritmul urmator.

Algoritmul 3.2 (TORT – Triangularizare ortogonala) (Se da A ∈Rm×n, cu m > n. Se calculeaza reflectorii Uk si matricea superior tri-unghiulara R astfel ıncat are loc relatia (3.18). Matricea R se memo-reaza peste triunghiul superior al lui A. Vectorii uk ce definesc reflectoriiUk, k = 1 : n, se memoreaza ıntr-un tablou U ; similar, scalarii βk sememoreaza ıntr-un vector β.)

1. Pentru k = 1 : n

1. σ ← sign(akk) ·√∑m

i=k a2ik

2. ukk = akk + σ3. Pentru i = k + 1 : m

1. uik = aik4. βk = σukk5. akk ← −σ6. Pentru j = k + 1 : n

1. τ ← (∑mi=k uikaij) /βk

2. Pentru i = k : m1. aij ← aij − τuik

Comentarii. Complexitatea algoritmului este NTORT = 2(mn2−n3/3) flopi.Deci, ın cazul m = n, algoritmul TORT este de doua ori mai complex decatalgoritmul G.

Se observa ca partea subdiagonala a matricei A ramane nefolosita. In eaar putea fi memorate valorile uik, cu i > k; ın acest fel ar mai fi necesar doarun vector pentru a memora separat valorile ukk, k = 1 : n. Cititorul esteinvitat sa modifice algoritmul pentru a optimiza astfel memoria. Necesarulefectiv de memorie este MTORT = mn+ 2n.

3.4 Rezolvarea problemei CMMP

Folosim relatia (3.18) ce defineste triangularizarea ortogonala pentru trans-formarea reziduului ın scopul rezolvarii problemei CMMP (3.2). Notand

d = Un . . . U2U1b (3.19)


si tinand seama de ortogonalitatea reflectorilor si de Propozitia 3.7, observamca

‖b−Ax‖ = ‖Un . . . U2U1(b−Ax)‖ = ‖d−Rx‖.

Evidentiind primele n linii, scriem

d =

[d′

d′′

], R =

[R′

0

],

cu d′ ∈ Rn, d′′ ∈ Rm−n, R′ ∈ Rn×n. Utilizand aceasta scriere si teorema luiPitagora, obtinem

‖b−Ax‖2 =

∥∥∥∥[d′ −R′xd′′

]∥∥∥∥2 = ‖d′ −R′x‖2 + ‖d′′‖2.

Al doilea termen nu depinde de x (este deci fix) iar primul este minim atuncicand d′ −R′x = 0. Egalitatea este posibila deoarece, atunci cand A are rangmaxim (si deci si R), rezulta ca matricea R′ este nesingulara. Asadar solutiaCMMP se poate calcula rezolvand sistemul superior triunghiular R′x = d′.Obtinem algoritmul urmator.

Algoritmul 3.3 (CMMP – Calculul solutiei CMMP a unui sistemsupradeterminat) (Se dau A ∈ Rm×n, m > n, cu rang maxim, si b ∈ Rn.Se calculeaza solutia CMMP a sistemului supradeterminat Ax = b.)

1. [U, β,R] = TORT(A)2. Pentru k = 1 : n

1. τ ← (∑mi=k uikbi) /βk

2. Pentru i = k : m1. bi ← bi − τuik

3. x = UTRIS(R, b)

Comentarii. Instructiunea 2 este implementarea pe loc a relatiei (3.19), ınfiecare iteratie efectuandu-se operatia b← Ukb folosind Propozitia 3.9b.

Efortul de calcul este dominat de triangularizarea ortogonala, restul opera-tiilor fiind O(mn) (pasul 2) sau O(n2) (pasul 3).

3.5 Factorizarea QR

Definitia 3.12 Se numeste factorizare QR a matricei A ∈ Rm×n scriereaacesteia ın forma

A = QR, (3.20)

unde Q ∈ Rm×m este ortogonala si R ∈ Rm×n este superior triunghiulara.

3.6 Triangularizare ortogonala cu rotatii 75

Observam ca egalitatea (3.20) poate fi scrisa si ın forma

A = QR = [Q′ Q′′]

[R′

0

]= Q′R′, (3.21)

cu Q′ ∈ Rm×n avand coloanele ortogonale si R′ ∈ Rn×n. Si aceasta relatieeste numita tot factorizare QR (sau uneori forma comprimata a factorizariiQR).

Calculul factorizarii QR se poate face folosind rezultatul algoritmuluiTORT. Din (3.18), stiind ca reflectorii sunt simetrici si ortogonali, rezultaca

A = U1U2 . . . UnR,

deciQ = U1U2 . . . Un. (3.22)

Asadar R rezulta din TORT iar Q poate fi calculat ulterior prin schema

1. Q← I2. Pentru k = n : −1 : 1

1. Q← UkQ

Implementarea este lasata cititorului ca exercitiu. In majoritatea cazurilor,calculul explicit al matriceiQ nu este necesar; este mai eficienta utilizarea pro-dusului (3.22), adica a secventei de reflectori rezultati din algoritmul TORT.

Observatia 3.13 Factorizarea QR este utila pentru calcularea unor bazeortogonale asociate cu subspatii generate de coloanele lui A.

Presupunand (ca si pana acum) ca rangA = n, putem demonstra cacoloanele matricei Q′ din (3.21) formeaza o baza ortogonala pentru ImA.Intr-adevar, matricea R′ este nesingulara si atunci Ax = Q′R′x = Q′z. Deci(revedeti Definitia (3.1)), pentru un y ∈ ImA exista (prin definitie) un x ∈ Rnastfel ıncat y = Ax, dar si z = (R′)−1x astfel ıncat y = Q′z. Rezulta caImA = ImQ′, ceea ce justifica afirmatia initiala.

Desigur, avand ın vedere ca subspatiile ImA si KerAT sunt ortogonalesi ImA ∪ KerAT = Rm, coloanele matricei Q′′ din (3.21) formeaza o bazaortogonala pentru KerAT .

3.6 Triangularizare ortogonala cu rotatii

La fel ca reflectorii, rotatiile Givens sunt matrice ortogonale elementarecare pot fi utilizate pentru triangularizare ortogonala, rezolvarea problemeiCMMP si, ın general, pentru orice probleme ale caror solutii pot fi aflate cutransformari ortogonale.


Definitia 3.14 O matrice Gij ∈ Rm×m se numeste rotatie Givens cu indiciii si j (prin conventie i < j) daca are forma

Gij =

Ii−1

c sIj−i−1

−s cIm−j

, (3.23)

unde scalarii c si s satisfac relatia

c2 + s2 = 1. (3.24)

O rotatie (vom omite mai departe numele lui Givens, care este cel carea propus utilizarea rotatiilor ın calculul numeric) este, de fapt, o matrice2 × 2 definita de elementele c si s, scufundata ın matricea identitate m ×m. Cazuri particulare de rotatii sunt matricea identitate (c = 1, s = 0) sipermutarile elementare (c = 0, s = 1), vezi Definitia 2.4. Similar cu Propozitia3.9, rezultatele de mai jos descriu proprietatile si modul de lucru cu rotatii.

Propozitia 3.15 a) O rotatie Gij este o matrice ortogonala.

b) Fie x ∈ Rm un vector dat, iar Gij o rotatie data. Atunci vectorultransformat y = Gijx are elementele

y` =

x`, pentru ` 6= i, j,

cxi + sxj , pentru ` = i,

−sxi + cxj , pentru ` = j.

(3.25)

c) Fie x ∈ Rm un vector dat. Atunci exista o rotatie Gij astfel ıncatvectorul transformat y = Gijx satisface yj = 0. Notand

r =√x2i + x2j ,

rotatia este definita de

c =xir, s =

xjr

(3.26)

iar vectorul transformat este

y` =

x`, pentru ` 6= i, j,

r, pentru ` = i,

0, pentru ` = j.

(3.27)

3.6 Triangularizare ortogonala cu rotatii 77

Demonstratie. a) Prin calcul direct se constata ca GTijGij = I (se poatelucra cu matrice 2× 2). Atentie, o rotatie nu e simetrica.

b) Prin structura ei, este evident ca o rotatie nu poate modifica decat douaelemente ale unui vector pe care ıl ınmulteste. Relatia (3.25) rezulta imediatprin simpla ınmultire.

c) Se verifica imediat introducand (3.26) ın (3.25).

Observatia 3.16 Pentru a anula mai multe elemente ale unui vector, estenevoie de mai multe rotatii, cate una pentru fiecare element anulat. De ex-emplu, pentru a anula elementele k + 1 : m ale unui vector x ∈ Rm, putemcalcula

Gkn . . . Gk,k+2Gk,k+1x,

sauGk,k+1Gk+1,k+2 . . . Gm−1,mx.

Sunt posibile si alte secvente de rotatii. Atentie la ordine, ca nu cumva orotatie sa strice un zero introdus de o rotatie precedenta. De exemplu, ınprimul exemplu de mai sus, ordinea inversa este si ea corecta, dar ın al doileanu.

Folosind prima dintre cele doua secvente de rotatii de mai sus, obtinemun analog cu rotatii pentru Algoritmul 3.1.

Algoritmul 3.4 (Se dau un vector x ∈ Rm si un indice k ∈ 1 : m−1.Se calculeaza rotatiile Gki, i = k + 1 : m, astfel ıncat vectorul y =Gkn . . . Gk,k+2Gk,k+1x satisface yi = 0, i = k + 1 : m. Rotatia Gki sememoreaza prin elementele ci, si, iar vectorul y se memoreaza peste x.)

1. Pentru i = k + 1 : m

1. r =√x2k + x2i

2. ci = xk/r3. si = xi/r4. xk ← r

Lucrand ca mai sus pentru fiecare coloana, obtinem o metoda de triangu-larizare similara celei din Teorema 3.10.

Teorema 3.17 Fie A ∈ Rm×n, cu m > n. Atunci exista rotatiile Gki, k =1 : n, i = k + 1 : m, astfel ıncat matricea

Gnm . . . Gn,n+1︸︷︷︸pt. col n

. . . G2m . . . G23︸︷︷︸pt. col 2

G1m . . . G13G12︸︷︷︸pt. col 1

A = R (3.28)



Demonstratia decurge ın modul cunoscut si nu o mai detaliem. Triangu-larizarea ortogonala este organizata conform schemei

1. Pentru k = 1 : n1. Pentru i = k + 1 : m

1. Calculeaza rotatia Gki (adica scalarii cik si sik), astfelıncat (GkiA)ik = 0

2. Calculeaza A← GkiA

Ca si ın algoritmul TORT, produsul A ← GkiA nu afecteaza coloanele1 : k − 1 (la stanga celei curente). Pentru coloanele k + 1 : n, produsuleste o simpla ınmultire rotatie-vector, descrisa de Propozitia 3.15b. Detaliindoperatiile din schema de mai sus obtinem algoritmul urmator. Din nou, nune preocupam de optimizarea memoriei, lasand neutilizata ın final parteasubdiagonala a matricei A.

Algoritmul 3.5 (TORT ROT – Triangularizare ortogonala cu rotatii)(Se da A ∈ Rm×n, cu m > n. Se calculeaza rotatiile Gki si matriceasuperior triunghiulara R astfel ıncat are loc relatia (3.28). Matricea Rse memoreaza peste triunghiul superior al lui A. Scalarii cik si sik ce de-finesc rotatiile Gki, k = 1 : n, i = k+ 1 : m, se memoreaza ın tablourileC si S, cu dimensiune m× n.)

1. Pentru k = 1 : n1. Pentru i = k + 1 : m

1. r =√a2kk + a2ik

2. cik = akk/r3. sik = aik/r4. akk ← r5. Pentru j = k + 1 : n

1. α← akj2. akj ← cikakj + sikaij3. aij ← −sikα+ cikaij

Comentarii. Bucla 1.1.5 implementeaza ınmultirea rotatiei curente cucoloana j a matricei A. Salvarea temporara a elementului akj ın instructiunea1.1.5.1 este necesara pentru actualizarea din 1.1.5.3, unde este necesara val-oarea initiala a elementului, nu cea modificata ın 1.1.5.2.

Complexitatea algoritmului TORT ROT este cu aproximativ 50% maimare decat cea a algoritmului TORT. Asadar, pentru matrice oarecare, esterecomandata utilizarea reflectorilor. Rotatiile sunt utile atunci cand matriceaA are deja o structura aproape superior triunghiulara, adica sunt putineelemente de anulat pe fiecare coloana.

Se observa ca este folosita doar partea subdiagonala a matricelor C si S.Deci, una dintre ele ar putea fi memorata peste A (de aceea am notat cik, desi

3.7 Rezolvarea problemei CMMP pentru sisteme subdeterminate 79

rotatia este Gki), cu modificari minime ale algoritmului. Necesarul efectiv dememorie este MTORT ROT = 2mn− n(n+ 1)/2.

Pentru rezolvarea problemei CMMP, abordarea este identica celei dinSectiunea 3.4, dar ınlocuind secventa de reflectori cu cea de rotatii, adicaınlocuind (3.19) cu

d = Gnm . . . G13G12b. (3.29)

Algoritmul 3.3 se modifica prin schimbarea procedurii de triangularizare laprimul pas si prin ınlocuirea pasului 2 cu operatiile corespunzand relatiei(3.29) (efectuate pe loc ın b). Pasul 3 este identic, anume rezolvarea sistemuluisuperior triunghiular la care a fost redusa problema CMMP.

Algoritmul 3.6 (CMMP ROT – Calculul solutiei CMMP a unui sis-tem supradeterminat folosind rotatii) (Se dau A ∈ Rm×n, m > n,cu rang maxim, si b ∈ Rn. Se calculeaza solutia CMMP a sistemuluisupradeterminat Ax = b.)

1. [C, S,R] = TORT ROT(A)2. Pentru k = 1 : n

1. Pentru i = k + 1 : m1. α← bk2. bk ← cikbk + sikbi3. bi ← −sikα+ cikbi

3. x = UTRIS(R, b)

Desigur, deoarece triangularizarea determina complexitatea, acest algo-ritm este recomandat atunci cand matricea A are putine elemente subdiago-nale nenule.

3.7 Rezolvarea problemei CMMP pentru sistemesubdeterminate

Ne ocupam aici de problema CMMP subdeterminata, a carei solutie estedescrisa ın Definitia 3.3. Reamintim ca matricea A ∈ Rm×n are mai multecoloane decat linii (m < n) si deci sistemul Ax = b are o infinitate de solutii.(Pastram ipoteza ca rangul matricei A este maxim, rangA = m.)

Teorema 3.18 Daca x0 ∈ Rn este o solutie particulara a sistemului liniarAx = b, atunci multimea tuturor solutiilor sistemului este subspatiul afin (sauvarietatea liniara)

x0 + KerA = {x = x0 + z | z ∈ KerA}.


0

KerA

x0 + KerAx0

x∗

Fig. 3.3 Geometria problemei CMMP subdeterminate.

Demonstratie. Daca Ax0 = b, atunci z = x − x0 satisface Az = 0, deciz ∈ KerA; reciproc, Az = 0 implica A(x0 + z) = b.

Teorema 3.19 Atunci cand matricea A ∈ Rm×n are rang maxim (si m <n), solutia problemei CMMP subdeterminate (3.5) este

x∗ = AT (AAT )−1b. (3.30)

Demonstratie. Luand x0 o solutie particulara a sistemului, Teorema 3.18spune ca solutiile se afla ın subspatiul KerA deplasat cu x0. Geometria prob-lemei e ilustrata ın Figura 3.3, ın care KerA este reprezentat ca o dreapta(adica are dimensiune 1). Solutiile x sunt vectori care pornesc din origine (caorice vectori) si au varful pe subspatiul afin x0 + KerA. Solutia x∗ cu normaminima este ortogonala pe acest subspatiu, deci ortogonala pe KerA, ceea ceimplica x∗ ∈ ImAT . Exista deci y ∈ Rm astfel ıncat x∗ = AT y. Dar Ax∗ = b,deci AAT y = b, adica

y = (AAT )−1b,

ceea ce duce imediat la (3.30).

Ca si ın cazul solutiei normale pentru problema CMMP supradeterminatatratata ın Sectiunea 3.1, formula (3.30) nu este recomandata pentru calcululnumeric, datorita numarului de conditionare al matricei AAT . Solutia sigurase obtine cu transformari ortogonale. Prezentam mai jos varianta bazata petriangularizarea cu reflectori.

Fie U ∈ Rn×n o matrice ortogonala. Notand y = Ux si observand ca‖x‖ = ‖y‖ si

Ax = b ⇔ AUTUx = b ⇔ AUT y = b,

problema se transforma ın aflarea solutiei CMMP pentru ultimul sistem, cunecunoscuta y. Pastram ideea de reducere la forma triunghiulara, dar avandın vedere ca AUT este o matrice ”lata”, alegem U astfel ıncat

AUT = L = [L′ 0],

3.7 Rezolvarea problemei CMMP pentru sisteme subdeterminate 81

unde L ∈ Rm×n este inferior triunghiulara, iar L′ ∈ Rm×m este inferior tri-unghiulara si nesingulara (deoarece rangA = m). Nu este necesar un algoritmspecial pentru aceasta triangularizare ortogonala la dreapta. Notand R = LT ,R′ = L′T , prin transpunere obtinem

UAT = R =

[R′

0

],

adica U = Um . . . U1 este formata din reflectorii care triangularizeaza AT , carese pot calcula cu algoritmul TORT. Comparati cu relatia (3.18). Partitionam

y =

[y′

y′′

],

cu y′ ∈ Rm si y′′ ∈ Rn−m, ceea ce, tinand seama de structura lui L, reducesistemul initial la sistemul inferior triunghiular L′y′ = b, care are solutieunica. Avand ın vedere ca

‖y‖2 = ‖y′‖2 + ‖y′′‖2

si ca y′ este determinat, solutia cu norma minima se obtine pentru y′′ = 0.In fine, pentru a calcula x∗, egalitatea y = Ux conduce la

x = UT y = U1U2 . . . Umy. (3.31)

Am tinut seama de faptul ca reflectorii Uk sunt matrice ortogonale si simet-rice. Asambland toate aceste operatii obtinem algoritmul urmator.

Algoritmul 3.7 (CMMP SUB – Calculul solutiei CMMP a unui sis-tem subdeterminat) (Se dau A ∈ Rm×n, m < n, cu rang maxim,si b ∈ Rn. Se calculeaza solutia CMMP a sistemului subdeterminatAx = b.)

1. [U, β,R] = TORT(AT )2. b = LTRIS(RT , b)3. x(1 : m)← b, x(m+ 1 : n)← 04. Pentru k = m : −1 : 1

1. τ ← (∑mi=k uikxi) /βk

2. Pentru i = k : n1. xi ← xi − τuik

Comentarii. Variabila x a fost folosita si pentru y. Instructiunea 2 esterezolvarea sistemului inferior triunghiular L′y′ = b. Instructiunea 3 contineinitializarea y′′ = 0. Instructiunea 4 implementeaza pe loc relatia (3.31),fiecare iteratie implementand operatia x ← Ukx folosind Propozitia 3.9b.Efortul de calcul este dominat de triangularizarea ortogonala.


3.8 O aplicatie: clasificare binara cu CMMP

Clasificarea pe baza de semnale de antrenare este o componenta importantaa tehnicilor de inteligenta artificiala. Ne ocupam aici de problema clasificariibinare, ın care avem exemple de semnale (vectori) din doua clase si dorimca, pe baza lor, sa cream o regula prin care sa clasificam semnale viitoare.De obicei, semnalele nu sunt cele brute culese de senzori (imagini, sunete),ci sunt obtinute prin extragerea (cu metode uneori complicate!) unor car-acteristici semnificative ale semnalelor brute; pentru imagini, acestea pot fi:histograme de culori, histograme de gradienti (directii semnificative), con-tururi, puncte cheie; pentru sunete: puteri pe diverse intervale de frecventa,frecventa fundamentala si armonice, numar de treceri prin zero, etc.

Notam N1, N2 numarul de vectori de antrenare apartinand claselor C1, re-spectiv C2; vectorii au dimensiune `. O categorie importanta de clasificatoareconstruieste functia liniara

f(v) = cT v + d (3.32)

unde v ∈ R` este un semnal de antrenare, iar c ∈ R` si d ∈ R sunt parametriiınca necunoscuti ai functiei. Ecuatia f(v) = 0 reprezinta suprafata de sep-arare a punctelor dintre cele doua clase; ın cazul (3.32), suprafata este unhiperplan. Scopul (ideal) este de a construi f(v) astfel ıncat{

f(v) > 0, v ∈ C1f(v) ≤ 0, v ∈ C2

(3.33)

Aceasta problema poate fi rezolvata ın mai multe feluri pentru ca pot existamai multe hiperplane de separare. De asemenea, este posibil ca o suprafatade separare sa nu existe, pentru ca exemplele nu sunt suficient de distincte.Cele doua situatii sunt exemplificate ın Figura 3.4. Pentru ` = 2 hiperplanulseparator devine o dreapta, cu verde ın figura. Vectorii de antrenare dincele doua clase sunt reprezentati cu rosu si albastru; cu negru sunt axelede coordonate. In exemplul din stanga, cele doua clase sunt bine separatede dreapta; se vede usor ca solutia nu e unica; aceasta este situatia ideala,mai rar ıntalnita ın practica. In exemplul din dreapta nu exista o dreaptaseparatoare; cea reprezentata, furnizata de metoda descrisa mai jos, este ınsaun compromis acceptabil pentru a generaliza din exemple clasele “rosie” si“albastra”.

Din motivele de mai sus, solutiile pentru gasirea functiei (3.32) se bazeazape diverse euristici care pun problema astfel ıncat sa se obtina o solutie unica,anume cea mai buna aproximare (ıntr-un anume sens) a unei suprafete deseparare. Un clasificator foarte simplu si robust este cel bazat pe CMMP.(Atragem atentia ca exista multe metode de clasificare mai eficiente, dar simai sofisticate.) Aproximam (3.33) cu cerinta

3.8 O aplicatie: clasificare binara cu CMMP 83

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1.5

-1

-0.5

0

0.5

1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Fig. 3.4 Exemple de clasificator CMMP.

f(v) =

{1, v ∈ C1−1, v ∈ C2

(3.34)

Desigur, aceasta relatie nu poate fi satisfacuta exact, de aceea cautam osolutie ın sens CMMP. Notand 1n un vector avand toate elementele egale cu1 si dimensiune n, relatia (3.34) conduce la sistemul CMMP[

V T1 1N1

V T2 1N2

] [cd

]=

[1N1

−1N2

], (3.35)

ın care V1 ∈ R`×N1 si V2 ∈ R`×N2 sunt matrice ale caror coloane sunt vectoriide antrenare din cele doua clase. Rezolvarea sistemului (3.35), a carui matriceare dimensiunea (N1 +N2)× (`+1), produce valorile optime ale parametrilorc si d pentru aceasta metoda.

In mod natural, clasificarea efectiva a unui semnal nou se face tinandseama de semnul functiei (3.32); sau, privind din punctul de vedere CMMP,clasificarea tine cont de care dintre valorile 1 sau −1 este mai aproape val-oarea functiei. Mai precis, daca v este acum un semnal oarecare, decizia declasificare este

v ∈

{C1, sign(cT v + d) = 1,

C2, altfel(3.36)

Asadar, antrenarea clasificatorului presupune rezolvarea unei problemeCMMP, ceea ce poate fi destul de costisitor daca numarul de semnale deantrenare si dimensiunea lor sunt mari; totusi, pentru antrenare ne putempermite alocarea unui timp suficient de mare, deoarece aceasta este o operatiecare se face o singura data. In schimb, utilizarea propriu-zisa a clasificatoruluipresupune doar calculul unui produs scalar, deci este foarte rapida; metodaeste deci potrivita aplicatiilor ın care exista restrictii de timp.

Rezultatele clasificarii semnalelor de antrenare pot fi cuantificate prin ma-tricea de confuzie, al carui element (i, j) reprezinta numarul de semnale dinclasa i care sunt clasate ın clasa j; pe diagonala se afla numarul de sem-


nale clasificate corect, iar pe celelalte pozitii numarul de semnale clasificateincorect. Pentru Figura 3.4 din dreapta (clasa rosie fiind C1 si fiecare clasaavand 20 de semnale), matricea de confuzie este[

20 01 19

].

Toate semnalele rosii sunt clasificate corect, dar unul albastru este clasificatincorect.

Capitolul 4

Valori si vectori proprii

Valorile proprii ale unei matrice patrate, reala sau complexa, contin foartemulta informatie despre matricea respectiva si caracterizeaza proprietatifundamentale ın automatica, precum stabilitatea. In acest capitol descriemmetode pentru calculul numeric al valorilor si vectorilor proprii.

Definitia 4.1 Fie A ∈ Cn×n. Numarul λ ∈ C se numeste valoare proprie amatricei A daca exista un vector v ∈ Cn, v 6= 0, astfel ıncat

Av = λv. (4.1)

In acest caz, v se numeste vector propriu al matricei A, asociat lui λ.

Se observa ca, daca v este vector propriu al matricei A, atunci si u = αveste vector propriu al lui A, corespunzand aceleiasi valori proprii, oricareα ∈ C, α 6= 0. Asadar, vectorii proprii sunt determinati numai ca directie,nu si ca marime. De aceea, de multe ori, se considera ca vectorii proprii suntnormati (adica au norma egala cu 1).

4.1 Fundamente matematice

Teorema 4.2 Fie A ∈ Cn×n. Scalarul λ ∈ C este valoare proprie a matriceiA daca si numai daca matricea λI −A este singulara. Mai mult, matricea Aare exact n valori proprii (numarand separat valorile proprii multiple), carecoincid cu zerourile polinomului caracteristic

p(λ) = det(λIn −A). (4.2)

Daca A ∈ Rn×n, valorile proprii apar ın perechi complex conjugate (iar vec-torii proprii asociati acestor perechi pot fi alesi complex conjugati).

85

86 4 Valori si vectori proprii

Demonstratie. Daca λ este valoare proprie a matricei A, atunci exista unvector v 6= 0 astfel ıncat Av = λv, deci (λIn − A)v = 0 si ın consecintaλIn − A este singulara. Reciproc, daca pentru un λ ∈ C matricea λIn − Aeste singulara, atunci exista v ∈ Cn, v 6= 0, astfel ıncat (λIn−A)v = 0, adicaAv = λv si deci λ este valoare proprie a matricei A.

Deoarece λIn − A este singulara daca si numai daca λ satisface det(λI −A) = 0, iar p(λ) = det(λI−A) este un polinom monic de grad n, din teoremafundamentala a algebrei rezulta ca (4.2) are n zerouri care coincid cu valorileproprii ale matricei A.

Daca A este reala, polinomul caracteristic p(λ) are coeficienti reali si atuncizerourile sale complexe apar ın perechi complex conjugate. Este usor de veri-ficat ca daca u+jw este vector propriu pentru valoarea proprie α+jβ, atunciu− jw este vector propriu pentru α− jβ.

Ca o consecinta directa a faptului ca valorile proprii sunt radacinile ecuatieidet(λIn−A) = 0, observam ca daca A este triunghiulara, ın particular diag-onala, atunci valorile sale proprii sunt chiar elementele diagonale.

Notam prin

λ(A) = {λ1, λ2, . . . , λn} = {λ ∈ C | det(λI −A) = 0} (4.3)

multimea valorilor proprii ale matricei A ∈ Cn×n (cu repetarea valorilorproprii multiple). Multimea σ(A) se numeste spectrul matricei A.

Transformari de asemanare. Suntem interesati de transformarile careconserva spectrul unei matrice date.

Definitia 4.3 Doua matrice A,B ∈ Cn×n sunt numite asemenea daca existao matrice nesingulara T ∈ Cn×n astfel ıncat

B = TAT−1. (4.4)

T se numeste transformare de asemanare. Daca T este unitara (ortogonala,ın cazul real), atunci A si B sunt ortogonal asemenea.

Reamintim ca o matrice Q este unitara daca QHQ = I, unde QH este Her-mitiana matricei Q (transpusa complex conjugata). Trecerea de la cazul realla cel complex, ın general (nu doar pentru ortogonalitate), se face ınlocuindQT cu QH .

Teorema 4.4 Doua matrice asemenea A,B ∈ Cn×n au acelasi spectru,adica λ(A) = λ(B). Mai mult, daca T este matricea de transformare din(4.4) si daca vA este vector propriu al matricei A corespunzator valorii pro-prii λ ∈ λ(A), atunci vectorul

vB = TvA (4.5)

este vector propriu al matricei B, corespunzator aceleiasi valori proprii.

4.2 Forma Schur complexa 87

Demonstratie. Daca A si B satisfac (4.4), atunci

det(λI −B) = det(λI − TAT−1) = det(T (λI −A)T−1) =

= det(T ) det(λI −A) det(T−1) = det(λI −A),

deci λ(A) = λ(B). Daca vA este vector propriu al matricei A corespunzatorvalorii proprii λ, atunci AvA = λvA, sau TAT−1TvA = λTvA, de undeBvB = λvB si deci (4.5) este adevarata.

Matrice simple. O matrice A ∈ Cn×n care are un set complet de n vectoriproprii liniar independenti se numeste simpla. Acesta este cazul generic sitotodata cel ın care proprietatile spectrale ale matricelor sunt mai ”vizibile”.Se poate demonstra ca daca matricea A are n valori proprii distincte, atunciea este simpla.

Teorema 4.5 Fie A ∈ Cn×n o matrice simpla si V ∈ Cn×n matricea alecarei coloane sunt vectorii proprii ai lui A. Atunci

V −1AV = Λ ∈ Cn×n (4.6)

este o matrice diagonala.

Demonstratie. V = [v1 v2 . . . vn] si Avj = λjvj , j ∈ 1 : n, unde λj suntvalorile proprii ale lui A. Atunci

AV = [Av1 Av2 . . . Avn] = [λ1v1 λ2v2 . . . λnvn] =

= [v1 v2 . . . vn] diag(λ1, λ2, . . . , λn) = V Λ.

Deoarece vectorii vj , j ∈ 1 : n, sunt liniar independenti, matricea V estenesingulara si (4.6) este adevarata.

Asadar o matrice simpla poate fi diagonalizata peste C printr-o trans-formare de asemanare. Totusi, diagonalizarea nu functioneaza pentru toatematricele si, ın plus, nu este o optiune buna din punct de vedere numeric.De asemenea, forma Jordan (ale carei detalii nu le mai amintim), desi val-abila pentru toate matricele (diagonalizarea este un caz particular), suferade aceeasi rea conditionare. Sunt preferabile transformarile de asemanareortogonale, care pot aduce matricea initiala la o forma (cvasi)-superior tri-unghiulara, asa cum vom vedea ın sectiunile urmatoare.

4.2 Forma Schur complexa

Discutam ıntai cazul matricelor complexe, care este mai simplu. Acest faptnu este neaparat o surpriza, data fiind relatia dintre valorile proprii ale unei


matrice si radacinile unui polinom. Pe scurt, orice matrice poate fi adusa printransformari de asemanare unitare la forma superior triunghiulara. Aceastaforma este ınsa complexa; elementele ei diagonale sunt valorile proprii alematricei initiale.

Teorema 4.6 (Forma Schur complexa) Oricare ar fi matricea A ∈ Cn×n,exista o matrice unitara Q ∈ Cn×n astfel ıncat matricea

QHAQ = S, (4.7)


Matricea S se numeste forma Schur complexa (FSC) a matricei A, iarcoloanele matricei de transformare Q se numesc vectori Schur ai matricei Aasociati formei Schur S.

Demonstratie. Pasul 1. Fie x1 ∈ Cn un vector propriu al matricei A, cu‖x1‖ = 1; deci Ax1 = λ1x1 pentru un λ1 ∈ σ(A). (Observam din start cademonstratia nu e constructiva, ci se bazeaza pe simplul fapt ca orice matriceare cel putin un vector propriu distinct – ca directie.) Atunci exista o matriceunitara

Q1 = [x1 X1].

Prima coloana este vectorul propriu; celelalte sunt alese (neunic) astfel ıncatQ1 sa fie unitara. Din QH1 Q1 = I rezulta imediat ca XH

1 x1 = 0. Prin calculdirect, constatam ca (am notat cu ∗ blocurile neinteresante)

S1 = QH1 AQ1 =

[xH1 Ax1 ∗XH

1 Ax1 A2

]=

[λ1x

H1 x1 ∗

λXH1 x1 A2

]=

[λ1 ∗0 A2

].

Am reusit deci sa aducem matricea A la forma superior triunghiulara ın primacoloana, printr-o transformare unitara de asemanare.

Pasul 2. Procedam similar cu matricea A2 ∈ C(n−1)×(n−1). Daca x2 esteun vector propriu al ei, adica A2x2 = λ2x2, cu λ2 ∈ σ(A2) ⊂ σ(A), atunciexista o matrice unitara Q2 astfel ıncat

S2 = QH2 A2Q2 =

[λ2 ∗0 A3

].

Notand

Q2 =

[1 0

0 Q2

],

constatam ca

QH2 QH1 AQ1Q2 =

λ1 ∗ ∗0 λ2 ∗0 0 A3

.

4.3 Forma Schur reala 89

Am gasit deci transformarea pentru a obtine forma superior triunghiulara ınprimele doua coloane.

Continuand similar pana la pasul n− 1, obtinem (4.7).

Modul de lucru din demonstratia de mai sus, ın care matricea de interesdevine din ce ın ce mai mica, se mai numeste deflatie.

4.3 Forma Schur reala

Atunci cand matricea A este reala, FSC are dezavantajul de a contine nu-mere complexe (nu doar pe diagonala, ci si ın alte pozitii). Calculul ei, desiınca neprecizat procedural, contine inevitabil operatii cu numere complexe.Este deci de dorit un rezultat similar cu cel din Teorema 4.6, dar care satransforme matricea A ıntr-una suficient de ”simpla” (dar mai ”complicata”decat superior triunghiulara). Avem nevoie de doua rezultate pregatitoare.

Propozitia 4.7 Fie S ∈ Rn×n o matrice bloc superior triunghiulara, custructura

S =

[S11 S12

0 S22

],

unde S11 ∈ Rp×p, S22 ∈ R(n−p)×(n−p). Atunci

σ(S) = σ(S11) ∪ σ(S22), (4.8)

adica valorile proprii ale matricei sunt cele ale blocurilor diagonale.

Demonstratie. Fie λ o valoare proprie a lui S si

x =

[x1x2

]vectorul propriu asociat (deci Sx = λx), cu x1 ∈ Rp, x2 ∈ Rn−p. Tinandseama de structura matricei, obtinem

Sx = λx ⇔{S11x1 + S12x2 = λx1

S22x2 = λx2

Daca x2 6= 0, rezulta ca S22x2 = λx2, deci λ ∈ σ(S22). Daca x2 = 0, atunciS11x1 = λx1, deci λ ∈ σ(S11). Obtinem σ(S) ⊂ σ(S11) ∪ σ(S22). Relatia(4.8) rezulta din egalitatea numarului de valori proprii din σ(S), respectivσ(S11) ∪ σ(S22).

Propozitia 4.8 Daca A ∈ Rn×n, B ∈ Rp×p si X ∈ Rn×p, cu p < n sirangX = p, satisfac relatia


AX = XB, (4.9)

atunci exista Q ∈ Rn×n ortogonala astfel ıncat

QTAQ = S =

[S11 S12

0 S22

], (4.10)

cu S11 ∈ Rp×p.

Demonstratie. Fie

X = QR = Q

[R′

0

],

factorizarea QR a lui X, cu R′ ∈ Rp×p. Deoarece rangX = p, matricea R′

este nesingulara. Notam

QTAQ =

[S11 S12

S21 S22

].

Pentru a demonstra (4.10) trebuie deci sa demostram ca S21 = 0. InlocuindX = QR ın (4.9), obtinem

QTAQR = RB ⇔[S11 S12

S21 S22

] [R′

0

]=

[R′B

0

], (4.11)

de unde rezulta ca S21R′ = 0. Matricea R′ fiind nesingulara, singura solutie

este S21 = 0.

Observatia 4.9 Relatia (4.9) poate fi privita ca o generalizare (formala)a definitiei valorilor proprii. Se constata imediat ca daca x ∈ ImX, atunciAx ∈ ImX. Subspatiul ImX este numit A-invariant.

Observam ın plus ca relatia (4.9) implica σ(B) ⊂ σ(A). Intr-adevar, din(4.11) rezulta si ca S11R

′ = R′B, adica S11 si B sunt asemenea, deci σ(S11) =σ(B). Tinand seama de Propozitia 4.7, rezulta imediat ca σ(B) ⊂ σ(A).

Putem da acum rezultatul omolog Teoremei 4.6 pentru matrice reale.

Teorema 4.10 (Forma Schur reala) Oricare ar fi matricea A ∈ Rn×n, existao matrice ortogonala Q ∈ Rn×n astfel ıncat matricea

QHAQ = S, (4.12)

este bloc superior triunghiulara, cu blocuri diagonale 1 × 1 (corespunzandvalorilor proprii reale) sau 2×2 (corespunzand valorilor proprii complex con-jugate). Matricea S se numeste forma Schur reala (FSR) a matricei A.

Demonstratie. Procedam tot prin deflatie, ca ın Teorema 4.6. Avem ınsadoua situatii. Daca λ este o valoare proprie reala, reducem dimensiunea cu

4.3 Forma Schur reala 91

1 exact ca ın demonstratia Teoremei 4.6. Daca λ = α + jβ ∈ C, notamx = u+ jv vectorul propriu asociat. Cu aceste notatii obtinem

Ax = λx ⇔{Au = αu− βvAv = βu+ αv

⇔ A[u v] = [u v]

[α −ββ α

], (4.13)

adica o egalitate de tip (4.9), cu X = [u v]. Vom demonstra mai jos cavectorii u si v sunt liniar independenti, adica rangX = 2. Putem deci aplicaPropozitia 4.8, obtinand o relatie (4.10) cu

S11 =

[α −ββ α

].

Se poate verifica imediat ca aceasta matrice are ıntr-adevar valorile propriiα± jβ. Am realizat astfel o deflatie cu dimensiune 2, punand ın evidenta unbloc diagonal 2× 2 cu valori proprii complex conjugate. Continuand similarse demonstreaza existenta FSR (4.12).

Mai ramane sa demonstram ca ıntr-adevar rangX = 2. Daca, prin absurd,rangX < 2, atunci u si v sunt vectori colineari: u = γv. Inlocuind ın (4.13),rezulta{γAv = (αγ − β)vAv = (βγ + α)v

⇒ γ(βγ + α) = αγ − β ⇒ β(γ2 + 1) = 0 ⇒ β = 0

Aceasta ar ınsemna ınsa ca valoarea proprie λ este reala, contrar ipotezei.Rezulta ca rangX = 2 si deci demonstratia este ıncheiata.

Observatia 4.11 Diagrama de mai jos ilustreaza doua FSR ale unei matrice7 × 7 cu trei valori proprii reale si patru complexe (doua perechi complexconjugate)

× × × × × × ×× × × × × × ×0 0 × × × × ×0 0 0 × × × ×0 0 0 × × × ×0 0 0 0 0 × ×0 0 0 0 0 0 ×

× × × × × × ×0 × × × × × ×0 0 × × × × ×0 0 0 × × × ×0 0 0 × × × ×0 0 0 0 0 × ×0 0 0 0 0 × ×

FSR (ca si FSC) nu este unica. Ordinea blocurilor diagonale (respectiv, avalorilor proprii) este arbitrara. Deoarece suntem interesati doar de calcululvalorilor proprii, ordinea acestora este irelevanta, deci oricare FSR este la felde buna.


4.4 Metodele puterii si puterii inverse

Metodele iterative prezentate ın aceasta sectiune calculeaza, ın anumiteconditii, un vector propriu al unei matrice date. Desi le discutam ın con-textul deflatiei, unde calculul unui vector propriu ar putea fi ın principiu util(atentie ınsa, deflatia conform demonstratiilor pentru FSR si FSC nu esteeficienta), aceste metode au si importanta de sine statatoare.

4.4.1 Metoda puterii

Metoda puterii calculeaza un vector propriu al unei matrice A ∈ Rn×n, cuconditia ca aceasta sa aiba o valoare proprie dominanta, adica valorile eiproprii sa respecte conditia

|λ1| > |λ2| ≥ . . . ≥ |λn|. (4.14)

Un exemplu suficient de des ıntalnit de matrice care nu respecta aceastaconditie apare cand valoarea proprie maxima ın modul este complexa, avandın vedere ca si complex conjugata ei este valoare proprie. Valoarea propriedominanta este deci obligatoriu reala.

Metoda puterii este definita de urmatoarea iteratie foarte simpla:

x(k+1) =Ax(k)

‖Ax(k)‖. (4.15)

Asadar, vectorul curent x(k) este ınmultit cu matricea A, iar rezultatul estenormat, astfel ıncat ‖x(k+1)‖ = 1.

Teorema 4.12 Fie A ∈ Rn×n o matrice care satisface (4.14). Notamv1 ∈ Rn vectorul propriu corespunzator valorii proprii dominante. Daca vec-torul x(0) are o componenta pe directia v1, atunci metoda puterii descrisa deiteratia (4.15) converge (ın directie) catre v1.

Demonstratie. Pentru simplitate, consideram doar cazul ın care matriceaA are un set de n vectori liniar independenti, v1, v2, . . . , vn. Atunci putemscrie

x(0) =

n∑i=1

αivi,

cu α1 6= 0, conform ipotezei. Tinand seama ca Avi = λivi, obtinem

x(1) =Ax(0)

‖Ax(0)‖= ρ1

n∑i=1

αiλivi,

4.4 Metodele puterii si puterii inverse 93

unde ρ1 e un factor de normare ce asigura ‖x(1)‖ = 1. Prin ınmultiri repetatecu A se ajunge la

x(k) = ρk

n∑i=1

αiλki vi = ρkλ

k1

(α1v1 +

n∑i=2

αiλkiλk1vi

),

Deoarece (λi/λ1)k → 0, i = 2 : n, atunci cand k →∞, rezulta ca x(k) → v1,ın directie, toate celelalte directii disparand.

Observatia 4.13 Doar conditia (4.14) de dominanta este esentiala ın func-tionarea practica a algoritmului. Chiar daca α1 = 0 (caz ın care directia v1nu are teoretic cum aparea ın componenta lui x(k)), faptul ca ın calculele ınvirgula mobila apar erori este de data aceasta benefic, producand de regulamici componente pe directia v1. O data aparuta, o astfel de componenta esteamplificata de metoda puterii si devine dominanta pe masura ce iteratiileavanseaza.

Fiind vorba despre o metoda iterativa ın care aproximatia urmatoareeste calculata folosind-o doar pe cea curenta, implementarea metodei put-erii urmeaza structura cunoscuta a algoritmului Jacobi 2.17. Micile diferentesunt comentate mai jos.

Algoritmul 4.1 (MP – Metoda puterii pentru calculul vectorului pro-priu asociat valorii proprii dominante) (Se dau A ∈ Rn×n; aproximatiainitiala x(0) ∈ Rn; toleranta ε pentru oprire; numarul maxim de iteratiiNmax. Se calculeaza vectorul propriu x.)

1. y ← x(0)/‖x(0)‖2. Pentru k = 1 : Nmax

1. x← Ay2. x← x/‖x‖3. Daca 1− |xT y| < ε atunci Stop


Comentarii. Relatia (4.15) este implementata de instructiunile 2.1 si 2.2,subliniind ca normarea se face eficient dupa ınmultirea matrice-vector. Testulde convergenta nu se bazeaza pe diferenta x− y deoarece un vector propriueste definit de directie; chiar si dupa normare, exista ın continuare doi vec-tori distincti caracterizand acelasi vector propriu, x si −x. Este deci posibilca, desi metoda puterii a convers, sa avem x ≈ −y, adica semnul vectoruluipropriu sa se schimbe la fiecare iteratie, caz ın care x− y nu tinde la zero. Inschimb, produsul scalar tinde la 1 sau −1, adica vectorii x si y devin (aproxi-mativ) colineari. Testul de oprire consta deci ın verificarea colinearitatii, cazın care |xT y| ≈ 1.


Atingerea numarului maxim de iteratii semnifica de regula inexistenta uneivalori proprii dominante, caz ın care metoda puterii nu converge.

Observatia 4.14 Dispunand de un vector propriu (aproximativ) x, de ex-emplu cel furnizat de metoda puterii, se pune ıntrebarea cum se calculeazavaloarea proprie asociata λ. Egalitatea (aproximativa) Ax = λx este de faptun sistem de n ecuatii cu o singura necunoscuta, λ. Fiind vorba despre unsistem supradeterminat, calculam solutia CMMP folosind (3.7) si obtinem

λ =xTAx

xTx, (4.16)

aproximatie (optima) numita si cat Rayleigh. Cand vectorul propriu este nor-mat, relatia devine λ = xTAx.

Metoda puterii este limitata la calculul vectorului propriu asociat valoriiproprii dominante. Pentru extinderea ei, plecam de la urmatoarele rezultate.

Propozitia 4.15 Fie o matrice A ∈ Cn×n cu spectrul σ(A) = {λ1, λ2, ..., λn }.Atunci matricele Ak, k ∈ N∗, A − µIn, µ ∈ C, si, ın cazul ın care A estenesingulara, A−1 au aceiasi vectori proprii cu matricea A si spectrele

σ(Ak) = {λk1 , λk2 , ..., λkn }, (4.17)

σ(A− µIn) = {λ1 − µ, λ2 − µ, ..., λn − µ }, (4.18)

σ(A−1) =

{1

λ1,

1

λ2, ...,

1

λn

}. (4.19)

Demonstratie. Plecand de la Ax = λx, observam ca A2x = λAx = λ2x.Procedand prin inductie se obtine imediat (4.17). Celelalte doua spectrerezulta observand trivial ca, pornind de la aceeasi definitie Ax = λx se ajungela (A− µI)x = (λ− µ)x si A−1x = λ−1x.

Deoarece matricele A si A − µI au aceiasi vectori proprii, putem calculaun vector propriu al lui A aplicand metoda puterii lui A − µI, unde µ ∈ Reste o deplasare aleasa a priori. O scurta analiza arata ca valoarea propriedominanta a lui A− µI poate fi ori λ1 − µ, ori λn − µ. Astfel, putem calculacel mult doi vectori proprii diferiti cu metoda puterii cu deplasare. Este decinevoie de modificari suplimentare daca dorim posibilitatea calcularii oricaruivector propriu. Acestea sunt descrise ın sectiunea urmatoare.

4.4.2 Metoda puterii inverse

Varianta de baza a metodei puterii inverse consta ın aplicarea metodei puteriimatricei inverse A−1 (ın ipoteza destul de generala ca inversa exista). O

4.4 Metodele puterii si puterii inverse 95

iteratie se rezuma deci la calculul

x(k+1) = A−1x(k), (4.20)

urmat de normarea rezultatului. Desigur, calculul din (4.20) se face rezolvandsistemul liniar Ax(k+1) = x(k). Iteratia (4.20) converge la valoarea propriedominanta a matricei A−1, adica λn (cea mai mica valoare proprie a lui A,ın modul, revedeti Propozitia 4.15).

O prima ımbunatatire consta ın adaugarea unei deplasari, adica ın ınlocui-rea lui A−1 cu (A−µI)−1 ın (4.20). In functie de valoarea deplasarii µ, acumorice valoare proprie (λi − µ)−1 poate fi dominanta, atunci cand deplasareaeste suficient de aproape de valoarea proprie (reala) λi. Asadar metoda puteriiinverse cu deplasare poate calcula, ın principiu, orice vector propriu real.Costul acestui avantaj evident fata de metoda puterii este complexitateasporita a iteratiei: o rezolvare de sistem ın loc de, la metoda puterii, o simplaınmultire matrice-vector.

Pentru accelerarea convergentei se foloseste o deplasare variabila, egala cucatul Rayleigh (4.16). Se obtine urmatorul algoritm.

Algoritmul 4.2 (MPI – Metoda puterii inverse cu deplasare Rayleighpentru calculul unui vector propriu) (Se dau A ∈ Rn×n; aproximatiainitiala x(0) ∈ Rn; toleranta ε pentru oprire; numarul maxim de iteratiiNmax. Se calculeaza vectorul propriu x.)

1. y ← x(0)/‖x(0)‖2. Pentru k = 1 : Nmax

1. Calculeaza deplasarea µ = yTAy2. Rezolva sistemul liniar (A− µI)x = y3. Normeaza: x← x/‖x‖4. Daca 1− |xT y| < ε atunci Stop


Comentarii. Chiar daca initializarea x(0) este departe de un vector propriu,deplasarea µ aleasa ın pasul 2.1 este cea mai buna aproximatie posibila aunei valori proprii, ceea ce face ca urmatorul vector calculat ın pasii 2.2–2.3sa fie o aproximatie mai buna a unui vector propriu. Acest proces prin caredeplasarea prin cat Rayleigh si iteratia puterii inverse se potenteaza reciprocconduce la o convergenta foarte rapida. Este tipica obtinerea unui vectorpropriu ın 5-6 iteratii (pentru matrice cu valori proprii bine separate). Spredeosebire, metoda puterii necesita cel putin zeci de iteratii, dar de obicei multmai multe pentru matrice mari.

Observam ın plus ca algoritmul MPI poate calcula orice vector propriuasociat unei valori proprii reale simple, ın functie de initializare. Totusi, estegreu de gasit o modalitate prin care sa se aleaga initializari care sa conducala obtinerea tuturor vectorilor proprii reali.


1

2

3

4

5

6

Fig. 4.1 Exemplu de graf orientat.

In concluzie, metodele prezentate ın aceasta sectiune sunt utile pentrucalculul unui singur vector propriu real al unei matrice reale, anume:

• metoda puterii, pentru cel asociat valorii proprii dominante;

• metoda puterii inverse, pentru cel asociat valorii proprii celei mai apropiatede zero (cea mai mica ın modul);

• metoda puterii inverse cu deplasare Rayleigh, pentru un vector propriureal oarecare.

4.5 Algoritmul PageRank

Prezentam aici un algoritm celebru, cel de stabilire a importantei unei paginiweb. El a fost propus initial ın scientometrie, pentru clasificarea revistelorstiintifice, dar a devenit celebru prin utilizarea ın motorul de cautare Google.In acest context a fost numit PageRank. Asa cum vom vedea ın continuare,el este o varianta a metodei puterii.

Ansamblul paginilor web poate fi modelat printr-un graf orientat, precumcel din Figura 4.1, ın care nodurile sunt pagini iar un arc reprezinta un linkaflat pe pagina sursa cu adresa paginii destinatie. Se pune problema ordonariipaginilor ın functie de importanta lor, fara a dispune de numarul efectiv devizitatori (care e o informatie locala si privata), ci doar pe baza topologieigrafului (care, ın principiu, poate fi descoperita de oricine are acces la web).Cel mai simplu model utilizat este cel al unui cititor care, aflat pe o pagina,sare aleator, cu probabilitati egale, la alta pagina catre care exista link ınpagina curenta.

Notam n numarul de pagini si x ∈ Rn vectorul necunoscut de ponderiasociate paginilor; o pondere mai mare ınseamna o importanta mai mare.Valorile ponderilor nu sunt interesante ın sine, ci doar ordinea lor, care da

4.5 Algoritmul PageRank 97

ordonarea paginilor dupa importanta. Ponderile sunt initializate cu valorioarecare, de exemplu egale, si modificate iterativ astfel: la fiecare moment detimp, un nod da ponderea proprie nodurilor catre care are arce, ın cantitatiegale. De exemplu, ın Figura 4.1, nodul 1 da 1/2 din ponderea sa lui 2 si1/2 lui 6. Privind invers, nodul 1 primeste 1/3 din ponderea lui 3, pondereaıntreaga a lui 4 si 1/2 din ponderea lui 5. Fie A matricea de adiacenta agrafului, ın care aij = 1 daca exista arc de la nodul i la nodul j si aij = 0daca nu exista. Fie D matricea diagonala ın care dii este gradul de iesire alnodului i (adica numarul de arce care ies din acel nod). Notand x(k) vectorulponderilor la pasul k, modelul de mai sus este sintetizat de relatia

x(k+1) = ATD−1x(k). (4.21)

De exemplu, pentru Figura 4.1, avem

ATD−1 =

0 0 1

3 1 12 0

12 0 0 0 0 00 1

2 0 0 0 00 0 1

3 0 12 1

0 0 0 0 0 012

12

13 0 0 0

.

Relatia (4.21) este de fapt metoda puterii aplicata matricei M = ATD−1.Asadar sirul x(k) converge catre vectorul propriu asociat valorii proprii celeimai mari a lui M .

Matricea M are suma elementelor de pe fiecare coloana egala cu 1 si toateelementele sunt nenegative. O astfel de matrice se numeste stochastica (pecoloane). Se observa imediat ca

n∑i=1

x(k+1)i =

n∑i=1

n∑j=1

mijx(k)j =

n∑j=1

x(k)j

n∑i=1

mij =

n∑j=1

x(k)j ,

deci suma ponderilor se conserva; deci nu este nevoie de normare explicitaın (4.21), cum este cazul ın general ın metoda puterii. Initializarea naturala

este x(0)i = 1/n, ceea ce face ca suma ponderilor sa fie permanent egala cu 1.

In plus, daca ponderile initiale sunt nenegative, atunci ele raman nenegativeın fiecare iteratie.

O matrice stochastica (pe coloane) are valoarea proprie maxima egalacu 1, iar celelalte valori proprii sunt de obicei clar mai mici, ceea ce faceca metoda puterii sa convearga relativ repede. In exemplul de mai sus, adoua valoare proprie are modulul egal cu 0.7889. Pentru grafuri mai marisi ın general foarte rare, deoarece o pagina contine de obicei putine link-uri,distanta poate fi si mai mare, ceea ce face ca zeci de iteratii din metodaputerii sa fie suficiente pentru calcularea aproximativa a vectorului propriual valorii proprii dominante, mai ales ca nu dorim decat precizie suficienta


pentru a ordona corect cele mai importante noduri (“coada” clasamentului eneinteresanta). Mai mult, cand graful se modifica datorita aparitiei unor noinoduri sau link-uri, o initializare utilizand ponderile asociate grafului anterioreste aproape de solutie, ceea ce favorizeaza convergenta rapida.

Algoritmul PageRank se bazeaza pe un model ceva mai complicat, carepresupune ca cititorul se poate si opri din navigat (pentru a reıncepe altadata de la o pagina arbitrara). Exista deci un factor de atenuare α ∈ (0, 1):cititorul continua cu un link la ıntamplare cu probabilitate α sau se oprestecu probabilitate 1−α. (Valoarea tipica este α = 0.85.) Presupunand ca sumaponderilor este 1, aceasta modifica iteratia (4.21) ın

x(k+1) =

(αM +

1− αn

E

)x(k),

unde E este o matrice cu toate elementele egale cu 1. Matricea M = αM +1−αn E este si ea stochastica pe coloane, iar convergenta metodei puterii este

mai rapida decat pentru M .

4.6 Algoritmul QR

Algoritmul QR calculeaza forma Schur reala a unei matrice A ∈ Rn×n. Elare o structura de baza extrem de simpla:

1. A0 = A2. Pentru k = 0, 1, 2, . . .

1. Calculeaza o deplasare µk2. Calculeaza factorizarea QR: Ak − µkI = QkRk (QR1)3. Calculeaza Ak+1 = RkQk + µkI (QR2)

Vom discuta aici cateva proprietati ale algoritmului QR ın aceasta formagenerala, sectiunea urmatoare fiind dedicata aspectelor legate de imple-mentarea efectiva.

4.6.1 Asemanare

O prima observatie elementara este ca algoritmul QR construieste o secventade matrice ortogonal asemenea cu A. Intr-adevar, extragand din (QR1)

Rk = QTk (Ak − µkI)

si ınlocuind ın (QR2) obtinem

4.6 Algoritmul QR 99

Ak+1 = QTk (Ak − µkI)Qk + µkI = QTkAkQk.

Aplicand aceasta relatie ıncepand de la k = 0, rezulta

Ak = QTk−1 . . . QT1Q

T0︸︷︷︸

QTk

AQ0Q1 . . . Qk−1︸︷︷︸Qk

. (4.22)

Asemanarea este o cerinta, dar ea nu explica convergenta la FSR. O explicatiepartiala este data ın urmatoarele doua sectiuni.

4.6.2 Legatura cu metoda puterii

Se poate spune ca algoritmul QR contine metoda puterii, ın sensul precizatmai jos.

Introducem (4.22) ın (QR1) si obtinem

QTk (A− µkI)Qk = QkRk, (4.23)

de unde rezulta

(A− µkI)Qk = QkQkRk = Qk+1Rk.

Extragem prima coloana din primul si al treilea termen al egalitatii de maisus, adica ınmultim cu vectorul unitate e1. Tinand seama ca Rk este superior

triunghiulara si deci Rke1 = r(k)11 e1 (indicele superior reprezinta numarul

iteratiei QR), obtinem

(A− µkI)q(k)1 = r

(k)11 q

(k)1 .

Asadar prima coloana a transformarii ortogonale Qk este de fapt obtinutaprin metoda puterii (revedeti (4.15)), cu deplasare.

Presupunand ca procesul converge (prin alegerea judicioasa a deplasarii),

rezulta ca q(k)1 → x, unde x este un vector propriu al matricei A, iar

Qk → [x X], matricea [x X] fiind evident ortogonala. Rezulta deci, ca ındemonstratia Teoremei 4.6, ca algoritmul QR realizeaza asimptotic deflatia

QTkAQk →[λ ∗0 ∗

],

unde λ este valoarea proprie asociata cu x.

In concluzie, algoritmul QR produce rezultatul dorit, adica FSR, ın primacoloana.


4.6.3 Legatura cu metoda puterii inverse

Intr-un sens similar, algoritmul QR contine metoda puterii inverse.

Transpunem relatia (4.23)

QTk (AT − µkI)Qk = RTkQTk ,

de unde obtinem imediat

(AT − µkI)Qk+1 = QkRTk .

Luam ultima coloana a egalitatii de mai sus, adica ınmultim cu en. Tinand

seama ca RTk e inferior triunghiulara si deci RTk en = r(k)nn en, rezulta ca

(AT − µkI)q(k+1)n = r(k)nn q

(k)n . (4.24)

Deci ultima coloana a matricei Qk este obtinuta conform metodei puteriiinverse (4.20).

Presupunand din nou convergenta, rezulta ca q(k)n → y, unde y este un vec-

tor propriu al matricei AT . (Atentie, matricea AT are aceleasi valori propriica A, dar vectori proprii ın general diferiti.) Deci Qk → [Y y] si

QTkAT Qk →

[Y TATY Y TAT yyTATY yTAT y

]=

[∗ 0∗ γ

].

Am tinut seama ca AT y = γy (deci γ e valoare proprie a lui AT si A) si caY T y = 0 (pentru ca [Y y] este ortogonala). Transpunem relatia de mai sussi obtinem

QTkAQk →[∗ ∗0 γ

].

Deci algoritmul QR produce rezultatul dorit (FSR) pe ultima linie.

Aceasta este doar o justificare partiala a convergentei, dar vom vedeamai departe ca implementarea eficienta a algoritmului tine seama si deconvergenta partiala a matricei la FSR.

4.6.4 Alegerea deplasarii

Legatura cu metoda puterii inverse sugereaza si alegerea deplasarii. In (4.24),deplasarea Rayleigh (revedeti algoritmul MPI) se calculeaza cu

q(k)T

n AT q(k)n = eTn QTkA

T Qken = eTnATk en = a(k)nn .

4.7 Implementarea algoritmului QR 101

Asadar, deplasarea optima este elementul din dreapta jos al matricei curenteproduse de algoritmul QR. Nu este necesar nici un calcul special.

4.7 Implementarea algoritmului QR

Deoarece o iteratie a algoritmului QR, atunci cand matricea A este plina,este foarte costisitoare (necesitaO(n3) operatii), matricea se aduce prin trans-formari de asemanare la o forma cat mai aproape de FSR. Aceasta este formasuperior Hessenberg, la care se poate ajunge cu un algoritm direct. Abia apoise trece la faza iterativa; si acolo exista o serie de artificii de implementaredintre care vom discuta doar o parte, fara a da detalii complete.

4.7.1 Aducerea la forma Hessenberg

Descriem aici algoritmul direct de aducere la forma Hessenberg prin trans-formari ortogonale de asemanare. Ideea este foarte simpla: se utilizeaza reflec-tori pentru introducerea de zerouri ın partea inferioara a matricei. Deoarecetransformarea trebuie sa fie de asemanare, reflectorii ınmultesc matricea sidin dreapta. (Reamintim: reflectorii sunt simetrici.)

Teorema 4.16 Fie o matrice A ∈ Rn×n. Atunci exista reflectorii elementariU2, . . . , Un−1 astfel ıncat

Un−1 . . . U2AU2 . . . Un−1 = H, (4.25)

cu H superior Hessenberg. (Deci A si H sunt ortogonal asemenea.)

Demonstratie. Reflectorii se calculeaza ca la triangularizarea ortogonala, cudiferenta ca zerourile sunt introduse sub prima subdiagonala. Toate calculelese fac pe loc ın A. Ilustram ordinea operatiilor si efectul lor printr-o diagrama,existenta reflectorilor fiind evidenta. Notam × un element oarecare si ∗ unelement care se modifica ın urma aplicarii transformarii precedente. Pentruo matrice 4× 4, operatiile sunt:

× × × ×× × × ×× × × ×× × × ×

U2·−→

× × × ×∗ ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗

·U2−→

× ∗ ∗ ∗× ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗

U3·−→

× × × ×× × × ×0 ∗ ∗ ∗0 0 ∗ ∗

·U3−→

× × ∗ ∗× × ∗ ∗0 × ∗ ∗0 0 ∗ ∗

Asadar, U2 se calculeaza pentru a introduce zerouri ın prima coloana ınpozitiile 3 : n. Apoi, ınmultirea A ← U2A modifica toate celelalte coloane


ale matricei, ın pozitiile 2 : n. Inmultirea A ← AU2 se efectueaza ınmultindfiecare linie cu reflectorul; deoarece reflectorul este simetric, operatia esteidentica cu ınmultirea reflector-vector descrisa de Propozitia 3.9b. Inmultireaafecteaza elementele 2 : n de pe fiecare linie. Deci, zerourile introduse ınprima coloana nu sunt afectate si procesul poate continua similar cu a douacoloana, etc. (Se observa totodata ca reducerea direct la forma triunghiularaeste imposibila.)

Demonstratia Teoremei 4.16 conduce la urmatoarea schema de calcul:

1. Pentru k = 1 : n− 21. Calculeaza reflectorul Uk+1 astfel ıncat (Uk+1A)ik = 0, pen-

tru i = k + 2 : n2. Calculeaza A← Uk+1A (se modifica coloanele k + 1 : n)3. Calculeaza A← AUk+1 (se modifica toate liniile)

Detaliind operatiile cu reflectori (revedeti Sectiunea 3.2), obtinem algorit-mul de mai jos.

Algoritmul 4.3 (HESS – Aducere la forma Hessenberg prin trans-formari ortogonale de asemanare) (Se da A ∈ Rn×n. Se calculeaza re-flectorii Uk, k = 2 : n−1, si matricea superior Hessenberg H astfel ıncatare loc relatia (4.25). Matricea H se memoreaza peste A. Vectorii uk cedefinesc reflectorii Uk+1, k = 1 : n− 2, se memoreaza ıntr-un tablou U ;similar, scalarii βk se memoreaza ıntr-un vector β.)

1. Pentru k = 1 : n− 2

1. σ ← sign(ak+1,k) ·√∑n

i=k+1 a2ik

2. uk+1,k = ak+1,k + σ3. Pentru i = k + 2 : m

1. uik = aik4. βk = σuk+1,k

5. ak+1,k ← −σ6. Pentru j = k + 1 : n

1. τ ←(∑n

i=k+1 uikaij)/βk

2. Pentru i = k + 1 : n1. aij ← aij − τuik

7. Pentru i = 1 : n

1. τ ←(∑n

j=k+1 ujkaij

)/βk

2. Pentru j = k + 1 : n1. aij ← aij − τujk

Comentarii. Pentru eficientizarea memoriei, vectorul uk (si nu uk+1) siscalarul βk, k = 1 : n − 2, sunt cei care definesc Uk+1. O folosire multmai eficienta a memoriei ar rezulta daca vectorii uk ar fi memorati pesteA, cu exceptia primelor lor elemente semnificative. Scrierea acestei forme aalgoritmului este un exercitiu pentru cititor.

4.7 Implementarea algoritmului QR 103

4.7.2 Pas QR explicit

Descriem acum o iteratie QR pentru o matrice de forma superior Hessenberg,pe care o notam ın continuare cu A. Principalul obiectiv este de a gasi o formaeficienta a operatiilor. In subsidiar, dar nu lipsit de importanta, este sa neasiguram ca iteratia QR nu afecteaza structura Hessenberg a matricei, pentruca altfel algoritmul HESS descris mai sus devine inutil.

Factorizarea QR a matricei A−µI, care are evident tot structura Hessen-berg, se face ın mod natural cu rotatii, deoarece pe fiecare coloana trebuieanulat un singur element. Nu este necesara calcularea explicita a matriceiortogonale Q, ci doar triangularizarea ortogonala (revedeti Teorema 3.17)

Gn−1,n . . . G23G12︸︷︷︸QT

A = R. (4.26)

Particularizarea algoritmului TORT ROT la cazul matricelor Hessenbergeste un exercitiu relativ simplu, de aceea ilustram procesul de calcul cu odiagrama:

× × × ×× × × ×0 × × ×0 0 × ×

G12·−→

∗ ∗ ∗ ∗0 ∗ ∗ ∗0 × × ×0 0 × ×

G23·−→

× × × ×0 ∗ ∗ ∗0 0 ∗ ∗0 0 × ×

G34·−→

× × × ×0 × × ×0 0 ∗ ∗0 0 0 ∗

A doua operatie cu complexitate importanta a unei iteratii QR este cal-culul produsului RQ+ µI. Ilustram calculul produsului

RQ = RGT12GT23 . . . G

Tn−1,n (4.27)

cu o diagrama ın stilul celei de mai sus:

× × × ×0 × × ×0 0 × ×0 0 0 ×

·GT12−→

∗ ∗ × ×∗ ∗ × ×0 0 × ×0 0 0 ×

·GT23−→

× ∗ ∗ ×× ∗ ∗ ×0 ∗ ∗ ×0 0 0 ×

·GT34−→

× × ∗ ∗× × ∗ ∗0 × ∗ ∗0 0 ∗ ∗

Inmultirea la dreapta cu rotatia Gk,k+1 afecteaza doar coloanele k si k + 1,combinand elemente de pe aceeasi linie a acestora. Daca ambele elementesunt zero, rezultatul ramane zero. Astfel, ınmultirea la dreapta cu GT12 stricadoar zeroul din pozitia (2, 1), dar le lasa neschimbate pe cele din liniile 3 : n.Pe scurt, ınmultirea cu toate rotatiile modifica doar zerourile subdiagonale,rezultatul fiind o matrice Hessenberg. Asadar, forma Hessenberg este conser-vata de o iteratie QR.

Prezentam mai jos algoritmul pentru o iteratie QR, asa cum rezulta dindiagramele de mai sus.


Algoritmul 4.4 (PAS QR EXPL – Iteratie QR explicita) (Se daA ∈ Rn×n, superior Hessenberg. Se calculeaza rotatiile Gk,k+1, k = 1 :n − 1, care efectueaza triangularizarea ortogonala (4.26). Se executaapoi ınmultirile (4.27). Transformarile se executa pe loc ın A.)

1. Calculeaza deplasarea µ si efectueaza A← A− µI2. Pentru k = 1 : n− 1

1. r ←√a2kk + a2k+1,k

2. ck = akk/r3. sk = ak+1,k/r4. akk ← r5. Pentru j = k + 1 : n

1. α← akj2. akj ← ckakj + skak+1,j

3. ak+1,j ← −skα+ ckak+1,j

3. Pentru k = 1 : n− 11. Pentru i = 1 : k + 1

1. α← aik2. aik ← ckaik + skai,k+1

3. ai,k+1 ← −skα+ ckai,k+1

4. Efectueaza A← A+ µI

Comentarii. Pasul 2 este implementarea triangularizarii ortogonale pen-tru matricea deplasata calculata la pasul 1. Pasul 3 este implementareaınmultirilor (4.27). Se observa usor ca, pentru o matrice superior Hessen-berg, complexitatea unei iteratii QR este O(n2), semnificativ mai mica decatO(n3) pentru o matrice generala.

Algoritmul de mai sus nu este folosit ın practica, ci unul echivalent ın careınsa deplasarea nu este efectuata explicit iar pasul QR este realizat printr-osingura trecere prin matricea A. Mai mult, din motive expuse mai jos, seefectueaza pasi QR dubli.

4.7.3 Testul de oprire si deflatia

Algoritmul QR aplicat unei matrice superior Hessenberg poate fi consideratca ıncheiat atunci cand nu exista doua elemente nenule consecutive pe primasubdiagonala, deci nu exista blocuri diagonale mai mari de 2 × 2. Blocurile2× 2 trebuie examinate: daca valorile proprii sunt reale, atunci blocul 2× 2trebuie triangularizat, ceea ce se poate face printr-o transformare (utilizandun vector propriu, ca ın demonstratia Teoremei 4.6).

Totusi, este suficienta aparitia unui singur zero pe prima subdiagonalapentru obtinerea unei matrice bloc superior triunghiulare

4.7 Implementarea algoritmului QR 105[A11 A12

0 A22

].

Propozitia 4.7 ne spune ca, ın acest caz, algoritmul QR poate continua separatın blocurile A11 si A22. Astfel, problema initiala este sparta ın probleme maimici si deci rezolvata mai rapid.

Desigur, elementele subdiagonale nu devin exact zero, ci sunt neglijatecand devin foarte mici. Testul prin care se decide “egalitatea” cu zero estegasirea unui indice i pentru care

|ai+1,i| < ε(|aii|+ |ai+1,i+1|), (4.28)

unde ε este o constanta mica (de exemplu de 100-1000 ori mai mare decatepsilon masina). Asadar, elementul subdiagonal trebuie sa fie mult mai micdecat elementele diagonale vecine cu el.

4.7.4 Pasul dublu

Toata discutia de pana acum a ignorat cazul valorilor proprii complexe.Rationamentul asupra convergentei bazat pe legatura cu metoda puteriifunctioneaza doar cand valorile proprii sunt reale. De asemenea, alegereadeplasarii este inspirata de metoda puterii inverse, deci este din nou valabiladoar pentru valori proprii reale. O solutie nesatisfacatoare este alegerea unordeplasari complexe, deoarece aceasta ar ınsemna calcule efectuate cu numerecomplexe, ceea ce ın general dorim sa evitam.

Exista totusi o iesire din acest impas, care permite deplasari complexe,dar calcule cu numere reale. Ideea este de a efectua doua iteratii QR cudeplasari egale cu valorile proprii ale blocului 2×2 din dreapta jos al matriceicurente. (Desigur, alegerea deplasarii este inspirata din cazul real discutatın Sectiunea 4.6.4.) Fie µ1, µ2 aceste valori proprii, pe care le presupunemcomplex conjugate; daca sunt reale metoda de mai jos ramane valabila. Douaiteratii QR succesive cu aceste deplasari sunt definite prin

A− µ1I = Q1R1

A1 = R1Q1 + µ1IA1 − µ2I = Q2R2

A2 = R2Q2 + µ2I

Observand ca A1 = QT1 AQ1 si ınlocuind ın a treia relatie de mai sus, seobtine

Bdef= (A− µ2I)(A− µ1I) = Q1Q2R2R1.

Observam ın plus ca matricea


B = A2 − (µ1 + µ2)A+ µ1µ2I

este reala. Asadar transformarea ortogonalaQ1Q2, care apare ın relatia dintrematricele initiala si dupa doi pasi, adica A2 = QT2Q

T1 AQ1Q2, este cea din

factorizarea QR a matricei reale B. Calcularea lui B, care ar implica O(n3)operatii, si a factorizarii QR sunt evitate prin tehnica deplasarii implicite,pe care n-o discutam aici. Complexitatea O(n2) se pastreaza si pentru pasuldublu, la fel ca si la cel simplu.

Toate aspectele de implementare practica, discutate mai sus sau nu, con-duc la o complexitate totala de O(n3) pentru algoritmul QR, deoarececonvergenta este atinsa de regula ın O(n) iteratii. Calculul valorilor propriiintra deci ın aceeasi categorie de complexitate cu rezolvarea sistemelor saucu problema CMMP, desi timpul de executie este, ın general, mai mare.

4.8 Matrice simetrice

Matricele simetrice au proprietati care permit simplificarea algoritmului QR.

Teorema 4.17 Daca matricea A ∈ Rn×n este simetrica, adica A = AT ,atunci ea este diagonalizabila printr-o transformare de asemanare ortogonala,adica exista Q ∈ Rn×n ortogonala si D ∈ Rn×n diagonala astfel ıncat

QTAQ = D.

In consecinta, o matrice simetrica are valori proprii reale.

Demonstratie. Teorema 4.10 spune ca exista o matrice ortogonala Q astfelıncat QTAQ = S este bloc superior triunghiulara, cu blocuri diagonale celmult 2 × 2. Transpunand aceasta relatie si tinand seama de simetria lui A,rezulta QTATQ = QTAQ = ST , deci S = ST . Asadar matricea S este blocdiagonala cu blocuri diagonale cel mult 2 × 2. Ramane deci sa demonstramca nu pot exista blocuri 2× 2. Consideram matricea simetrica

S1 =

[α ββ γ

].

Valorile ei proprii sunt radacinile polinomului de gradul doi

det(λI − S1) = λ2 − (α+ γ)λ+ αγ − β2,

al carui discriminant este

(α+ γ)2 − 4(αγ − β2) = (α− γ)2 + 4β2 ≥ 0.

4.9 Calculul vectorilor proprii 107

Deci o matrice simetrica 2 × 2 nu poate avea decat valori proprii reale. Inconcluzie, FSR este simetrica si nu are blocuri 2 × 2 pe diagonala, deci estediagonala.

O consecinta practica este ca ın algoritmul QR pentru matrice simetricesunt suficienti pasii simpli. Partea directa a algoritmului consta ın aducereala forma tridiagonala (ın loc de Hessenberg, ın cazul general). Partea iterativapoate fi si ea implementata mai eficient, tinand seama ca se lucreaza cu omatrice tridiagonala simetrica. In concluzie, algoritmul QR simetric este multmai rapid decat cel general.

4.9 Calculul vectorilor proprii

Pana acum ne-am concentrat atentia mai ales asupra valorilor proprii. Pen-tru calculul vectorilor proprii ai unei matrice A ∈ Rn×n este necesar calculultransformarii ortogonale Q din FSR (4.12). Aceasta implica operatii supli-mentare, deoarece pentru valorile proprii este necesara doar matricea blocsuperior triunghiulara S.

Un vector propriu x al lui A este definit prin Ax = λx. Putem deci scrie

QTAQQTx = λx ⇔ SQTx = λQTx.

Deci, daca v este un vector propriu al lui S, atunci x = Qv. Problema sereduce deci la calculul vectorilor proprii ale unei FSR.

Consideram aici doar cazul cel mai simplu, ın care matricea S este superiortriunghiulara (deci toate valorile proprii sunt reale si se gasesc pe diagonalalui S) si valorile proprii sunt distincte. Fie λ = sjj o valoare proprie si vvectorul propriu (necunoscut) asociat. Fie T = λI − S; vectorul propriu veste solutie a sistemului liniar Tv = 0, care poate fi partitionatT11 y T130 wT

T33

v′vjv′′

= 0 ⇔

T11v′ + vjy + T13v

′′ = 0wT v′′ = 0T33v

′′ = 0, (4.29)

unde v′ ∈ Rj−1, v′′ ∈ Rn−j . Valorile proprii fiind distincte, T11 si T33 suntnesingulare. Din (4.29) rezulta T33v

′′ = 0, deci v′′ = 0. Ecuatia wT v′′ = 0este ıntotdeauna satisfacuta iar vj poate avea o valoare nenula oarecare, sa

zicem vj = 1. (In mod evident, sistemul nu are solutie unica, matricea T

fiind singulara; solutiile sunt ınsa toate colineare.) In fine, sistemul superiortriunghiular T11v

′ = −y are solutie unica, care se poate calcula cu algoritmulUTRIS.

Atunci cand valorile proprii sunt reale, dar exista multiplicitati, matriceaT11 este singulara si are (cel putin) un zero pe diagonala. Procedeul de mai


sus trebuie repetat; pot rezulta mai multe solutii liniar independente (celputin una, cel mult un numar egal cu multiplicitatea valorii proprii).

Cand exista valori proprii complexe, rezolvarea sistemului nedeterminat(S − λI)v = 0 trebuie sa apeleze ıntai la o triangularizare a matricei S − λI,vectorul rezultat fiind desigur complex.

Capitolul 5

Descompunerea valorilor singulare

Descompunerea valorilor singulare (DVS) este unul dintre cele mai utilizaterezultate din algebra liniara. Valorile si vectorii singulari caracterizeaza pro-prietati fundamentale ale unei matrice A ∈ Rm×n. In acest capitol, nu maiexista nicio restrictie asupra matricei A: ea poate fi patrata sau dreptunghiu-lara (cu mai multe linii sau mai multe coloane), iar rangul ei poate avea oricevaloare. DVS exista pentru orice matrice.

5.1 Existenta DVS

Teorema 5.1 Fie o matrice A ∈ Rm×n. Atunci exista matricele U ∈ Rm×msi V ∈ Rn×n ortogonale astfel ıncat

A = UΣV T . (5.1)

Matricea Σ ∈ Rm×n este diagonala si are forma

Σ =

[Σ1 00 0

], (5.2)

undeΣ1 = diag(σ1, σ2, . . . , σr) ∈ Rr×r, (5.3)

siσ1 ≥ σ2 ≥ . . . ≥ σr > 0. (5.4)

Numerele pozitive σi, i = 1 : r, se numesc valori singulare ale matricei A.Coloanele matricei ortogonale U se numesc vectori singulari la stanga, iarcoloanele lui V se numesc vectori singulari la dreapta ai matricei A.

Observam ın plus ca, separand primele r coloane ın U si V , adica scriind

109

110 5 Descompunerea valorilor singulare

U = [U1 U2], cu U1 ∈ Rm×r,V = [V1 V2], cu V1 ∈ Rn×r, (5.5)

relatia (5.1) este echivalenta cu

A = U1Σ1VT1 . (5.6)

Aceasta relatie este uneori numita forma redusa a DVS. De asemenea, dacanotam ui, i = 1 : m, coloanele lui U si vi, i = 1 : n, coloanele lui V , atuncirelatia (5.1) se poate scrie

A =

r∑i=1

σiuivTi . (5.7)

Atat (5.6) cat si (5.7) sunt deseori numite tot DVS.

Demonstratie. Daca A = 0, atunci U = I, Σ = 0 si V = I satisfac conditiileteoremei. Daca A 6= 0, atunci procedam prin deflatie, similar demonstratiilorde la forma Schur (Teoremele 4.6 si 4.10). Prezentam aici doar primul pas.Din definitia normei induse (2.55) aplicata pentru ‖A‖2, exista un vectorx ∈ Rn, cu ‖x‖ = 1, astfel ıncat

σ1def= ‖A‖2 = ‖Ax‖. (5.8)

Definim y = Ax/σ1 si observam ca ‖y‖ = 1. Construim acum matriceleortogonale U1 = [y Y ], V1 = [x X] (calculul ın sine nu e necesar, e suficientca aceste matrice exista). Folosind relatiile de mai sus, obtinem

A1def= UT1 AV1 =

[yTAx yTAXY TAx Y TAX

]=

[σ1 w

T

0 B1

]. (5.9)

Am tinut seama ca Y T y = 0, datorita ortogonalitatii lui U1, si am notatwT = yTAX, B1 = Y TAX. Transformarile ortogonale nu schimba norma 2,deci ‖A1‖2 = ‖A‖2 = σ1. Notam

z =1√

σ21 + ‖w‖2

[σ1w

]si observam ca ‖z‖ = 1. Din definitia normei induse (2.55) rezulta ca

‖A1z‖2 ≤ ‖A1‖2 = σ1.

Pe de alta parte

‖A1z‖2 =1√

σ21 + ‖w‖2

∥∥∥∥[σ1 wT0 B1

] [σ1w

]∥∥∥∥ =1√

σ21 + ‖w‖2

∥∥∥∥[σ21 + ‖w‖2B1w

]∥∥∥∥ ≥ σ1.

5.2 Semnificatii si utilizari ale DVS 111

Cele doua inegalitati de mai sus sunt satisfacute doar daca w = 0. Relatia(5.9) devine

A1 = UT1 AV1 =

[σ1 00 B1

].

Am gasit deci transformari ortogonale care diagonalizeaza A ın prima liniesi coloana. Procedand similar cu B1, obtinem

A2 = UT2 B1V2 =

[σ2 00 B2

]unde σ2 = ‖B1‖2, si asa mai departe pana la diagonalizarea completa amatricei A. Am demonstrat deci (5.1).

Mai ramane sa demonstram (5.4). E suficient sa aratam ca σ1 ≥ σ2, ceeace e evident avand ın vedere ca σ1 = ‖A1‖2, σ2 = ‖B1‖2, iar B1 este osubmatrice a lui A1.

5.2 Semnificatii si utilizari ale DVS

5.2.1 Rangul

Asa cum se vede din (5.6), unde fiecare matrice din produsul din dreaptaare rangul r, matricea A are rangul r. Asadar, rangul este egal cu numarulde valori singulare. DVS este modalitatea cea mai stabila numeric pentrucalculul rangului.

In practica, de cele mai multe ori, matricele care teoretic ar trebui saaiba deficienta de rang, au ın fapt rangul maxim, adica min(m,n), datoritazgomotului (de masurare, de cuantificare, de reprezentare) care afecteazaelementele matricei. Deci calculul DVS produce r = min(m,n). Totusi, unelevalori singulare sunt mici, la nivelul zgomotului. Este datoria utilizatoruluisa decida un prag θ sub care valorile singulare sunt neglijabile. Impunereapragului conduce la calcularea unui rang efectiv r, astfel ıncat σi ≥ θ pentrui = 1 : r si σi < θ pentru i > r.

Exemplul 5.2 O situatie tipica este cea ilustrata ın figura 5.1. Matricea Aeste obtinuta prin A = BCT , unde B,C ∈ R8×4 sunt luate aleator. Asadar,rangul lui A este 4. Totusi, datorita calculului numeric al lui A printr-un pro-dus matriceal, precum si al calculului DVS, se obtin 8 valori singulare. Esteevident ca 4 dintre ele sunt foarte mici, de ordinul erorii de reprezentare ındubla precizie. Fixand un prag convenabil, de exemplu θ = 10−14 (linia rosieorizontala din figura), valorile singulare mici pot fi usor neglijate (automat,fara vizualizare). De regula pragul se fixeaza relativ la valoarea singulara


0 1 2 3 4 5 6 7 8 910-20

10-15

10-10

10-5

100

105

Fig. 5.1 Valori singulare tipice calculate pentru o matrice cu rangul teoretic egal cu 4.

maxima σ1 si tine seama de dimensiunea matricei. Functia rank din MAT-LAB calculeaza rangul aplicand un astfel de prag, egal cu max(m,n)σ1εm,unde εm ≈ 2.2·10−16 este epsilon masina (vezi sectiunea 1.1). Desigur, praguleste relativ la valoarea singulara maxima.

Mai frecvent, matricea este si afectata de zgomot, adica elementele ei suntcunoscute cu precizie mai mica decat cea oferita de formatul de reprezentare.De exemplu, matricea are forma A = BCT + Z, unde Z este o matricealeatoare cu elemente mici; cand acestea au ordinul 10−8, valorile singulareau aspectul din Figura 5.2. In continuare, primele 4 valori singulare suntmult mai mari decat celelalte, dar valorile singulare ”mici” sunt de ordinulde marime al zgomotului. Un prag rezonabil pentru trunchiere este acumθ = 10−6, iar acesta trebuie decis de utilizator fie prin cunoasterea prealabilaa nivelului zgomotului, fie prin observarea empirica a distantei mari ıntrevalorile singulare ”mari” si cele ”mici”.

In cazul ın care se doreste o matrice A cu rangul egal cu r care sa fie ceamai buna aproximatie a matricei A cu acest rang, ea este intuitiv definitade cele mai mari r valori singulare si vectorii singulari asociati. Matricea Apoate fi interpretata drept cea mai buna aproximatie a ”adevaratei” matriceA, obtinuta dupa eliminarea zgomotului.

Propozitia 5.3 Fie A ∈ Rm×n o matrice. Cea mai buna aproximare A curangul r, adica solutia problemei

minA

‖A− A‖2

a.ı. rangA = r(5.10)


0 1 2 3 4 5 6 7 8 910-10

10-8

10-6

10-4

10-2

100

102

Fig. 5.2 Valori singulare tipice pentru o matrice afectata de zgomot, cu rangul efectiv

egal cu 4.

este

A =

r∑i=1

σiuivTi . (5.11)

Nu prezentam demonstratia, dar observam imediat ca aproximatia optimaeste trunchierea DVS (5.7) la primele r valori singulare.

5.2.2 Legatura cu valorile proprii

Tinand seama ca V este ortogonala, din (5.1) si (5.2) rezulta ca

AAT = UΣV TV ΣTUT = U

[Σ2

1 00 0

]UT .

Aceasta este descompunerea valorilor proprii pentru matricea simetrica sipozitiv semidefinita AAT (vezi Sectiunea 4.8). Rezulta deci ca

σ(AAT ) = {σ21 , . . . , σ

2r , 0, . . . , 0}. (5.12)

Deoarece AAT are dimensiune m×m, numarul de valori proprii egale cu zeroeste m− r.

Similar, deoarece U este ortogonala, obtinem


ATA = V ΣTUTUΣV T = V

[Σ2

1 00 0

]V T

siσ(ATA) = {σ2

1 , . . . , σ2r , 0, . . . , 0}. (5.13)

Acum, deoarece ATA are dimensiune n× n, numarul de valori proprii egalecu zero este n− r.

Din cele de mai sus rezulta ca, pentru a obtine valorile singulare, amputea calcula fie AAT , fie ATA (care este mai mica!), calculand apoi valorileproprii ale matricei alese. Valorile singulare sunt radacinile patrate (pozitive)ale valorilor proprii nenule. Aceasta metoda nu este indicata decat eventualpentru matrice mici, deoarece formarea matricelorAAT sauATA ınrautatestenumarul de conditionare, asa cum am vazut la problema CMMP.

5.2.3 Calculul unor norme matriceale

Dupa cum reiese din demonstratia Teoremei 5.1, valoarea singulara maximaeste de fapt norma 2 (numita si norma spectrala) a matricei A:

σ1 = ‖A‖2. (5.14)

Asadar, calculul normei 2 este mai complicat decat cel al altor norme discu-tate ın Sectiunea 2.10.

Deoarece norma Frobenius este invarianta la transformari ortogonale, din(5.1) rezulta ca

‖A‖F = ‖Σ‖F =

√√√√ r∑i=1

σ2i . (5.15)

Aceasta relatie este utila ın calcule doar daca DVS este deja disponibila.Altfel, se foloseste direct definitia (2.54).

5.2.4 Numarul de conditionare

Pentru o matrice A ∈ Rn×n nesingulara, numarul de conditionare este definitde (2.60). Pentru norma 2, el este κ2(A) = ‖A‖2‖A−1‖2. Din (5.1) rezulta ca

A−1 = V Σ−1UT . (5.16)

Atentie, aceasta nu este DVS a matricei A−1, deoarece valorile diagonale alematricei Σ−1 sunt ordonate crescator. Este clar ınsa ca valorile singulare ale


lui A−1 sunt σ−1i . Rezulta ca ‖A−1‖2 = σ−1n , aceasta fiind cea mai marevaloare singulara a lui A−1. Numarul de conditionare este deci

κ2(A) =σ1σn. (5.17)

In concluzie, o matrice este rau conditionata (κ2(A)� 1) daca σ1 � σn.

Observatia 5.4 Fie J ∈ Rn×n matricea cu 1 pe antidiagonala si zero ın rest.Observam ca, dat un vector v ∈ Rn, vectorul Jv contine elementele lui v ınordine inversa. Matricea J este ortogonala si simetrica. Relatia (5.16) esteechivalenta cu

A−1 = V J(JΣ−1J)JUT . (5.18)

Aceasta este DVS a matricei A−1, deoarece JΣ−1J are elementele diagonaleordonate descrescator, iar V J si JUT sunt ortogonale.

5.2.5 Baze ortogonale pentru subspatii asociate uneimatrice

Din definitia (3.1) a imaginii unei matrice A ∈ Rm×n, rezulta ca aceasta esteformata din vectori de forma Ax, cu x ∈ Rn. Dimensiunea subspatiului ImAeste r, adica rangul matricei (numarul de coloane liniar independente). Din(5.6) putem scrieAx = U1Σ1V

T1 x. DeoareceΣ1 este inversabila iar V1 ∈ Rn×r

(reamintim ca n ≥ r) are coloanele ortogonale, atunci cand x acopera Rn,vectorul Σ1V

T1 x acopera Rr. Rezulta ca ImA = ImU1, deci coloanele lui U1

formeaza o baza ortogonala pentru ImA. Deoarece U2 este completarea lui U1

la o matrice ortogonala, rezulta imediat ca coloanele lui U2 formeaza o bazaortogonala pentru KerAT .

Nucleul matricei A, vezi (3.6), este format din vectori x ∈ Rn pentrucare Ax = 0. Dimensiunea lui este n− r (numarul de coloane ”redundante”,care pot fi exprimate liniar ın functie de cele r liniar independente). RelatiaAx = U1Σ1V

T1 x ne spune acum ca V T1 x = 0; ıntr-adevar, deoarece Σ1 este

inversabila si U1 ∈ Rm×r (cu m ≥ r) are coloanele liniar independente, dacaV T1 x 6= 0, atunci si U1Σ1V

T1 x 6= 0. Rezulta ca x ∈ ImV2 (deoarece V2 este

completarea lui V1 la o matrice ortogonala). Conchidem deci ca coloanele luiV2 formeaza o baza ortogonala pentru KerA. Automat, rezulta ca coloanelelui V1 formeaza o baza ortogonala pentru ImAT .


5.3 Rezolvarea problemei CMMP generale

In capitolul 3 am discutat problema CMMP pentru cazul ın care matricea A ∈Rm×n are rang maxim; au rezultat doua situatii, dupa cum m > n (sistemesupradeterminate) sau m < n (sisteme subdeterminate), ambele rezolvatepornind de la triangularizarea ortogonala a lui A, respectiv AT . Presupunemacum ca rangul poate avea orice valoare. Desigur, dorim ca reziduul ‖b−Ax‖2sa fie minim. Din cauza rangului deficitar sau a subdeterminarii sistemului,reziduul poate fi minimizat de mai multi vectori x. Dintre acestia, solutiaCMMP este vectorul cu norma cea mai mica.

Definitia 5.5 Solutia ın sens CMMP a sistemului Ax = b este acel vectorx∗ pentru care

‖x∗‖2 = min‖b−Ax‖2este minim

‖x‖2. (5.19)

Se observa imediat ca, atunci cand rangul este maxim, aceasta definitiecoincide cu Definitiile 3.1, pentru m > n, respectiv 3.3, pentru m < n.De asemenea, daca m = n si rangul este maxim, adica matricea A este in-versabila, solutia problemei (5.19) este solutia (unica a) sistemului determinatAx = b.

Rezolvarea problemei CMMP generale utilizeaza DVS (5.1). (Deoarecefolosim doar norma euclidiana ın aceasta sectiune, o vom nota simplu cu‖ · ‖.) Folosind ortogonalitatea lui U , transformam

‖b−Ax‖ = ‖b− UΣV Tx‖ = ‖d−Σy‖,

unde d = UT b si y = V Tx. Deoarece V este ortogonala, rezulta ca ‖y‖ = ‖x‖,deci minimizarea lui ‖x‖ poate fi ınlocuita cu minimizarea lui ‖y‖; am efectuato schimbare de variabila care nu modifica solutia CMMP. Tinand seama de(5.2), (5.5) si notand

y =

[y′

y′′

], d = UT b =

[UT1 bUT2 b

],

cu y′ ∈ Rr, observam ca

‖d−Σy‖ =

∥∥∥∥[UT1 b−Σ1y′

UT2 b

]∥∥∥∥ .Deoarece Σ1 este nesingulara, minimul este atins cand prima componenta avectorului de mai sus este zero, adica atunci cand

y′ = Σ−11 UT1 b.

Asasar, minimizarea reziduului determina unic y′. Componenta y′′ rezultaimediat minimizand

5.4 Problema CMMP totale 117

‖y‖2 = ‖y′‖2 + ‖y′′‖2,

ceea ce produce y′′ = 0. Soutia problemei CMMP generale este deci

x∗ = V y = V1y′ = V1Σ

−11 UT1 b.

Am demonstrat deci urmatoarea teorema. (A doua parte a formulei esteimediata prin evidentierea coloanelor matricelor U1 si V1.)

Teorema 5.6 Solutia problemei CMMP generale (5.19) este

x∗ = V1Σ−11 UT1 b =

r∑i=1

uTi b

σivi. (5.20)

Observatia 5.7 Cand matricea A are rangul maxim, formula de mai sus esteechivalenta cu cele deduse ın capitolul 3. De exemplu, pentru cazul suprade-terminat (m > n), atunci cand rangul e maxim (r = n), descompunerea SVD(5.6) are forma A = U1Σ1V , adica V1 = V . Solutia CMMP (3.7) devine

x∗ = (ATA)−1AT b = (V Σ21V

T )−1V Σ1UT1 b.

Simplificari imediate duc la (5.20).

Observatia 5.8 Desi solutia descrisa ın aceasta sectiune este valabila pen-tru toate problemele CMMP, ea trebuie folosita doar atunci cand nu existainformatii a priori despre rang. In multe aplicatii se stie (prin constructie sauprin natura datelor) ca rangul este maxim si atunci se aplica algoritmii dinsectiunea 3, care au complexitate mai mica.

5.4 Problema CMMP totale

Problema CMMP standard descrisa de Definitia 3.1, cu A ∈ Rm×n avandrang maxim si m > n, poate fi reformulata astfel:

minx,f

‖f‖

a.ı. Ax = b+ f(5.21)

Solutia CMMP este cea pentru care termenul liber b este minim modificatcu vectorul f astfel ıncat sa aiba loc egalitatea Ax = b+ f , adica sa obtinemun sistem compatibil. (Ignoram temporar problema unicitatii solutiei CMMPx.) Asa cum am vazut ın Exemplul 3.2, de multe ori doar termenul liber esteafectat de erori, de aceea aceasta formulare este extrem de utila ın practica.

Pot ınsa aparea cazuri ın care si matricea A contine erori, nu doar termenulliber b. Pentru acestea este adecvata problema CMMP totale (CMMPT), for-mulata astfel:


minx,E,f

‖[E f ]‖Fa.ı. (A+ E)x = b+ f

(5.22)

Problema este deci de a gasi acel x pentru care perturbatii minime ale ma-tricei A si termenului liber b fac sistemul consistent.

Gasirea solutiei porneste de la rescrierea restrictiei problemei (5.22) ınforma

([A b] + [E f ])

[x−1

]= 0. (5.23)

Putem interpreta aceasta egalitate ca fiind obtinuta de cea mai mica perturbatie[E f ] care adunata la [A b] produce o matrice cu rang mai mic decat n + 1(adica numarul de coloane al lui [A b]; presupunem cazul nebanal ın careb 6∈ ImA). Altfel spus, [A+E b+ f ] este matricea cu rang mai mic (adica n)cea mai apropiata de [A b]. Propozitia 5.3 ne furnizeaza aceasta matrice. Fie

[A b] =

n+1∑i=1

σiuivTi

DVS a matricei [A b] ın forma (5.7). Cea mai buna aproximare de rang n estetrunchierea DVS la primii n termeni:

[A b] + [E f ] =

n∑i=1

σiuivTi . (5.24)

Mai mult, avand ın vedere (5.24), unde vectorii ui sunt liniar independenti,egalitatea (5.23) conduce la

vTi

[x−1

]= 0, i = 1 : n.

Fiind perpendicular pe primele n coloane ale lui V , vectorul din dreapta estecoliniar cu ultima coloana a lui V , adica[

x−1

]= αvn+1 = α

[vn+1

vn+1,n+1

],

unde vn+1 ∈ Rn contine primele n elemente ale lui vn+1. Rezulta imediat caα = −1/vn+1,n+1 si deci solutia CMMPT este

x = −vn+1/vn+1,n+1.

Asadar, solutia CMMPT depinde doar de vectorul singular dreapta core-spunzator celei mai mici valori singulare.

5.5 Calculul DVS 119

5.5 Calculul DVS

Desi conexiunea dintre valorile singulare ale lui A si cele proprii ale lui AAT

sau ATA nu poate fi exploatata direct, ea sugereaza o legatura cu algoritmulQR pentru matrice simetrice. Fara a intra ın detalii, spunem ca algoritmulDVS lucreaza direct pe matricea A, desi este inspirat de algoritmul QR si-metric si are o structura similara cu a acestuia.

Prezentam aici doar etapa directa a algoritmului, ın care matricea A esteadusa prin transformari ortogonale la forma superior bidiagonala, aceastafiind forma cea mai apropiata de diagonala care poate fi obtinuta printr-unalgoritm direct. Desigur, algoritmul DVS are natura iterativa, problema fiindechivalenta cu una de valori proprii. Pentru simplitate, consideram doar cazulm > n, ideea fiind aceeasi indiferent de dimensiuni.

Teorema 5.9 Fie o matrice A ∈ Rm×n, cu m > n. Atunci exista reflectoriielementari U1, . . . , Un si V2, . . . , Vn−1 astfel ıncat

Un . . . U1︸︷︷︸UT

AV2 . . . Vn−1︸︷︷︸V

= B, (5.25)

cu B superior bidiagonala.

Demonstratie. Utilizarea reflectorilor este naturala ın acest context, sin-gura subtilitate fiind ordinea aplicarii lor. Ilustram transformarea matriceiA printr-o diagrama, detaliile de calcul fiind prezentate mai tarziu. Reflec-torul U1 introduce zerouri subdiagonale pe prima coloana. Apoi, reflectorulV2 este calculat pentru a introduce zerouri pe prima linie, ın coloanele 3 : n.Tranformarile au urmatorul efect (calculele se fac pe loc ın A; reamintim ca× semnifica un element oarecare, iar ∗ un element care se modifica ın urmaaplicarii transformarii precedente):

× × × ×× × × ×× × × ×× × × ×× × × ×

U1·−→

∗ ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗

·V2−→

× ∗ 0 00 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗

Calculul si aplicarea U1 (adica operatia A← U1A) sunt la fel ca ın primul pasal triangularizarii ortogonale. Pentru a calcula V2, tinem seama ca aplicareala dreapta a unei transformari revine la ınmultirea fiecarei linii a matricei cuacea transformare. Notand aT1 prima linie din A, trebuie sa calculam reflec-torul V2 astfel ıncat aT1 V2 are elementele 3 : n nule. Deoarece reflectorul estesimetric, prin transpunere problema revine la a calcula V2 astfel ıncat V2a1are elementele 3 : n nule; este vorba deci de problema standard de calcul alunui reflector, rezolvata de Algoritmul 3.1. Singura diferenta, pur formala,


este ca acum lucram cu un vector linie. Calculul A← AV2 se face pe linii siafecteaza toate elementele matricei, mai putin cele din prima coloana. Devineacum evident motivul aplicarii alternative a transformarilor din stanga si dindreapta: efectul aplicarii mai ıntai a mai multor transformari din stanga arfi anulat de aplicarea lui V2.

Aplicarea alternativa continua pana la Un−2 si Vn−1; pentru exemplulnostru, aceasta ınseamna:

× ∗ 0 00 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗0 ∗ ∗ ∗

U2·−→

× × 0 00 ∗ ∗ ∗0 0 ∗ ∗0 0 ∗ ∗0 0 ∗ ∗

·V3−→

× × 0 00 × ∗ 00 0 ∗ ∗0 0 ∗ ∗0 0 ∗ ∗

In acest moment, partea superioara a matricei are forma dorita. Pentru aelimina si elementele din partea inferioara mai e nevoie de doi reflectoriaplicati la stanga, Un−1 si Un.

Comentarii. Daca m ≤ n, principiul aplicarii alternative a reflectorilorramane, dar finalul este diferit. Daca m = n, e nevoie de un singur reflectorla stanga. Daca m = n− 1, nu mai este nevoie de nici un alt reflector. Dacam < n− 1, e nevoie de un reflector la dreapta.

Demonstratia Teoremei 5.9 conduce la urmatoarea schema de calcul pentrubidiagonalizarea ortogonala:

1. Pentru k = 1 : n− 21. Calculeaza reflectorul Uk astfel ıncat (UkA)ik = 0, pentrui = k + 1 : m

2. Calculeaza A← UkA3. Calculeaza reflectorul Vk+1 astfel ıncat (AVk+1)kj = 0, pen-

tru j = k + 2 : n4. Calculeaza A← AVk+1

2. Pentru k = n− 1 : n1. Calculeaza reflectorul Uk astfel ıncat (UkA)ik = 0, pentrui = k + 1 : m

2. Calculeaza A← UkA

Detaliind operatiile cu reflectori (revedeti Sectiunea 3.2), obtinem algo-ritmul de mai jos.

Algoritmul 5.1 (BIDIAG – Bidiagonalizare prin transformari or-togonale bilaterale) (Se da A ∈ Rm×n, cum > n. Se calculeaza reflectoriiUk, k = 1 : n, si Vk, k = 2 : n−1, si matricea bidiagonala B astfel ıncatare loc relatia (5.25). Matricea B se memoreaza peste A. Reflectorii Uksunt definiti de vectorii uk si scalarii βk, k = 1 : n, memorati ıntr-untablou U , respectiv un vector β. Reflectorii Vk sunt definiti de vectorii

5.5 Calculul DVS 121

vk−1 si scalarii γk−1, k = 2 : n−1, memorati ıntr-un tablou V , respectivun vector γ.)

1. Pentru k = 1 : n

1. σ ← sign(akk) ·√∑m

i=k a2ik

2. ukk = akk + σ3. Pentru i = k + 1 : m

1. uik = aik4. βk = σukk5. akk ← −σ6. Pentru j = k + 1 : n

1. τ ← (∑mi=k uikaij) /βk

2. Pentru i = k : m1. aij ← aij − τuik

7. Daca k < n− 1 atunci

1. σ ← sign(ak,k+1) ·√∑n

j=k+1 a2kj

2. vk+1,k = ak,k+1 + σ3. Pentru j = k + 2 : n

1. vjk = akj4. γk = σvk+1,k

5. ak,k+1 ← −σ6. Pentru i = k + 1 : m

1. τ ←(∑n

j=k+1 vjkaij

)/γk

2. Pentru j = k + 1 : n1. aij ← aij − τvjk

Comentarii. Pentru eficientizarea memoriei, vectorul vk si scalarul γk,k = 1 : n − 2, sunt cei care definesc Vk+1. Desigur, o folosire mult maieficienta a memoriei ar rezulta daca vectorii uk, vk ar fi memorati pesteA, cu exceptia primelor lor elemente semnificative. Scrierea acestei formea algoritmului este un exercitiu pentru cititor. Instructiunea Daca 1.7 estefolosita pentru scurtarea codului; ın general, scrierea explicita a buclelor, caın schema de calcul ce precede algoritmul, aduce un usor castig de viteza sieste preferabila.

Date post:	30-Dec-2019
Category:	Documents
Upload:	others
View:	9 times
Download:	0 times

Bogdan Dumitrescubine precizate de operat˘ii elementare (cum ar adunarea, ^ nmult˘irea, etc.) prin...

Documents