Poznámky z numerické matematiky

matematiky

Petr Tichý

11. listopadu 2021

1 Chování algoritmů v počítačové aritmetice 7

1.1 Počítačová aritmetika . . . 8

1.1.1 Reprezentace čísel v počítači . . . . 8

1.1.2 Operace s čísly . . . 10

1.1.3 Škálování dat a kancelace . . . 12

1.2 Analýza chyb . . . 14

1.2.1 Výpočet skalárního součinu . . . 14

1.2.2 Přímá a zpětná analýza chyb . . . . 18

1.2.3 Stabilita algoritmu . . . 20

1.3 Podmíněnost problému. . . 22

1.3.1 PodmíněnostAx . . . 26

1.3.2 Číslo podmíněnosti matice. . . 27

1.3.3 PodmíněnostAx=b . . . 27

1.4 Přesnost výpočtu . . . 28

2 Schurův rozklad 31 2.1 Úloha nalezení spektraA . . . 32

2.1.1 Podobnostní transformace . . . 33

2.1.2 Unitární transformace . . . 34

2.1.3 Tvar transformované matice . . . 36

2.2 Schurův rozklad. . . 37

2.2.1 Reálný Schurův rozklad . . . 41

2.3 Důsledky Schurovy věty . . . 42

2.3.1 Normální matice . . . 42

2.3.2 Diagonalizovatelné matice . . . 44

3 Ortogonální transformace 46 3.1 Ortogonální transformace . . . 47

3.2 Givensovy rotace . . . 48

3.2.1 Rotace vektoru . . . 48

3.2.2 Nulování složek vektoru . . . 53

3.3 Householderovy reflexe . . . 55

3.3.1 Zrcadleníxnay . . . 57

3.3.2 Nulování prvků vektoru . . . 58

3.4 Matice rotace a reflexe vC . . . 59

3.5 Maticové rozklady . . . 62

3.5.1 QR pomocí reflexí . . . 62

3.5.2 Převod na horní Hessenbergův tvar. 65 3.5.3 Převod na bidiagonální matici. . . . 67

3.5.4 Praktické použití rozkladů. . . 68

4 Gram-Schmidtův proces 70

4.1 Gram-Schmidt a QR rozklad . . . 71

4.2 Algoritmická realizace . . . 77

4.2.1 Projekce a projektory . . . 77

4.2.2 Klasický Gram-Schmidt . . . 78

4.2.3 Modifikovaný Gram-Schmidt . . . . 79

4.2.4 Iterované verze . . . 81

5 Výpočet QR rozkladu 84 5.1 Cena výpočtu . . . 85

5.2 Numerická stabilita. . . 86

5.2.1 Aplikace rotací a reflexí . . . 87

5.2.2 Ortogonalita spočtené maticeQˆ . . 88

5.2.3 Stabilita výpočtu faktoruRˆ . . . 89

5.2.4 Norma reziduaA−QˆRˆ . . . 90

6 LU rozklad 91 6.1 Gaussova eliminace a LU rozklad . . . 94

6.1.1 LU rozklad a řešeníAx=b . . . 97

6.1.2 Proveditelnost a silná regularita . . 99

6.1.3 Choleského rozklad . . . 100

6.2 GE s částečnou pivotací . . . 102

6.2.1 Algoritmus a proveditelnost . . . 103

6.2.2 Maticový zápis GEPP . . . 104

6.3 Numerická stabilita. . . 109

6.4 Praktické otázky . . . 113

6.4.1 Citlivost na změnu pravé strany . . 113

6.4.2 Kvalita aproximace a norma rezidua 114 6.4.3 Škálování . . . 115

6.4.4 Iterační zpřesnění. . . 116

6.4.5 GE a velké řídké matice . . . 120

6.4.6 Řídkost, zaplnění, stabilita . . . 122

6.4.7 ŘešeníAx=b pomocí LU a QR . . 123

7 Singulární rozklad 124 7.1 Singulární rozklad . . . 127

7.1.1 Spektrální a singulární rozklad . . . 130

7.2 Zápisy singulárního rozkladu . . . 133

7.3 Aplikace . . . 134

7.3.1 Normy . . . 134

7.3.2 Číslo podmíněnosti κ(A). . . 135

7.3.3 Pseudoinverze matice . . . 136

7.4 Aproximace matice . . . 139

7.4.1 Numerická hodnost matice . . . 143

7.5 Výpočet SVD . . . 145

7.5.1 Spektrální rozkladAA^∗a A^∗A . . . 145

7.5.2 Golub-Kahan-Reinsch algoritmus . . 146

7.5.3 Výpočetní náročnost SVD . . . 147

8 Problém nejmenších čtverců 148 8.1 Problémy nejmenších čtverců . . . 149

8.1.1 Existence řešení. . . 151

8.1.2 Jednoznačnost řešení . . . 152

8.2 Metody řešení LS . . . 155

8.2.1 Soustava normálních rovnic . . . 156

8.2.2 QR rozklad . . . 157

8.2.3 Rozšířená soustava rovnic . . . 159

8.2.4 Singulární rozklad . . . 162

9 Stacionární iterační metody 164

9.1 Stacionární iterační metody . . . 166

9.1.1 Obecné iterační schéma . . . 166

9.1.2 Podmínky konvergence. . . 168

9.2 Příklady stacionárních metod . . . 173

9.2.1 Richardsonova metoda . . . 173

9.2.2 Jacobiova iterační metoda . . . 174

9.2.3 Gauss-Seidelova metoda . . . 175

9.2.4 Další možnosti štěpení . . . 176

9.3 Věty o konvergenci . . . 179

10 Mocninná metoda 187 10.1 Algoritmus . . . 190

10.2 Konvergence. . . 191

10.3 Modifikace mocninné metody . . . 193

11 Krylovovské metody 195 11.1 Idea krylovovských metod . . . 197

11.2 DimenzeKk(A, v). . . 198

11.3 Arnoldiho algoritmus . . . 201

11.3.1 Souvislost s QR rozkladem . . . 205

11.3.2 Arnoldiho metoda . . . 206

11.4 Lanczosův algoritmus . . . 208

11.4.1 Lanczosova metoda. . . 210

11.4.2 Analýza Lanczosovy metody . . . . 211

11.4.3 Jacobiho matice . . . 213

11.4.4 Lanczos a ortogonální polynomy? . . 214

11.5 Další krylovovské metody . . . 216

Numerická (výpočtová) matematika

dimenzen

problém

diskrétní

∞dimenze

problém

spojitý

diskretizace

relevance

Snaha nalézt vhodnou numerickou aproximaci řešení za použitídostupnýchprostředků.

• Nepřesnostiv datech, stačí přibližné řešení.

• Vývojalgoritmů, numerickáanalýza.

Ze zadaných dat se snažíme spočíst dostatečně dobrou číselnou (nume- rickou) aproximaci řešení. Nemá smysl hledatpřesné řešenídiskrétní úlohy. Stačí dostatečně dobrá aproximace.Zdroje chyb: neurčitá (ne- přesná) data, diskretizace někdy také truncation (useknu Taylorův roz- voj), zaokrouhlovací chyby.

Numerická(numerus, vyčíslení pomocí čísel)analýza(analyzovat, po- chopit, detailně rozebrat). Numerická analýza se zabývástudiem algo- ritmů a metodpro nalezení (přibližného) numerického řešení různých matematických problémů.Numerická lineární algebrastuduje výpo- četní algoritmy lineární algebry, obzvláště maticové výpočty.

1

Chování algoritm ˚ u

v po ˇcíta ˇcové aritmetice

1.1 Počítačová aritmetika

1.1.1 Reprezentace čísel v počítači

double precision (64 bitů)

1 11 52

± exponent mantisa

• znaménko(sign) - 1 bit

• exponent11 bitů - zakóduji čísla0. . .2047nebo

−1023,−1022

| {z }

emin

, . . . ,0, . . . ,1023

| {z }

e_max

,1024

• mantisa- binární číslo například01110011011

• reprezentuječíslo v normalizovaném tvaru

±1,mantisa·2^e je-lie_min≤e≤e_max.

• reprezentacespeciálních hodnot exponent mantisa hodnota emin−1 0 0

emin−1 6= 0 ±0,mantisa·2^e

emax+ 1 0 ±∞

emax+ 1 6= 0 NaN(not-a-number)

• nejmenší a největší vyjádřitelné kladné číslo xmin≈2,2·10⁻³⁰⁸, xmax≈1,8·10³⁰⁸ realminarealmaxv Matlabu

• V[2^e,2^e+1]jsou čísla

±1,mantisa·2^e

rovnoměrně rozložena, vzdálenost 2^e−52 ⇒ vzdále- nost mezi 1 a nejbližším větším reprezentovatelným číslem je2⁻⁵²≈2,2·10⁻¹⁶.

• OznačeníF(floating point numbers).

Formátdouble precision(binary64)nejpoužívanější. Dnešní standardní procesory podporují ještě single precision(binary32)(e = 8, t = 23).

Dříve různé aritmetiky, dnes zafixováno standardem IEEE 754 (Institute of Electrical and Electronic Engineers) dvojková (binární)aritmetika s pohyblivou řádovou čárkou. Další formáty binary16(half precision)→ GeForce FX (nVidia),binary128(quadruple precision)→malá hardwa- rová podpora, softwarová emulace pomocí2×double precision.

1.1.2 Operace s čísly

Výrobci hardwaru splňující standard IEEE 754 garantují, že základní operace mezi čísly zFsplňují určitá pravidla.

• round-to-nearest strategy, označení round : x∈R→x¯∈F,

x_min ≤ |x| ≤x_max, kdex¯ je nejbližší číslo k xzF. xmůže být například výsledek operace mezi dvěma číslyFspočtený s vyšší přesností. Garance

round(x) =x(1 +δ), |δ| ≤u,

kdeuse nazývástrojová přesnostči zaokrouhlovací jednotka. Je-lit počet bitů mantisy, pak

u=1 22^−t.

Pro double precisiont= 52, tj.u≈1,1·10⁻¹⁶.

• Strojové epsilon

ε_mach = 2u ≈ 2,2·10⁻¹⁶

je vzdálenost mezi 1 a nejbližším větším číslem zF.

Pro počítání v FP aritmetice (FPA) platí určitá pravidla, která jsou důležitá z matematického hlediska (garantovaná přijatými standardy a výrobci hardware). Každé reálné číslo (v intervalu dané uvažované repre- zentace) může být reprezentováno s relativní přesností menší než strojová přesnostu(zaokrouhlovací jednotka).

• Standardní model aritmetiky s konečnou přesností (floating point arithmetic - FPA) předpokládá

fl(x ◦ y) = (x ◦ y)(1 +δ), |δ| ≤u, prox∈Fa y∈F, kde◦ →operace +,−,∗, /,√ označenífl(·)→výpočet proveden v FPA

Uvažované vztahy říkají, že daný výpočet je možné provést s relativní přesnostíu. Obecně tedy není možné očekávat, že algoritmus spočte ře- šení s lepší relativní přesností než jeu.

1.1.3 Škálování dat a kancelace

• Sčítání malých a velkých čísel, příklad

∞

X

k=1

1

k ≈ 22,04

∃k0:∀k≥k0 již1/k nezvětší součet.

Velká čísla →velké mezery, provádí senormalizace vektorůw=v/kvk,škálovánímatic diagonálníD.

• Cancellation (vyrušení, ztráta informace). Nechť ˆ

a=a(1 + ∆a), ˆb=b(1 + ∆b), a, bpřesná data,a,ˆ ˆbperturbovaná. Uvažujme

x=a−b a xˆ= ˆa−ˆb.

Jak daleko jexodxˆ v relativním smyslu?

|x−x|ˆ

|x| = | −a∆a+b∆b|

|a−b|

≤ max(|∆a|,|∆b|)|a|+|b|

|a−b| . Je-li|a−b| |a|+|b| ⇒ velká relativní chyba.

Při odečítání dvou blízkých čísel může dojít k vyru- šení obsažené informace a k velké relativní nepřes- nosti spočteného výsledku.

Perturbace v datech může být způsobena například nepřesnými daty úlohy nebo tím, že jsou to výsledky výpočtu v FPA. Kan- celace při odečítání zesílí předchozí nepatrné chyby.

Vyhodnoceníf(x) = ^1−cos_x2 ^x vx= 1,2·10⁻⁸

fl(cosx) = 0.

15

z }| {

9. . .9 89, cosx= 0.

15

z }| { 9. . .9 928 fl

1−cosx x²

≈0,77.

Čitatel spočten relativně nepřesně, x² zesílí chybu (a= ˆa= 1,b= cosx,ˆb= fl(cosx) =b+ ∆b).

cosx= 1−2 sinx

2 2

, f(x) =1 2

sin^x₂

x 2

²

= 0,5.

Na co si dávat pozor pří počítání v konečné aritmetice. Při návrhu algoritmu snaha eliminovat tyto dva jevy - normalizace, vyhni se odečítání blízkých čísel - ztráta přesnosti.

Počítačová aritmetika !

• Reprezentace čísel v počítači

• Operace s čísly

• Škálování dat a kancelace

1.2 Analýza chyb

Jak popsatvliv zaokrouhlovacích chybna výsledek?

1.2.1 Výpočet skalárního součinu

s=x^Ty, x, y∈Fⁿ, 01 s= 0

02 for j= 1 :ndo 03 s=s+xjyj

04 end

Výsledek spočtený algoritmem v FPA→s_n = fl(x^Ty).

Použijeme model aritmetiky

fl(x ◦ y) = (x ◦ y)(1 +δ), |δ| ≤u.

• Nechťsj je mezivýsledek na koncij-tého cyklu, sj= fl(s_j−1+ fl(xjyj)), s0= 0.

• Model aritmetiky⇒

sj = [sj−1+xjyj(1 +εj)](1 +δj),

|εj| ≤u,|δj| ≤u,δ1= 0, indukcí sn =

n

X

j=1

xjyj(1 +εj)

n

Y

i=j

(1 +δi)

| {z }

1+µ_j

=

n

X

j=1

x_jy_j(1 +µ_j).

• Jak omezit velikost|µj|?

LemmaNechť |αi| ≤u,i= 1, . . . , n a nechťnu<0.01.

Potom platí

n

Y

i=1

(1 +αi) = 1 +βn, kde|βn| ≤1.01nu.

Důkaz. Platí

|β_n|=

n

Y

i=1

(1 +α_i)−1

≤(1 +u)ⁿ−1.

Jelikož je1 +x≤e^x prox≥0, platí (1 +u)ⁿ−1 ≤ e^nu−1

= nu+(nu)²

2! +(nu)³ 3! +. . .

< nu

1 +nu 2 +nu

2 2

+. . .

= nu

1−nu/2

≤ nu

0.995 <1.01nu.

Za předpokladů předchozího lemmatu platí

|µj| ≤ 1.01nu.

Celkem s_n =

n

X

j=1

x_jy_j(1 +µ_j)

| {z }

ˆ y_j

=x^Ty+

n

X

j=1

x_jy_jµ_j

| {z }

δˆn

,

kde|µ_j| ≤1.01nu.

• Dva pohledy, jak interpretovat spočtený výsledek sn = xb^Tyb = x^Ty+ ˆδn,

bx=xaybmá prvkybyj =yj(1 +µj).

• První interpretacesn =bx^Tby:Vypočtenésn je přes- ným skalárním součinem perturbovaných vektorů, relativní chyby perturbovaných dat omezeny

kx−xkb

kxk = 0 a ky−ykb

kyk ≤ 1.01nu

| {z }

O(u)

.

• Druhá interpretaces_n=x^Ty+ ˆδ_n:Vypočtenés_n je rovno přesnému skalárnímu součinu původních vek- torů plus chybaˆδn,

|δˆn| ≤ 1.01nu

| {z }

O(u)

|x|^T|y|.

1.01n... nejhorší možný scénář, většinou nadhodnocení.

1.2.2 Přímá a zpětná analýza chyb

Použijeme-li algoritmusa(x) k výpočtu v FPA, můžeme vlivem zao- krouhlovacích chyb získat nepřesné výsledky. Úkolem numerické analýzy je vhoně interpretovat spočtené výsledky a porozumnět chování daného algoritmu.

algoritmus: x7→a(x)

Označmea(x) přesný a fl(a(x))spočtený výsledek. Jaké (největší) chyby se algoritmus může při výpočtu dopustit?

a(ˆx) = fl(a(x)) a(x) fl(a(x))

x ˆx x

• Přímá analýza chyb- popis šíření zaokrouhlova- cích chyb v algoritmu, odhadpřímé chyby

ka(x)−fl(a(x))k. Efektivní odhad možný zřídka.

• Zpětná analýza chyb - snaha interpretovat zao- krouhlovací chyby vzniklé při výpočtu pomocí změn vstupních dat→vyjádření a odhadzpětné chyby.

kx−xkˆ .

Hledáme modifikovaná vstupní data úlohy tak, aby spočtené ře- šení bylopřesnýmřešením té samé úlohy, avšak s modifikovanými vstupními daty.

a(ˆx) a(x)

x

ˆ

zpětnáchyba x přímáchyba výp

očet přesná

aritmetik a

přesná

aritmetik a

1.2.3 Stabilita algoritmu

OznačeníO(u)→neurčité číslo,

|O(u)| ≤Ku,

kdeK nezávisí na datech, může záviset na dimenzi úlohy (lineárně, kvadraticky, kubicky), viz skalární součin.

• Algoritmusa(x)jezpětně stabilní, pokud∀x∃x:ˆ fl(a(x)) =a(bx) a kx−xkb

kxk =O(u).

• Algoritmusa(x)je zpětně stabilní, jestliže se chyby výpočtu způsobené zaokrouhlováním promítnou do malých změn vstupních dat.

Analýza chyb !

• Výpočet skalárního součinu

• Přímá a zpětná analýza chyb

• Stabilita algoritmu

James Hardy Wilkinson

backward error analysis 1919 – 1986

1970 ACM Turing Award For outstanding scientific contributions in the computing

community

• collaborator of Alan Turing

• fellow of the Royal Society

• “In order to give error bounds Wilkinson had tocreate a theoryof error analysis for floating-point computation.”

[G. Forsythe]

• “His unique genius has been to balance mathematical analysisso intimately withpractical computing.”

[G. Forsythe]

• https://nla-group.org/james-hardy-wilkinson/

• Wilkinson Prize - numerical software

1.3 Podmíněnost problému

Soustava lineárních rovnicAx=b.

1 0.99 0.99 0.98

x1

x2

= b1

b2

b=

1.9900 1.9700

→ x=

1 1

˜b=

1.9902 1.9704

→ x˜= 3

−1.02

kb−˜bk

kbk ≈10⁻⁴ vede na kx−xk˜ kxk ≈2.

Wilkinsonův příklad „The Perfidious Polynomial“

p(x)≡

20

Y

k=1

(x−k)

p(x) = a0+a1x+· · ·+a19x¹⁹+x²⁰

˜

p(x) = ˜a₀+a˜₁x+· · ·+˜a₁₉x¹⁹+x²⁰

˜

a_k ≡ a_k(1 + 10⁻¹⁰·randn(−1,1))

Vedou malé změny dat na velké změny řešení?

Fundamentální otázka, kterou je dobré si vždy položit. Hrozí úplné zne- hodnocení řešení kvůli nepřesnostem v datech a zaokrouhlovacím chy- bám. Vlastnost problému.

• Matematický problém

f :X→Y,

X, Y normovanélineární prostory,f spojitézobra- zení,X vstupní data,Y řešení.

• Vlastnosti

– existence,jednoznačnostřešení

– citlivostna malé změny dat (podmíněnost)

Matematický problém lze vidětjako zobrazení, které vstupním datům z normovaného vektorového prostoruXpřiřazuje řešení z normovaného vektorového prostoruY.f je obvykle nelineární, ale bývá alespoň spo- jité. Pro spolehlivé řešení je potřeba znát vlastnosti problému. Z nume- rického hlediska je důležitá citlivost. Geometrická představa: graf funkce s velkým gradient.

• Označení∆x→malá změna dat (perturbace)

• ∆f(x) =f(x+ ∆x)−f(x)změna řešení Číslo podmíněnostiproblémuf vx:

κf(x)≡lim

δ→0 sup

k∆xk≤δ

k∆f(x)k kf(x)k /k∆xk

kxk

Špatně či dobře podmíněnýproblém. Špatně podmíněný problém→malé změnyxvedou na velké změny vf(x).

Ill-conditioned nebo well-conditioned. Význam malý a velký může záviset na úloze a na použité počítačové aritmetice. Například u double precision aritmetiky jsme schopni počítat s relativní přesností10⁻¹⁶, velký tedy může znamenat například řádově10⁶a více.

Dána vektorová norma k · k na Cⁿ a C^m. Generovanou maticovou normounazveme funkcionál g:C^n×m→R,

g(A) ≡ max

x6=0

kAxk

kxk = max

kxk=1kAxk.

• g(A)je norma, značeníkAk.

• NavícmultiplikativitakABk ≤ kAkkBk.

• Z definicekAxk ≤ kAkkxk.

• Je-liA čtvercová, pakkIk= 1.

Tato tvrzení viz. cvičení. Příklady vektorových norem, 2-norma, Frobe- niova (není generovaná).

1.3.1 Podmíněnost Ax

Problémf :x7→Ax, A∈C^n×mdána.

κf(x) = lim

δ→0 sup

k∆xk≤δ

kA(x+ ∆x)−Axk kAxk /k∆xk

kxk

= kxk kAxk lim

δ→0 sup

k∆xk≤δ

kA(∆x)k k∆xk

= kxk kAxkkAk.

Aregulární, pak zkxk=kA⁻¹Axk ≤ kA⁻¹kkAxkplyne κf(x)≤ kAkkA⁻¹k.

1.3.2 Číslo podmíněnosti matice

• NechťAregulární ak · kje maticová norma, pak κ(A)≡ kAkkA⁻¹k

nazvemečíslo podmíněnosti maticeA.

• Progenerované maticové normy je

1 =kIk=kA⁻¹Ak ≤ kAkkA⁻¹k=κ(A).

Nejlépe podmíněné matice→násobek unitární.

• Z předchozího, podmíněnostAxsplňuje κf(x)≤κ(A) ∀x .

1.3.3 Podmíněnost Ax = b

Problémf :b→A⁻¹b,A∈C^n×n regulární dána, pak κ_f(b) =kA⁻¹k kbk

kA⁻¹bk =kA⁻¹kkAxk

kxk ≤κ(A).

Podmíněnost problému !

• PodmíněnostAx

• Číslo podmíněnosti matice

• PodmíněnostAx=b

1.4 Přesnost výpočtu

Algoritmusa(x), (relativní)přesnost výpočtu ka(x)−fl(a(x))k

ka(x)k . Lze o ní něco říci?

algoritmus a(x) stabilita problém

f(x) podmíněnost

• Je-lialgoritmus a řešící problém f v F zpětně sta- bilní, potomfl(a(x)) =a(ˆx)pro nějakéxˆ

kx−xkˆ

kxk =O(u)≡δ.

• Algoritmus řeší daný problém znamenáf(x) =a(x) (v přesné aritmetice). Jelikožfl(a(x)) =a(ˆx), platí

ka(x)−fl(a(x))k

ka(x)k = kf(x)−f(ˆx)k kf(x)k .

• Potom

kf(x)−f(ˆx)k kf(x)k kx−ˆxk kxk

≤ sup

kx−exk kxk ≤δ

kf(x)−f(ex)k kf(x)k kx−xke kxk

≈κ_f(x)

pro „rozumně“ spojitý problémf ka(x)−fl(a(x))k

ka(x)k .κ_f(x)O(u).

• Interpretace: přesnost zpětně stabilního algoritmu může být ovlivněna podmíněností problému.

2

Schur ˚ uv rozklad

Řešíme-li matematickou úlohu, je vhodné hle- dat jejíekvivalentní formulacitak, aby se řešení úlohy stalo zjevným→transformace úlohy.

2.1 Úloha nalezení spektra A

Hledámevlastní číslaA∈C^n×n.

• Přes charakteristický polynom?

• Podobnostní transformacepomocí regulárníS, B=S⁻¹A S ,

A∼B relaceekvivalence, stejný Jordanův tvar.

Praktickéotázky

• Jaké transformaceS používat?

• Na jaký tvarA převést?

Označení

• kxk. . . Euklidovská (dvojková) norma vektorux

• kAk. . . spektrální (dvojková) norma maticeA

• κ(S) =kSkkS⁻¹k. . . číslo podmíněnostiS.

2.1.1 Podobnostní transformace

B =S⁻¹A S Úvaha→jaké maticeS používat?

• Vstupnídata(maticeA) zatížena chybamiE, A= ˆA+E,

Aˆje matice původních dat.

• Hledáme podobnostní transformaci B=S⁻¹A S=S⁻¹A Sˆ +S⁻¹E S

• B∼Aˆaž natransformovanouchybuS⁻¹ES, kS⁻¹E Sk ≤ κ(S)kEk

• Pokudκ(S)1, pakmůžedojít k zesílení chyb ob- sažených v datech. SpektrumB se můžepodstatně lišitod spektraA.ˆ

• Vímeκ(S)≥1, jaké matice splňují κ(S) = 1 ?

2.1.2 Unitární transformace

ČtvercováU jeunitární, jestliže U^∗U =U U^∗=I.

U reálná→ortogonální.

• Pro unitárníU a vektorxplatí kxk=√

x^∗x=√

x^∗U^∗U x=kU xk.

• Pro unitárníU platí kUk= max

kxk=1kU xk = 1 =kU^∗k=kU⁻¹k

⇒ κ(U) = 1.

• Unitární podobnostní transformace (S unitární) kS⁻¹ESk = max

kxk=1kS^∗ESxk

= max

kSxk=1kESxk=kEk.

Poznámky

• NechťU je označuje unitární matici.

• Unitární invariantnost spektrální normy kU Ak=kAk=kAUk.

• Unitární invariantnost Frobeniovy normy:

kAk²_F =

n

X

i,j=1

|a_i,j|²=

n

X

i=1

ka_ik²= trace(A^∗A) implikuje

kU AkF =kAkF =kAUkF.

• Je-liS =γ U, kdeγ6= 0a U je unitární, pak κ(S) = 1.

Platí i opačná implikace, důkaz později.

2.1.3 Tvar transformované matice

B =S⁻¹A S . ZB chceme „vyčíst“ vlastní čísla.

• Jordanův tvar→S může být špatně podmíněna.

• Postačí, když B bude v horním trojúhelníkovém tvaru→vlastní čísla jsou na diagonále.

Dva požadavky:

1. Vynulovatvšechny poddiagonální prvky.

2. Transformace nesmí zvětšovat podstatně velikosti chyb→unitárnípodobnostní transformace.

Je to proveditelné?

Úloha nalezení spektra A !

• Podobnostní transformace

• Unitární transformace

• Tvar transformované matice

2.2 Schurův rozklad

Schurova věta:Pro libovolnou čtvercovou maticiAexis- tuje unitární matice U tak, že matice

R=U^∗A U

je horní trojúhelníková s vlastními čísly maticeA na dia- gonále v libovolném předepsaném pořadí.

Důkaz:Indukcí podle řádu maticeA.

• Pro matice řádun= 1věta platí.

• n→n+ 1:

• A ∈C(n+1)×(n+1), nechť je dáno uspořádání vlast- ních čísel. Označme λ první vlastní číslo v tomto uspořádání,xje příslušný vlastní vektor,

Ax=λx, kxk= 1.

• Doplňmexna čtvercovou unitárníH = [x, X], H^∗A H =

x^∗Ax x^∗AX X^∗Ax X^∗AX

=

λ b^∗ 0 C

, neboťx^∗Ax=λx^∗x=λaX^∗Ax=λX^∗xje nulový.

• Indukční předpoklad→existuje unitární V:V^∗CV je horní trojúhelníková matice s vlastními čísly na diagonále v předepsaném pořadí.

• Položíme-li

U = [x, XV] =H

1 0 0 V

, pak

U^∗AU =

1 0 0 V^∗

λ b^∗ 0 C

1 0 0 V

=

λ b^∗V 0 V^∗CV

=R

aA=U R U^∗ je hledaný rozklad.

• Rozklad

A=U R U^∗

nazvemeSchurovým rozklademmaticeA.

• Základ algoritmů pro výpočet vlastních čísel matic, funkcí matic aj. (QR algoritmus)

Issai Schur

Schur decomposition 1875 – 1941

• Student of F. G. Frobenius, worked in Germany (Berlin) for most of his life.

• group representations, matrix theory, combinatorics, number theory

• “Schur was a superb lecturer. His lectures were meticu- lously prepared... [and] were exceedingly popular.”

[W. Ledermann]

• Charismatic, Schur built a research group of major im- portance, successors.

• Trouble times in Berlin 1933 - 1939.

Lze Schurův rozklad vypočítat konečným algoritmem?

hledání kořenů polynomu problém

vlastních čísel

p(λ) =λⁿ+a_n−1λⁿ⁻¹+· · ·+a₁λ+a₀ je charakteristickým polynomem (companion) matice

C=



















0 −a₀

1 0 ...

1 . .. ... . .. 0 −a_n−2

1 −a_n−1

















 ,

tj. platí

p(λ) = det(λI−C).

• Věta (Abel, Galois) Obecná polynomiální rovnicen- tého stupně není pron≥5 řešitelná v radikálech.

• Důsledek. Vlastní čísla matice řádu n ≥ 5 nelze obecně spočíst v konečném počtu kroků⇒ani Schu- rův rozklad obecně vypočítat konečným algoritmem.

2.2.1 Reálný Schurův rozklad

• Areálná symetrická→ ∃ reálný Schurův rozklad.

• Areálná→Schurův rozklad obecněkomplexní.

• Snaha přechod do komplexní aritmetiky co nejvíce oddálit. Je-liAreálná, snažíme sepřiblížitSchurovu rozkladu a zároveňzachovat reálnéfaktory.

Definice:ČtvercováT jehorní kvazi-trojúhelníková, po- kud je blokově horní trojúhelníková, kde každý blok na hlavní diagonále je buď1×1 nebo2×2,

T =











T1,1 T1,2 . . . T1,m

0 T2,2 T2,m

... . .. . .. ... 0 . . . 0 Tm,m









 .

Věta. NechťA je reálná čtvercová matice. Potom existují ortogonálníU areálná kvazi-trojúhelníkováT takové, že

T =U^TAU.

Navíc vlastní čísla každého2×2diagonálního bloku matice T tvoří komplexně sdružený pár.

2.3 Důsledky Schurovy věty

2.3.1 Normální matice

Definice:Čtvercová maticeAse nazývánormální, platí-li A^∗A=A A^∗.

Věta (Spektrální věta pro normální matice)A∈C^n×n je normální⇔ ∃unitárníU a diagonálníD tak, že

U^∗A U =D, neboli A=U D U^∗.

Normální matice lze pomocí unitárních podobnostních transformací převést na diagonální tvar.

Ortogonální bázeCⁿ složená z vlastních vektorů Anormální⇐⇒

• ∃unitárníU a diagonálníD= diag(λ₁, . . . , λ_n):

A U =U D, rozepsáno po sloupcích

Auj =λjuj, j= 1, . . . , n.

Vlastní čísla, vlastní vektoryA.

• ∃ortogonální bázeCⁿsloženáz vlastních vektorůA.

Dyadický rozvoj normální matice

• Součin lze zapsat pomocívnějšího součinu A=U D U^∗=

n

X

i=1

λiuiu^∗_i, . . .dyadický rozvojnormální maticeA.

• Z definice 2-normy

kλ_iu_iu^∗_ik=|λ_i|, i= 1, . . . , n.

Příklady normálních matic

• Ajeunitární⇔je-li normální a její vlastní čísla leží najednotkové kružnici.

• A je hermitovská ⇔ je-li normální a všechna její vlastní čísla jsoureálná.

2.3.2 Diagonalizovatelné matice

A nazveme diagonalizovatelnou, pokud ∃ regulární S a diagonálníD:

S⁻¹A S=D.

• AS=SD ⇒vlastní vektorytvoří bázi Cⁿ.

• Tato báze může být špatně podmíněna.

Lze obecnou matici libovolně blízko aproximovat pomocí matice diagonalizovatelé?

Věta (o hustotě diagonalizovatelných matic)∀ A∈C^n×n a pro libovolně malé >0∃diagonalizovatelnáA∈C^n×n s navzájem různými vlastními čísly tak, že kA−Ak< .

Důkaz:

• Schurova věta⇒A=U R U^∗.

• Není-li A diagonalizovatelná → alespoň jedno ná- sobné vlastní číslo. Stačí porušit diagonáluRpomocí diagonální maticeD

R=R+D

aby na diagonáleRbyly různé hodnoty akDk< .

• Uvažujme

A≡U RU^∗.

• A různá vlastní čísla⇒diagonalizovatelná a kA−Ak = kU R U^∗−U(R+D)U^∗k

= kU DU^∗k=kDk< .

Důsledky Schurovy věty !

• Normální matice

• Diagonalizovatelné matice

3

Ortogonální transformace

3.1 Ortogonální transformace

• Unitární transformacenezvětšují chyby ve vstup- ních datech.

• Základní dva typy – Givensovy rotace – Householderovy reflexe

• Budeme uvažovat v Rⁿ (geometrický význam), lze zobecnit doCⁿ.

• Využití: (numericky stabilní) transformace matice na matici s předem zvolenou strukturou.

3.2 Givensovy rotace

3.2.1 Rotace vektoru

Úloha v R²

Chceme sestrojitG(ϕ)∈R^2×2, kterárealizuje pootočení vektoruxo úhelϕproti směru hodinových ručiček.

• Zapíšeme-li vektorxv bázi {e1, e2}, x=

ξ₁ ξ₂

=ξ1

1 0

+ξ2

0 1

=ξ1e1+ξ2e2, je možné vektorG(ϕ)xvyjádřit ve tvaru

G(ϕ)x=ξ₁(G(ϕ)e₁) +ξ₂(G(ϕ)e₂).

• Otáčí-liG(ϕ)bázové vektorye₁ a e₂o úhelϕ, otáčí i libovolný vektorxo úhelϕ.

0

G(ϕ) 1

0

=

cosϕ sinϕ

, G(ϕ) 0

1

=

−sinϕ cosϕ

, tj.

G(ϕ) =

cosϕ −sinϕ sinϕ cosϕ

MaticeG(ϕ)se nazývámatice Givensovy rotace.

Použití

Jsou-li dány vektoryx a y, kxk =kyk 6= 0, potom lze y získat pootočenímx,

y=G(ϕ)x.

Stačí určitcos(ϕ)asin(ϕ),

cos(ϕ) = x^Ty/kxkkyk sin(ϕ) = p

1−cos²(ϕ)

a sestrojitG(ϕ).cos(ϕ)∈[−1,1]odpovídá úhluϕ∈[0, π]

(viz fcearccos) a tudíž jesin(ϕ)≥0.

Úloha v Rⁿ

Rotaci v rovině dané dvojicí jednotkových vektorů {ei, ej}, i < j,o úhelϕve směru odei kej:

Gi,j(ϕ)=































































1 . ..

1

cosϕ −sinϕ 1

. .. 1

sinϕ cosϕ

1 . ..

1





























































 .

Elementární Givensova rotace. Vlastnosti

• je ortogonální,det(Gi,j(ϕ)) = 1(cvičení).

• G_i,j(ϕ)xmění jeni-tý aj-tý prvekx= [ξ₁, . . . , ξ_n]^T,































































1 . ..

1

cosϕ −sinϕ 1

. .. 1

sinϕ cosϕ

1 . ..

1































































Gi,j(ϕ)x =



































ξ1

... ξ_i−1 ξicosϕ−ξjsinϕ

...

ξisinϕ+ξjcosϕ ξ_j+1

... ξ_n



































· · · i-tý řádek

· · · j-tý řádek .

3.2.2 Nulování složek vektoru

Nulování složek vektoru vR²

• Nechťx6= 0.

• Požadujeme, aby platilo cosϕ −sinϕ

sinϕ cosϕ ξ₁ ξ₂

= p

ξ²₁+ξ²₂ 0

neboli

ξ1cosϕ−ξ2sinϕ = q

ξ²₁+ξ₂², ξ₁sinϕ+ξ₂cosϕ = 0,

či

−ξ₂ ξ₁ ξ1 ξ2

sinϕ cosϕ

= p

ξ²₁+ξ²₂ 0

• Řešením dostaneme

sinϕ = − ξ2

pξ₁²+ξ²₂, cosϕ = ξ₁

pξ²₁+ξ₂².

Nulování složek vektoru vRⁿ

• Chceme vynulovatn−1složek vektoru x∈Rⁿ, x −→ ±kxke1.

• Opakovaněaplikujme elementární Givensovy rotace

x=















•

∗ ...

∗

•















→















•

∗ ...

• 0















→ . . . →















±kxk 0

... 0















=y .

Nulujeme prvky na pozicích n, n−1, . . . ,2 (volíme např. roviny rotace span{e1, en},. . ., span{e1, e2})

• Označíme-li jednotlivé elementární Givensovy ro- tace jakoG_1,n, G_1,n−1, . . . , G_1,2, potom

y = Γ1x , kde Γ1 ≡ G1,2. . . G_1,n−1G1,n. Γ1 . . . složená Givensova rotace.

3.3 Householderovy reflexe

• Druhá základní unitární transformace (zrcadlení).

• Nechť je vRⁿdána nadrovina dimenzen−1, kterou popíšeme jejímnormálovým vektorem q,kqk= 1,

H(q)≡ {z∈Rⁿ: z⊥q}.

• Chceme naléztzrcadlový obraz libovolnéhox∈Rⁿ podle nadrovinyH(q)(nadrovina zrcadlení).

• Připomeňme, ortogonální projekcexdo prostoruQ s ortonormální bazíq1, . . . , qk:

x_Q = hx, q1iq1+· · ·+hx, qkiqk

= (q^T₁x)q₁+· · ·+ (q_k^Tx)q_k.

x

x−xq

y = x−2xq

xq

q

0

xq = (qq^T)x, y=x−2xq = (I−2qq^T)x .

• Nechťq∈Rⁿ a kqk= 1. Pak matici H(q) =I−2qq^T ∈R^n×n

nazýváme maticí Householderovy reflexe vzhledem k nadroviněH(q)určené normálovým vektoremq.

• Vlastnosti:H(q)je ortonormální a symetrická, platí H²(q) =I,det(H(q)) =−1 (cvičení).

3.3.1 Zrcadlení x na y

Nechť jsou dány dva různé vektory x ∈ Rⁿ a y ∈ Rⁿ, kxk=kyk, a nechť

q1≡ x−y

kx−yk, q2≡ x+y kx+yk. Potom

H(q₁)x=y, H(q₂)x=−y.

Důkaz:Vektorx−y jekolmýk nadrovině zrcadlení vek- toruxna vektory. Podobně, vektorx+yje kolmý k nadro- vině zrcadlení vektoruxna vektor−y.

3.3.2 Nulování prvků vektoru

Dánx∈Rⁿ a y=±kxke1, hledámeH(q).

• Dle předchozího, pro q1= x− kxke1

kx− kxke₁k, q2= x+kxke1

kx+kxke₁k je

H(q₁)x=kxke₁, H(q₂)x=−kxke₁.

• Je lepší zvolitq₁neboq₂?

x± kxke1 =











ξ₁± kxk ξ₂

... ξn











• Je-liξ1≥0,ξ1≈ kxk, může dojít kekancelaci ξ1− kxk ≈0

a následnémudělení malou normoukx−kxke1k ≈0.

• Kvůli numerické stabilitě, proξ₁≥0volíme q≡ x+kxke1

kx+kxke₁k a proξ₁<0 použijemeq= _kx−kxke^x−kxke1

1k.

3.4 Matice rotace a reflexe v C

Ztrácí geometrické vlasnosti, ale výpočetně OK.

Givensovy rotace vC G=

c −s

s c

, |s|²+|c|²= 1.

Nulování druhé složkyx= [ξ1, ξ2]^T pomocí s=− ξ₂

kxk, c= ξ₁ kxk. Platí

c −s

s c

ξ₁ ξ₂

= 1 kxk

ξ₁ξ1+ξ₂ξ2

−ξ2ξ1+ξ1ξ2

= kxk

0

.

Householderovy reflexe vC

Nadrovina určená normálovým vektoremq,kqk= 1, H(q) ≡ {z∈Cⁿ : z^∗q= 0}.

• Uvažujme unitární hermitovskouH(q) H(q) ≡ I−2qq^∗∈C^n×n.

• Umíme zobrazovatxnay stejné délky→nulování.

Zobrazeníxna násobeke1, například nakxke1, pro q= x− kxke1

kx− kxke1k.

• Abychome se vyhnulikancelaci, volme pro ξ1=|ξ1|e^iα, 0≤α <2π , vektor

q= x+e^iαkxke₁ kx+e^iαkxke1k. První prvek vektoruqje pak roven

(|ξ1|+kxk)e^iα kx+e^iαkxke1k.

Alston Scott Householder

Householder transformation 1904 – 1993

American mathematician mathematical biology

numerical analysis

• A pioneer of numerical linear algebra.

• “In the 1950s our knowledge of this topic was in a rather chaotic state. A large number of algorithms had been developed but no systematic study of their inter- relationships had been undertaken. It is primarily due to the work of Householder that order has emerged from this chaos.” [J. Wilkinson]

• The theory of matrices in numerical analysis, 1964

• Householder (Gatlinburg) Symposia, Householder Prize (best Ph.D. thesis)

3.5 Maticové rozklady

3.5.1 QR pomocí reflexí

A∈C^n×m, A= [a₁, . . . , a_m].

• Analogicky pro složené Givensovy rotace.

• Cíl: vynulování všech poddiagonálních prvků.

• H1 reflexe realizující transformaci

a₁=









 a_1,1 a_2,1 ... an,1











−→









 r_1,1

0 ... 0











• AplikaceH₁ naA, dostaneme (r_1,1= 0⇔a₁= 0)

A=













• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗













−→













∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗













≡A⁽¹⁾,

• AplikaceH2 naA⁽¹⁾=H1A

A⁽¹⁾=













∗ ∗ ∗ ∗ 0 • ∗ ∗ 0 • ∗ ∗ 0 • ∗ ∗ 0 • ∗ ∗













−→













∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 0 ∗ ∗ 0 0 ∗ ∗ 0 0 ∗ ∗













≡A⁽²⁾.

a⁽¹⁾₂ =



















 a⁽¹⁾_1,2 a⁽¹⁾_2,2 a⁽¹⁾_3,2 ...

a⁽¹⁾_n,2





















−→



















 a⁽¹⁾_1,2 r_2,2 0

... ... 0





















• Jak vypadáH2? Nezmění strukturu 1 sloupce.

• Prok≤min{m, n}nulujemen−kpoddiagonálních prvků sloupce a^(k−1)_k . Hk je blokově diagonální se dvěma bloky, první blok jeI_k−1.

• rk,k = 0 ⇔ je-li k-tý sloupec A lineární kombinací předchozíchk−1sloupců A.

• Případn > m, platí

A^(m−1)=













∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 0 ∗ ∗ 0 0 0 • 0 0 0 •













−→













∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 0 ∗ ∗ 0 0 0 ∗ 0 0 0 0













=A^(m)=Hm. . . H1A≡R.

• Hi jsou unitární⇒součin je unitární, pro Q≡H₁^∗H₂^∗. . . H_m^∗ je A=QR .

Tvary QR rozkladu

• A∈C^n×mobecně obdélníková,Qunitární:

A

=

Q R

@

@

@

0

A

=

Q R

@

@

@ 0

• EkonomickýQRrozklad→méně paměťového místa A

=

Qm R_m

@

@

@ 0

ekonomický→méně paměťového místa

3.5.2 Převod na horní Hessenbergův tvar

Cíl: PřevéstA∈C^n×n pomocí unitárních podobnostních transformacína matici v horním Hessenbergově tvaru:

C=















• • • · · · •

• • • · · · • 0 • • · · · • ... . .. . .. . .. ... 0 · · · 0 • •















= [c_i,j]

Definice:C∈C^n×n taková, žeci,j = 0proi > j+ 1, se nazýváhorní Hessenbergova matice.

Hledáme unitárníQtak, že Q^∗AQ=C .

Aplikujeme stejné Householderovy reflexe zleva i zprava.

V prvním kroku aplikujmeH1zleva naA, H1=

1 0 0 H˜1

A=













∗ ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗













−→













∗ ∗ ∗ ∗

∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗













=H₁A .

Násobíme-liH₁^∗ zprava,nenulová struktura se nezmění

H1AH₁^∗=













∗ ∗ ∗ ∗

∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗











 .

Pon−1krocích je

H_n−1. . . H1

| {z }

Q^∗

A H₁^∗. . . H_n−1^∗

| {z }

Q

horní Hessenbergova. Získáme rozklad A = Q C Q^∗.

• Startovací krok pro výpočetSchurova rozkladu.

3.5.3 Převod na bidiagonální matici

Cíl: Převést A ∈ C^n×m pomocí unitárních transformací na matici v bidiagonálním tvaru:

B=



















• • 0 · · · 0 0 • • . .. ... 0 0 • . .. 0 ... . .. . .. . .. • 0 · · · 0 0 •



















NaAaplikujeme střídavě vhodné reflexe zleva a zprava.













• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗

• ∗ ∗ ∗













−→













∗ ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗













=H₁^(L)A .













∗ • • • 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗













−→













∗ ∗ 0 0 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ ∗ ∗













=H₁^(L)AH₁^(R)

• Startovací krok pro výpočetsingulárního rozkladu.

3.5.4 Praktické použití rozkladů

• ÚlohyAx≈b,Ax=b, nalezení spektraA.

• PříkladAx=b,kdeA∈C^n×n je regulární, Ax=b ⇐⇒ QRx=b,

⇐⇒ QRx=QQ^∗b,

⇐⇒ Rx=Q^∗b.

Při výpočtu QR-rozkladu aplikujeme aktuální ro- tace či reflexe nab, nemusímeQuchovávat, jenR.

• Hledání Schurova rozkladu A=QRQ^∗ pomocí tzv.QR algoritmu:

– 1. krok: převod na horníHessenbergův tvar.

– 2. krok:A1=A,iterační procesi= 1,2, . . . Ai=QiRi, Ai+1=RiQi. – Platí

Ai+1=Q^∗_iAiQi=Q^∗_i. . . Q^∗₁A Q1. . . Qi. – Důkaz konvergence např. v knize Watkins.

Ortogonální transformace !

• Givensovy rotace

• Householderovy reflexe

• Matice rotace a reflexe vC

• Nulování složek vektorů

• Ortogonální transformace matice

• Příklad praktického použití rozkladů

4

Gram-Schmidt ˚ uv proces

4.1 Gram-Schmidt a QR rozklad

Opakování:

Ortogonální projekce

vektoruxdo prostoruQs ortonormální bazíq₁, . . . , q_k x_Q = hx, q1iq1+· · ·+hx, qkiqk, x−x_Q ⊥ Q.

Gram-Schmidtův proces

• Jak naléztortonormální báziq1, . . . , qm prostoru span{a1, . . . , am}

takovou, že proi= 1, . . . , k

span{a1, . . . , ai}= span{q1, . . . , qi}?

• a₁, . . . , a_mjsou lin. nezávislé (pro jednoduchost).

• k= 1:span{a1}= span{q1}a kq1k= 1 ⇒ q1= a1

ka₁k

| {z }

r1,1

⇒ a1=r1,1q1.

• k >1:span{a1, . . . , a_k−1}= span{q1, . . . , q_k−1}.

Odak odečtiortogonální projekcia normalizuj z = ak−

k−1

X

i=1

hak, qii

| {z }

r_i,k

qi

qk = z/kzk

|{z}r_k,k

.

• Po rozepsání a1 = r1,1q1

a₂ = r_1,2q₁+r_2,2q₂ ...

a_m = r_1,mq₁+r_2,mq₂+· · ·+r_m,mq_m.

• Označíme-liA= [a₁, . . . , a_m],Q= [q₁, . . . , q_m]a

R≡











r_1,1 r_1,2 · · · r_1,m 0 r_2,2 r_2,m ... . .. . .. ... 0 · · · 0 rm,m











∈C^m×m,

platíA = QR.

Gram-Schmidtův proces počítá ekonomickýQR rozklad A

=

Qm Rm

@

@

@ 0

Proveditelnost

a1, . . . , am lineárně nezávislé ⇒z vždy nenulový, a tudíž rk,k6= 0,k= 1, . . . , m.

Existence a jednoznačnost

• ExistenceQR (i ekonomického) plynez konstrunkce.

• Neníobecnějednoznačný:∀diagonálníD,D^∗D=I A=QR= (QD^∗)(DR) = ˜QR.˜

• diag(R)jen kladná čísla⇒jednoznačnost.

Gram-Schmidt process

Erhard Schmidt 1879 – 1959

German mathematician functional analysis and applied mathematics

• Schmidt studied in Göttingen under Hilbert. In 1917 he got a position at the University of Berlin, he star- ted the Institute of Applied Mathematics.

• Zur Theorie der linearen und nichtlinearen Integralglei- chungen. I. Teil: Entwicklung willkürlicher Funktionen nach Systemen vorgeschriebener. Math. Ann. 1907.

• Schmidt used what is now known as the classical Gram-Schmidt process.

• Schmidt acknowledged that the algorithm was es- sentially the same as that previously used by Gram.

Gram-Schmidt process

Jørgen Pedersen Gram 1850 – 1916

Danish mathematician probability and numerical analysis

• Gram worked for Hafnia Insurance Company.

• Ueber die Entwickelung reeller Funtionen in Reihen mit- telst der Methode der kleinsten Quadrate. J. Reine An- gew. Math. 1883.

• Gram’s theorem (invariant theory, algebraic geome- try) and the Gramian matrix are named after him.

• The orthogonalization algorithm had been used much earlier by other mathematicians, e.g., Laplace or Cauchy!

Gram-Schmidt process

Pierre-Simon, Marquis de Laplace 1749 – 1827

French polymath mathematics, statistics,

physics, philosophy

• One of the most influential scientists of his time.

• He wanted to estimate the mass of Jupiter and Sa- turn from astronomical data of 6 planets,Ax≈b.

• The method which Laplace introduced consists in successively projecting the system of equations orthogonally to a column of the matrixA.

• Laplace uses what is now known as the modified Gram-Schmidt process.