Faculty of Nuclear Sciences and Physical Engineering Czech Technical University in Prague

Tomáš Oberhuber Vlastní cˇ ísla matic ˇ Cásteˇ cný problém vlastních cˇ ísel Mocninná metoda Redukˇcní metoda

Výpoˇcet kompletního spektra matice

Tomáš Oberhuber

Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces

Faculty of Nuclear Sciences and Physical Engineering Czech Technical University in Prague

Householderovy transformace Givensovy rotace QR iterace Hessenbergovy QR iterace

Shrnutí a otázky

1 / 70

Tomáš Oberhuber

1 Vlastní cˇ ísla matic

Vlastní cˇ ísla matic ˇ Cásteˇ cný problém vlastních cˇ ísel Mocninná metoda Redukˇcní metoda

ˇ 2 Cásteˇ cný problém vlastních cˇ ísel Mocninná metoda Redukˇcní metoda

Výpoˇcet kompletního spektra matice Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace QR iterace Hessenbergovy QR iterace

Shrnutí a otázky

3 Výpoˇcet kompletního spektra matice

Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace QR iterace Hessenbergovy QR iterace

4 Shrnutí a otázky 2 / 70

Tomáš Oberhuber

Vlastní cˇ ísla matic - aplikace

Vlastní cˇ ísla matic ˇ Cásteˇ cný problém vlastních cˇ ísel

Vlastní cˇ ísla matic (eigenvalues and eigenvectors) Modelování a analýza vibrací

Mocninná metoda Redukˇcní metoda


Shrnutí a otázky

2 / 70

Tomáš Oberhuber

Vlastní cˇ ísla matic - aplikace


Modelování a analýza vibrací – Tacoma bridge



Shrnutí a otázky

3 / 70

Tomáš Oberhuber Vlastní cˇ ísla matic

Vlastní cˇ ísla matic aeroelasticita

ˇ Cásteˇ cný problém vlastních cˇ ísel Mocninná metoda Redukˇcní metoda

Výpoˇcet kompletního spektra matice Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces

Aeroelastic Phenomena and Related Research - Part 2 – https://www.youtube.com/watch?v=jW8I2hX4GSs


Shrnutí a otázky

4 / 70


Vlastní cˇ ísla matic - eigenfaces Zpracování obrazu – rozpoznání obliˇceju˚



Shrnutí a otázky

5 / 70


Vlastní cˇ ísla matic - spektra atomu˚



Shrnutí a otázky

6 / 70


Výpoˇcet kompletního spektra matice Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace

Vlastní cˇ ísla matic Theorem 1 P n

Bud’ pn (x) = k =0 ak x k polynom stupneˇ n. Potom nalezení jeho koˇrenu˚ je ekvivalentní výpoˇctu spektra Frobeniovy matice ∈ R n,n svázané s polynomem pn (x) a definované jako  −(an−1 /an ) −(an−2 /an ) . . . −(a1 /an ) −(a0 /an )  1 0 ... 0 0   0 1 . . . 0 0 =  .. .. . .. . .. ..  . . .

C

C

0

Givensovy rotace QR iterace

0

...

1



     

0

Hessenbergovy QR iterace

Shrnutí a otázky

Remark 2 ˇ pak plyne, že nelze sestrojit pˇrímou metodu Z Abelovy vety pro výpoˇcet kompletního spektra matice pro n ≥ 5. 7 / 70



Vlastní cˇ ísla matic Remark 3 • Všechny metody pro výpoˇcet spekter matic jsou proto

iterativní. • Je tedy nutné znát nejaký ˇ aposteriorní odhad chyby,

abychom dokázali stanovit vhodné zastavující kritérium pˇri napoˇcítávání iterací. • Apriorní odhady chyb aproximace vlastních cˇ ísel v tomto kurzu vynecháme.


Shrnutí a otázky

8 / 70

Tomáš Oberhuber

Aposteriorní odhad chyby



Theorem 4

A C

ˆ a ~xˆ 6= 0 Necht’ ∈ n,n je hermitovská matice a necht’ λ jsou napoˇcítané aproximace vlastního cˇ ísla λ a vlastního vektoru ~x . Pro residuum

Trojúhelníková metoda

~r =

LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace QR iterace Hessenbergovy QR iterace

Shrnutí a otázky

A~xˆ − λˆ~xˆ,

pak platí ˆ min λ − λi ≤

λi ∈σ(A)

~r

2 .

~ˆ

x 2

9 / 70


ˇ Cásteˇ cný problém vlastních cˇ ísel


Výpoˇcet kompletního spektra matice Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace

• výpoˇcet kompletního spektra je velmi nároˇcná úloha • cˇ asto nám staˇcí nalézt tˇreba je jedno, v absolutní

ˇ vlastní cˇ íslo hodnoteˇ nejvetší • mluvíme pak o cˇ ásteˇcném problému vlastních cˇ ísel • pro tento úˇcel se používá hlavneˇ mocninná metoda

QR iterace Hessenbergovy QR iterace

Shrnutí a otázky

10 / 70

Tomáš Oberhuber

Mocninná metoda


Výpoˇcet kompletního spektra matice Trojúhelníková metoda

a cˇ íselnou posloupnost {ρk }k =1

• postupujeme podle následujícího pˇredpisu (~ x0 lze volit

ˇ libovolne)


∞

• tato metoda konstruuje posloupnost vektoru˚ ~ x (k ) k =1 ∞

~y (k +1) =

A~x (k ),

ρk +1 = ~e1T ~y (k +1) , 1 (k +1) ~y ~x (k +1) = . ρk +1

Shrnutí a otázky

11 / 70

Tomáš Oberhuber

Mocninná metoda



Example 5 ˇ nejvetší ˇ vlastní cˇ íslo Pomocí mocninné metody naleznete matice   1 . 2 = 3

A

Givensovy rotace QR iterace Hessenbergovy QR iterace

Shrnutí a otázky

12 / 70

Tomáš Oberhuber

Mocninná metoda



Shrnutí a otázky

Remark 6 Metoda vlastneˇ funguje tak, že opakovaným násobením ˇ vektoru˚ ~x (k ) do smeru ˇ maticí se zesiluje "prum ˚ et" vlastního vektoru pˇríslušejícího v absolutní hodnoteˇ ˇ ˇ byl nejvetšímu vlastnímu cˇ íslu. Pokud by ale tento "prum ˚ et" od zaˇcátku nulový, metoda se "nemá od cˇ eho odrazit". ˇ není jednoduché. Naštestí ˇ nás Zajistit pˇredpoklady vety zachrání zaokrouhlovací chyby, které zpusobí, ˚ že cˇ asem ˇ do vhodného vlastního vždy vznikne nenulový "prum ˚ et" vektoru a metoda ho následneˇ zaˇcne zesilovat. V praxi je tak nejlepší volit za ~x (0) vektor ze samých jedniˇcek.

A

13 / 70

Tomáš Oberhuber

Mocninná metoda


Theorem 7 Vektor ~x (k ) z mocninné metody lze vyjádˇrit takto


~x (k ) =

Trojúhelníková metoda LR algoritmus

1 ρ1 ρ2 . . . ρk

Ak ~x (0).

QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace QR iterace Hessenbergovy QR iterace

Theorem 8 Až na "normování" jde tedy o aplikaci matice poˇcáteˇcní vektor ~x (0) .

Ak na

Shrnutí a otázky

14 / 70


Výpoˇcet kompletního spektra matice Trojúhelníková metoda LR algoritmus

Mocninná metoda Theorem 9

A = X−1JAX

Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces


Shrnutí a otázky

X

A

QR algoritmus

Householderovy transformace

A C

Necht’ matice ∈ n,n má jedno v absolutní hodnoteˇ ˇ vlastní cˇ íslo λ jednonásobné nebo vícenásobné se nejvetší stejnou algebraickou i geometrickou násobností a ~y (1) , . . . , ~y (r ) jsou pˇríslušné vlastní vektory. Necht’ matice pˇrevádí do Jordanova tvaru, tj.

J

tak, že λ se vyskytuje na prvních r rˇádcích A . Pak pro libovolné ~x (0) takové, že alesponˇ jedna z prvních r složek vektoru ~x (0) je nenulová, platí, že v mocninné metodeˇ konverguje posloupnost ρk k vlastnímu cˇ íslu λ a posloupnost ~x (k ) k pˇríslušnému vlastnímu vektoru.

X

15 / 70



Mocninná metoda Theorem 10

A C

ˇ Necht’ matice ∈ n,n má dveˇ v absolutní hodnoteˇ nejvetší vlastní cˇ ísla λ1 , −λ1 s opaˇcným znaménkem. Necht’ ~y (1) , ~y (2) jsou pˇríslušné vlastní vektory. Necht’ matice pˇrevádí do Jordanova tvaru, tj.

X

A

A = X−1JAX


J


Shrnutí a otázky

tak, že ±λ se vyskytuje na prvních dvou rˇádcích A . Pak pro libovolné ~x (0) takové, že alesponˇ jedna z prvních dvou složek vektoru ~x (0) je nenulová, platí, že v mocninné √ metodeˇ konverguje posloupnost ρ2k ρ2k +1 k absolutní ˇ hodnoteˇ nejvetších vlastních cˇ ísel a posloupnost ~x (2k ) + λ1~x (2k ) k vlastnímu vektoru pˇríslušejícímu k vlastnímu cˇ íslu λ1 a posloupnost ~x (2k ) − λ1~x (2k ) k vlastnímu vektoru pˇríslušejícímu k vlastnímu cˇ íslu −λ1 .

X

A

A

16 / 70

Tomáš Oberhuber

Mocninná metoda



Remark 11 K napoˇcítání samotného vlastního vektoru staˇcí konstruovat ˇ tzv. Krylovovu posloupnost ~x (0) , ~x (0) , 2~x (0) , . . . a delení cˇ íslem ρk není nutné, je dobré jen pro lepší numerickou stabilitu.

A

A


Shrnutí a otázky

Remark 12 ˇ λλ1i . Konvergenci Rychlost konvergence závisí na pomeru mužeme ˚ urychlit vhodným posunem spektra o cˇ íslo λ∗ úpravou puvodní ˚ matice na matici − λ∗ . Tím mužeme ˚ λi −λ∗ ∗ zmenšit podíl λ1 −λ∗ . K ρ pak musíme pˇriˇcíst λ .

A

I

17 / 70

Tomáš Oberhuber

Mocninná metoda



Remark 13 Pokud chceme najít v absolutní hodnoteˇ nejmenší vlastní cˇ íslo matice , mužeme ˚ použít mocninnou metodu na −1 matici . Nejmenší vlastní cˇ íslo pak bude 1ρ . Pokud chceme najít vlastní cˇ íslo matice , které je nejblíže ˇ vlastní cˇ íslo zvolenému cˇ íslu λ0 budeme hledat nejvetší matice −1 . − λ0

A

A

A

A

I


Hledané cˇ íslo pak získáme jako

1 ρ

+ λ0 .

Shrnutí a otázky

18 / 70



Shrnutí a otázky

Mocninná metoda Remark 14 V pˇrípadeˇ z pˇredchozí poznámky napoˇcítáváme ~y (k +1) =

A−1~x (k ).

Abychom se vyhnuli výpoˇctu inverze matice vztah na

A, pˇrevedeme

A~y (k +1) = ~x (k ), a rˇešíme soustavu lineárních rovnic. S výhodou lze využít iterativní metody nebot’: • pro k malé je zbyteˇcné rˇešit lineární soustavu s

vysokou pˇresností a staˇcí proto poˇcítat menší poˇcet iterací • pro k velké se po sobeˇ jdoucí vektory ~ y (k ) a ~y (k +1) pˇríliš neliší a ~y (k ) je dobrým odhadem rˇešení pro soustavu s neznámým vektorem ~y (k +1) ˇ inverze matice, jež Tento trik, jak se zbavit výpoctu ˇ pouze aplikujeme na jeden nebo nekolik vektoru, ˚ se v ˇ numerice používá velmi casto. 19 / 70

Tomáš Oberhuber

Redukˇcní metoda



• tato metoda slouží k odseparování již nalezených

vlastních cˇ ísel, pokud chceme napˇr. napoˇcítat více než ˇ vlastní cˇ íslo jedno v absolutní hodnoteˇ nejvetší • nehodí se ale pro výpoˇcet kompletního spektra, kvuli ˚

rychlé ztráteˇ pˇresnosti


Shrnutí a otázky

Remark 15

A C

A

ˇ Bud’ tedy ∈ n,n a λ1 ∈ σ ( ) a ~x k nemu pˇríslušný vlastní vektor. Odvodíme metodu pro výpoˇcet matice ∈ n−1,n−1 , která má všechna vlastní cˇ ísla stejná jako až na λ1 , u kterého se sníží násobnost.

B C

A

20 / 70

Tomáš Oberhuber

Redukˇcní metoda



• matici

• v této bázi bude matice vypadat takto

LR algoritmus

P−1AP =

QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace QR iterace

A

pˇrevedeme do báze, tvoˇrené vektory ~x , ~e2 , . . . , ~en

kde

λ1 ~q t θ~

B

,

P je matice pˇrechodu mezi bázemi.


Shrnutí a otázky

21 / 70


Redukˇcní metoda

Tvar matice

P je zˇrejmý, její sloupce tvoˇrí bazické vektory 


P=   



Shrnutí a otázky

P



x1 x2 1 .. .. . .

   

xn

1

ˇ Inverzní matici k matici je taková, která první rˇádek podelí x1 a vynuluje všechny prvky pod diagonálou v prvním sloupci. Takovou matici ale snadno najdeme, známe ji z GEM.  

P

−1

  =  

1 x1 − xx12

.. . − xxn1

1 ..

. 1

     22 / 70


Redukˇcní metoda Roznásobením snadno získáme tvar matice



Shrnutí a otázky



B=   

B

a22 − xx21 a12 a23 − xx21 a13 . . . a2n − xx12 a1n a32 − xx31 a12 a33 − xx31 a13 . . . a3n − xx13 a1n .. .. .. . . . an2 − xxn1 a12 an3 − xxn1 a13 . . . ann − xxn1 a1n

   . 

ˇ Její tvar lze také snadno odvodit, pokud si uvedomíme, že ˇ pravý dolní cˇ tverec matice o rozmerech (n − 1) × (n − 1) odpovídá pˇresneˇ prvkum ˚ matice . Pˇri −1 násobení maticí zleva pak v tomto cˇ tverci od každého ˇrádku i odeˇcítáme xx1i násobek prvního ˇrádku matice . To je práveˇ význam matice −1 .

AP

P

P

A

A

23 / 70

Tomáš Oberhuber

Redukˇcní metoda



Podobneˇ lze odvodit, že ~q T =

B

1 (a12 , a13 , . . . , a1n ) x1

ˇ Nyní známe matici a mužeme ˚ najít nekteré její vlastní ˇ cˇ íslo λ2 . Bud’ ~z = (z2 , . . . , zn )T k nemu pˇríslušný vlastní vektor. Jak najít odpovídající vlastní vektor matice ?

A


Shrnutí a otázky

24 / 70

Tomáš Oberhuber Vlastní cˇ ísla matic ˇ Cásteˇ cný problém vlastních cˇ ísel

Redukˇcní metoda Logicky jako lineární kombinaci bazických vektoru˚ daných sloupci matice . K tomu ale musíme dopoˇcítat složku z1 . Musí platit

P



λ1 ~q T θ~

Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace

B

z1 ~z

λ1 z1 + ~q T ~z = ~z λ 2 z1 = λ2~z

B


=

λ1 z1 + ~q T ~z λ2~z

a tedy


λ1 z1 + ~q T ~z = λ2 z1 ⇒ z1 =

~q T ~z . λ2 − λ1

Shrnutí a otázky

25 / 70



Redukˇcní metoda Pokud by bylo λ1 = λ2 , pak musí být ~q T ~z = 0 a z1 lze volit ˇ libovolne. Vlastní vektor matice je pak dán jako 0 z1 , ~y = . = z1~x + ~z ~z

A

P


Shrnutí a otázky

26 / 70




A C

Bud’ ∈ n,n , rˇešíme úlohu nalezení všech vlastních cˇ ísel a vlastních vektoru˚ této matice. Trojúhelníková metoda • konstruujeme dveˇ posloupnosti (k ) ∞ • ∈ n,n horní trojúhelníkové s jedniˇckami na k =1 diagonále (k ) ∞ • ∈ n,n dolní trojúhelníkové k =1

L R

• •

C C

L(0) volíme libovolneˇ (napˇr. I) L(1) a R(1) získáme rozkladem matice AL(0) tj. platí L(1)R(1) = AL(0) • •

L(1) je dolní trojúhelníková s jedniˇckami na diagonále R(1) je horní trojúhelníková

• obecneˇ poˇcítáme iterace tvaru

L(k +1)R(k +1) = AL(k )

Shrnutí a otázky

• •

L(k(k+1) je dolní trojúhelníková s jedniˇckami na diagonále R +1) je horní trojúhelníková 27 / 70

Tomáš Oberhuber




Remark 16

L(k ) → L a R(k ) → R, pak platí AL = LR ⇒ A = LRL−1

Pokud ukážeme, že

a jde tedy a podobnostní transformaci a vlastní cˇ ísla matice

A najdeme na diagonále matice R.


Shrnutí a otázky

28 / 70


Trojúhelníková metoda Remark 17 ˇ Jak vypoˇcítat vlastní vektory? Rešíme nejprve rovnici

R − λI) ~y i = 0.



(

R

I

Z tvaru matice − λ (dolní trojúhelníková s 0 na diagonále na i-tém rˇádku) snadno vidíme, že   0 j i k k =i rjj−λ nebot’ pro j > i musí platit


j−1 X


Shrnutí a otázky

yki rjk + yjj rjj − λ = 0.

k =i

R

A

Jelikož matice je matice vyjádˇrená v bázi dané sloupci matice , platí, že vlastní vektory matice jsou

L

A

L

~x i = ~y i . 29 / 70

Tomáš Oberhuber




•

L

ˇ nemusí být ani dolní 0 lze volit libovolne, trojúhelníková • díky tomu má metoda samoopravující schopnost • pokud bychom nekdy ˇ ˇ lze ho brát napoˇcítali (k ) špatne,

jako nové

L(0)

L


Shrnutí a otázky

30 / 70



Existence LU rozkladu Remark 18 Obeˇ metody jsou postaveny na poˇcítání LU rozkladu. Ten ale existuje jen pro silneˇ regulární matice. V první iteraci ˇ ˇ musíme pˇredpokládat, že toto je splneno. Pˇripomenme si navíc vzorce pro výpoˇcet LU rozkladu:

Trojúhelníková metoda LR algoritmus QR algoritmus

= aij −

lij


QR iterace

= (aij −

uij

lik ukj )/lii proi < j

k =1


Shrnutí a otázky

lik ukj proj ≤ i

k =1 i−1 X

Householderovy transformace Givensovy rotace

j−1 X

L U

Z nich vidíme, že prvky matic a závisí spojiteˇ na ˇ prvcích matice . Pokud tedy matici zmeníme jen málo, ˇ existovat. bude LU rozklad opet

A

A

31 / 70

Tomáš Oberhuber

Existence LU rozkladu



Theorem 19

A

Necht’ matice je silneˇ regulární tj. existuje její LU rozklad. Bud’ matice taková, že k k je dostateˇcneˇ malé. Pak existuje i LU rozklad matice + .

E

E A E

Theorem 20

Necht’ A = I + E, kde kEk je malé. Potom existuje rozklad A = LR, kde L je dolní trojúhelníková s jedniˇckami na diagonále a R je horní trojúhelníková. Platí-li, že pokud kEk → 0, pak L → I a R → I.

Shrnutí a otázky

32 / 70

Tomáš Oberhuber

Konvergence trojúhelníkové metody



Theorem 21 Pokud existuje trojúhelníkový rozklad matice k (0) = L(k ) R(k ) , pak platí

AL

L(k ) =


R(k ) =

LR algoritmus QR algoritmus

L(k ), R(k )R(k −1) . . . R(1).

Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace

Remark 22

Givensovy rotace

V dukazu ˚ konvergence budeme zkoumat matici k (0) . Odvodíme, kdy existuje její LU rozklad k (0) = L(k ) R(k ) . Jelikož platí, že (k ) = L(k ) , dokážeme-li, že L(k ) → , bude platit i (k ) → . Pak již lze i dokázat, že R(k ) → , cˇ ímž máme vyšetˇrenou konvergenci trojúhelníkové metody.


Shrnutí a otázky

L

L

AL

L

AL L

R

33 / 70

Tomáš Oberhuber

Konvergence trojúhelníkové metody



Theorem 23

A C AL

Necht’ matice ∈ n,n je regulární a diagonalizovatelná a pˇredpokládejme, že pro dostateˇcneˇ velké k existují LU (k ) a necht’ dále platí pˇredpoklady, které rozklady matic 1 . Pak posloupnosti matic (k ) a (k ) z vyplynou z dukazu ˚ trojúhelníkové metody konvergují a na diagonále matice je spektrum matice seˇrazené sestupneˇ podle velikosti v ˇ absolutní hodnote.

L

R

A

R


Shrnutí a otázky

Dukaz. ˚ Pouze pro pˇrípad všech vlastních cˇ ísel jednonásobných a ˇ ruzných ˚ v absolutní hodnote. 1

Pro pˇrípad všech vlastních cˇ ísel jednonásobných a ruzných ˚ v absolutní hodnoteˇ je to existence LU rozkladu matic a −1 (0)

X X L

34 / 70



LR algoritmus LR algoritmus • konstruujeme tˇri posloupnosti • •

A(k )o∞k∞=1 ∈ Cn,n Lˆ (k ) k =1 ∈ Cn,n dolní trojúhelníkové s jedniˇckami na

n

diagonále n o∞ ˆ (k ) • ∈

R

k =1

Cn,n horní trojúhelníkové

• na poˇcátku volíme


Shrnutí a otázky

A(1) Lˆ (1)Rˆ (1) A(2)

= = =

A, A(1), Rˆ (1)Lˆ (1)

• obecneˇ mají jednotlivé iterace tvar

Lˆ (k )Rˆ (k ) A(k +1)

= =

A(k ), Rˆ (k )Lˆ (k ). 35 / 70

Tomáš Oberhuber

LR algoritmus

Vlastní cˇ ísla matic ˇ Cásteˇ cný problém vlastních cˇ ísel Mocninná metoda

Remark 24 Z úvahy

Redukˇcní metoda

A(k +1)


R L L

ˆ (k ) ˆ (k ) = −1 (k ) (k ) ˆ ˆ ˆ (k ) ˆ (k ) =

LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces

=


Shrnutí a otázky

Lˆ (k )

L

−1

R L

A(k )Lˆ (k ).

A

(k ) ∞ plyne, že pokud posloupnost konverguje k horní k =1 trojúhelníkové matici, budou na její diagonále vlastní cˇ ísla matice .

A

36 / 70

Tomáš Oberhuber

LR algoritmus



Shrnutí a otázky

• má menší nároky na pamet’ ˇ než trojúhelníková metoda • trojúhelníková metoda potˇrebuje pracovní pole pro (k ) napoˇcítání LU rozkladu, pro výsledek souˇcinu a pro matici • LR algoritmus potˇrebuje pracovní pole pro LU rozklad a napoˇcítání souˇcinu ˆ (k ) ˆ (k ) • LR algoritmus nemá tak dobrou samoopravující schopnost, takže vlivem numerických chyb mužeme ˚ získat špatné spektrum

AL

A

R L

A

• nepamatuje si matici • nekonverguje pro libovolnou poˇcáteˇcní matici, tou je zde matice na rozdíl od trojúhelníkové metody, kde ji byla libovolná matice 0

A

L

37 / 70

Tomáš Oberhuber

Konvergence LR algoritmu



Theorem 25 Existuje-li trojúhelníkový rozklad matice platí

Ak = L(k )R(k ), pak

Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace

L(k ) = R(k ) =

Lˆ 1Lˆ 2 . . . Lˆ (k ) Rˆ (k )Rˆ (k −1) . . . Rˆ (1).


Shrnutí a otázky

38 / 70

Tomáš Oberhuber




Remark 26 Pokud v trojúhelníkové metodeˇ volíme k (0) = k a musí platit

AL

L(k ) R(k ) = R(k ) R(k −1) . . . R(1) L(k ) =


Shrnutí a otázky

A

= =

L(0) = I, pak je

Lˆ 1Lˆ 2 . . . Lˆ (k ), Rˆ (k )Rˆ (k −1) . . . Rˆ (1),

tj.

L(k ) R(k )

= =

Lˆ 1Lˆ 2 . . . Lˆ (k ), Rˆ (k ).

39 / 70

Tomáš Oberhuber




Remark 27 Z konvergence trojúhelníkové metody (k ) → , plyne

R

A(k ) = Lˆ (k )Rˆ (k ) = L(k−1) L(k )R(k ) → L−1LR = R,


Shrnutí a otázky

R

L(k ) → L a

−1

A

tj. matice (k ) z LR algoritmu skuteˇcneˇ konverguje k horní trojúhelníkové matici. Na diagonále má vlastní cˇ ísla matice seˇrazené sestupneˇ podle velikosti v absolutní hodnoteˇ (to už jsme ukázali).

A

40 / 70

Tomáš Oberhuber



ˇ Nyní již mužeme ˚ vše shrnout ve formeˇ následující vety:

Redukˇcní metoda

Výpoˇcet kompletního spektra matice Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace QR iterace

Theorem 28

A C

Necht’ matice ∈ n,n je regulární a diagonalizovatelná a pˇredpokládejme, že trojúhelníková metoda konverguje s volbou 0 = . Pak LR algoritmus konverguje také, a platí, že posloupnost (k ) konverguje k horní trojúhelníkové matici s hlavními cˇ ísly matice na diagonále seˇrazenými ˇ sestupneˇ podle velikosti v absolutní hodnote.

L

I

A

A


Shrnutí a otázky

41 / 70

Tomáš Oberhuber

QR algoritmus



Remark 29 Velkou nevýhodou LR algoritmu je jeho špatná numerická stabilita zejména pˇri aplikování na velké matice. Proto byl v roce 1961 J.G.F. Francisem navržený QR algoritmus. Funguje stejneˇ jako LR algoritmus, ale místo LU rozklad napoˇcítává QR rozklad, tj. rozklad na unitární a horní trojúhelníkovou matici. Omezíme se nyní pouze na reálné matice.


Shrnutí a otázky

42 / 70

Tomáš Oberhuber

QR algoritmus



Theorem 30

A R A QR

Bud’ ∈ n,n regulární matice. Pak existuje rozklad = , kde matice je unitární a je horní trojúhelníková. Pokud budeme pˇredpokládat, že diagonální prvky matice jsou kladné, pak je tento rozklad jednoznaˇcný.

Q

R

R


Shrnutí a otázky

43 / 70

Tomáš Oberhuber

Výpoˇcet QR rozkladu



Existují tˇri zpusoby, ˚ jak napoˇcítat QR rozklad: • Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces • Householderovy transformace • Givensovy rotace


Shrnutí a otázky

44 / 70



Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces • jde o algoritmus, který pˇrevede lineárneˇ nezávislé

vektory ~x (1) , . . . , ~x (n) na množinu vektoru˚ ~q (1) , . . . , ~q (n) , které jsou vzájemneˇ ortonormální a tvoˇrí bázi stejného prostoru jako puvodní ˚ vektory


Shrnutí a otázky

45 / 70

Tomáš Oberhuber



1: 2: 3:

Redukˇcní metoda


4: 5: 6: 7:


Shrnutí a otázky

8: 9: 10: 11: 12:

procedure G RAMM S CHMIDT K LASICKÝ(~x (1) , . . . , ~x (n) ) (n) ~ ~r (1) := ~ r (2) :=

. . . := ~r := 0 (1) (1)

~ r1 := x 2 ~q (1) := r1 ~x (1) 11 for k = 2, . . . , n do ˜ (k ) := ~x (k ) ~q for i = 1, . . . , k − 1 do (k ) ri := ~q (i) , ~x (k ) ˜ (k ) := ~q ˜ (k ) − r (k )~q (i) ~q i end for

˜ (k ) (k ) rk := ~q

2 1 ˜ (k (k ) ~q ) ~q := (k )

rk

end for return [~q (1) , . . . , ~q (n) ], [~r (1) , . . . , ~r (n) ] 15: end procedure 13:

14:

46 / 70

Tomáš Oberhuber




• pro lepší numerickou stabilitu se používá modifikovaný

G.-S. ort. proces procedure G RAMM S CHMIDT M ODIFIKOVANÝ(~x (1) , . . . , ~x (n) ) (n) ~ ~r (1) := ~ 2: r (2) :=

. . . := ~r := 0 (1) 3: r := ~x (1) 1:

1

4: 5: 6: 7: 8:


Shrnutí a otázky

9: 10: 11: 12:

2

~q (1) := r1 ~x (1) 11 for k = 2, . . . , n do ˜ (k ) := ~x (k ) ~q for i = 1, . . . , k − 1 do (k ) ˜ (k ) ri := ~q (i) , ~q ˜ (k ) := ~q ˜ (k ) − r (k )~q (i) ~q i

end for

˜ (k ) (k ) rk := ~q

2 ˜ (k ) ~q (k ) := 1 ~q (k )

rk

end for 14: return [~q (1) , . . . , ~q (n) ], [~r (1) , . . . , ~r (n) ] 15: end procedure 13:

47 / 70

Tomáš Oberhuber




Shrnutí a otázky

Theorem 31

A R

Bud’ ∈ n,n regulární matice. Položme v ˇ Gramove-Schmidtov eˇ ortonormalizaˇcním procesu vektory ~x (1) , . . . , ~x (n) rovny jednotlivým sloupcum ˚ matice , tj. ~x (1) = a·1 , . . . , ~x (n) = a·n . Pak lze psát   r11 r12 . . . r1n  r22 . . . r2n    (a·1 , . . . , a·n ) = ~q (1) , . . . , ~q (n)  ..  , . .  . .  rnn

A

A QR

Q

R

tj. = , kde je unitární a je horní trojúhelníková matice (s kladnými prvky na diagonále).

48 / 70

Tomáš Oberhuber



Pˇripomeneme si:



Definition 32 ˇ maticí (elementární Householderovou reflekcní unitární ~ w tvaru maticí) nazveme každou matici ~ = − 2w ~w ~ ∗, w

H

H

I

~ jeq kde w Householderuv ˚ vektor, pro který platí

w

~ ,w ~ = 1. ~ 2= w

Shrnutí a otázky

49 / 70

Tomáš Oberhuber




Theorem 33

H


Shrnutí a otázky

H

C

QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces

C

~ je Householderova reflekˇcní matice a ~v ∈ n w Necht’ ~ ~v je zrcadlový obraz je libovolný vektor. Pak vektor w vektoru ~v podle nadroviny o n ~ H ~x = w ~ , ~x = 0 L ≡ ~x ∈ n | w ˇ v tom smyslu, že splnuje následující podmínky:

~ ~v = ~v • w ~ ~v + ~v ∈ L • w ~ ~v − ~v ⊥ L. • w

H H H

50 / 70

Tomáš Oberhuber




Remark 34 Chceme-li pomocí Householderovy transformace

transformovat vektor ~x na vektor ~y , kde ~x 2 = ~y 2 , pak volíme ~x − ~y

~ = w

~x − ~y 2

(1)

Pokud zvolíme ~y = ± ~x 2 ~e , pak získáme unitární transformaci, která nuluje všechny složky vektoru ~x kromeˇ první.

Shrnutí a otázky

51 / 70

Tomáš Oberhuber




Remark 35 Pokud by byla první složka vektoru ~x silneˇ pˇrevládající tj.

|x1 | ≈ ~x 2 , pak bychom pˇri špatné volbeˇ znaménka vektoru ~y mohli ˇ témeˇ ˇ r nulou. Proto volíme dostat ~x − ~y velmi malé a delit

~y = −signx1 ~x 2 ~e(1)


Shrnutí a otázky

tj.

~x + signx1 ~x 2 ~e1

~ = w

~x + signx1 ~x ~e1 . 2 2

52 / 70



Householderovy transformace • obecne, ˇ pokud pro k ≥ 1 chceme zachovat prvních

k − 1 složek vektoru ~x a nulovat složky k + 1, . . . n, použijeme matici (k ) θ (k ) = ˜ (k ) , θ

Q

Q ˜ (k ) ∈ Rn−k +1,n−k +1 , I(k ) ∈ Rk −1,k −1 a kde Q T Q˜ (k ) = I − 2w~ (k ) w~ (k ) ,


Shrnutí a otázky

I

~ (k ) w

R

~x (k ) + signx1(k ) ~x (k ) 2 ~e(k )

, =

~ (k )

(k )

x + signx1 ~x (k ) 2 ~e(k ) 2

kde ~x (k ) ∈ n−k +1 je vektor tvoˇrený posledními n − k + 1 složkami vektoru ~x a vektor ~e(k ) je první vektor standardní báze prostoru n−k +1,n−k +1 .

R

53 / 70

Tomáš Oberhuber



• pro vektor ~ y=



yj =

  

Q(k )~x pak platí,

x

j (k ) (k )

signx1 ~x 2

pro pro 0 pro

j = 1, . . . , k − 1 j =k j = k + 1, . . . n


Shrnutí a otázky

• vhodnou aplikací Householderových transformací na

sloupce matice pod diagonálou

A tak lze eliminovat nenulové prvky

• jelikož se vše deje ˇ pomocí unitárních transformací,

získáme unitarní pˇrevod na matici v horním trojúhelníkovém tvaru, tj. QR rozklad

54 / 70



Householderovy transformace 1: 2: 3: 4: 5: 6:

LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace

7: 8:


Shrnutí a otázky

9: 10: 11: 12: 13:

A

procedure H OUSEHOLDER( ) (1) = (1) = for k = 1, . . . , n − 1 do ) ~x (k ) := (k k ...n,k (= složky k . . . n k -tého sloupce ) ~x (k ) +signx1(k ) k~x (k ) k ~e(k ) n−k +1 2

∈ ~ (k ) :=

w

(k )

~x (k ) +signx1 k~x (k ) k2~e(k ) 2 T ˜ (k ) := − 2w ~ (k ) w ~ (k ) ∈ n−k +1,n−k +1 (k ) (k ) := ∈ n,n ˜ (k )

R Q

A I

R

Q I I Q Q Q(k +1) := Q(k(k))Q(k ) R(k +1) := Q R(k ) end for return Q(n) , R(n)

R

R

R

end procedure

55 / 70



Givensovy rotace • Givensovy rotace jsou ortogonální transformace, které

R

ˇ eliminovat jednotlivé prvky vektoru ~x ∈ n umožnují • pro dvojici indexu˚ i, j a úhel θ je Givensova rotace definována jako  1  1   ..  .   cos θ sin θ   . .. (i, j, θ) =    − sin θ cos θ   ..  .   1

                

G

1

Shrnutí a otázky

• pro vektor ~ x∈

Rn odpovídá souˇcin ~y = G (i, j, θ) ~x

ˇ hodinových rotaci složek (xi , xj ) o úhel θ po smeru ruˇciˇcek 56 / 70


Givensovy rotace • platí tedy


  xj cx + sxj yk =  i −sxi + cxj



Shrnutí a otázky

pro k 6= i, j pro k = i pro k = j

kde jsme použili znaˇcení s = sin θ a c = cos θ. • volbou

xj xi c=q , s=q , xi2 + xj2 xi2 + xj2 q pak bude yi = xi2 + xj2 a yj = 0

• to odpovídá rotaci o úhel θ = arctan −xj /xi • pomocí Givensových rotací lze postupneˇ eliminovat

ˇ tak získat QR všechny prvky pod diagonálou a opet rozklad 57 / 70

Tomáš Oberhuber

QR rozklad


Ukázali jsme si tˇri zpusoby, ˚ jak získat QR rozklad

Redukˇcní metoda


• Gramuv-Schmidt ˚ uv ˚ ortogonalizaˇcní proces • je numericky nestabilní • používá se v modifikované podobeˇ v nekterých ˇ metodách pro ˇrešení soustav lineárních rovnic • Householderovy a Givensovy transformace jsou

ˇ numericky stabilnejší • ve všech pˇrípadech je složitost n3


Shrnutí a otázky

58 / 70

Tomáš Oberhuber

QR rozklad


• ztrátu ortogonality lze pomeˇ ˇ rovat pomocí hodnoty

I QQ∗k

k −



• lze ukázat, že platí ( oznaˇcuje strojovou pˇresnost

aritmetiky)


Shrnutí a otázky

Algoritmus Householderuv ˚ QR rozklad Givensuv ˚ QR rozklad Klasický GS Modifikovaný GS Iteraˇcní GS

I QQ∗k

k −

A A

κ( )2 κ( )

59 / 70



QR algoritmus QR algoritmus • pro zadanou matici posloupnosti •

• •

A ∈ Rn,n konstruujeme tˇri

(k ) ∞ T ∈ Rn,n (k ) k∞=1 Q k =1 ∈ Rn,n ortogonální matice (k ) ∞ R k =1 ∈ Rn,n horní trojúhelníkové

• na poˇcátku volíme


Shrnutí a otázky

T(1) Q(1)R(1) T(2)

= = =

A, T(1), R(1)Q(1)

• obecneˇ mají jednotlivé iterace tvar

Q(k )R(k ) T(k +1)

= =

T(k ), R(k )Q(k ). 60 / 70

Tomáš Oberhuber

QR algoritmus


• platí tedy

R(k ) = Q(k )



a

QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces Householderovy transformace Givensovy rotace

Shrnutí a otázky

−1

T(k )Q(k ),

resp.

T(k ) = Q(0)Q(1) . . . Q(k ) A Q(0)Q(1) . . . Q(k )


T(k )

T(k +1) = R(k )Q(k ) = Q(k )


−1

pro k > 0 a tedy matice matici .

A

T

T(k ) je ortogonálneˇ podobná

61 / 70

Tomáš Oberhuber

Konvergence QR algoritmu



Remark 36 ˇ spojitou funkcí prvku˚ matice QR rozklad je opet ˇ napˇr. z Grammova-Schmidtova algoritmu. videt

Lemma 37 Existuje-li QR rozklad matice

A, jak lze

Ak = Q(k )R(k ), pak platí

LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces

Q(k ) =


R(k ) =


Q(1)Q(2) . . . Q(k ), R(k )R(k −1) . . . R(1).


Shrnutí a otázky

Dukaz. ˚ Stejný jako u LR algoritmu.

62 / 70

Tomáš Oberhuber

Konvergence QR algoritmu



Theorem 38 Je-li matice

ˇ A ∈ Rn,n taková, že její vlastní cˇ ísla λi splnují |λ1 | > |λ2 | > . . . > |λn | ,


T


pak posloupnost (k ) konverguje k horní trojúhelníkové matici s vlastními cˇ ísly λi na diagonále seˇrazené podle velikosti. Je-li symetrická, pak (k ) konverguje k diagonální matici.

A

T

Shrnutí a otázky

63 / 70

Tomáš Oberhuber

QR algoritmus



• klasický pˇrístup ke QR rozkladu vyžaduje v každém

kroku n3 operací • ukážeme si tzv. Hessenbergovy QR iterace se

složitostí n2


Shrnutí a otázky

64 / 70

Tomáš Oberhuber

Hessenberguv ˚ tvar matice



Definition 39 Matice v Hessenbergoveˇ tvaru je horní trojúhelníková matice, která má navíc nenulové prvky bezprostˇredneˇ pod diagonálou.


Remark 40 Hessenberguv ˚ tvar je zˇrejmeˇ nejjednodušší tvar, na který dokážeme libovolnou regulární matici pˇrevést podobnostní transformací, kterou lze získat pˇrímým algoritmem.

Shrnutí a otázky

65 / 70

Tomáš Oberhuber




Hessenbergovy QR iterace • metoda využívá Hessenberguv ˚ tvar matice • na zaˇcátku metody pˇrevedeme matici

A

podobnostní transformací do Hessenbergova tvaru se složitost O(n3 )

• ukážeme, že následující iterace lze provádet ˇ se

složitostí O(n2 )


Shrnutí a otázky

66 / 70



Shrnutí a otázky

Hessenbergovy QR iterace • v prvním kroku nulujeme prvky v prvním sloupci na

ˇrádcích 3 až n pomocí Householderovy transformace 1 (1) := ˜ (1)

Q

Q

A∗ = Q(1)A dostaneme matici, jejíž první sloupec odpovídá Hessenbergovu tvaru, není ale podobná matici A (1) • trik je v tom, že matice Q pˇri násobení ˇ první sloupce matice A∗ , A(1) = A∗(Q(1))T nemení takže ten stále odpovídá Hessenbergovu tvaru a A(1) je navíc podobná matici A • transformací

• takto po celkem n − 2 krocích získáme matici v

Hessenbergoveˇ tvaru, kterou oznaˇcíme podobná puvodní ˚ matici

A

H(0), a která je 67 / 70



Hessenbergovy QR iterace • Hessenbergovy QR iterace generují posloupnost matic

H(k ), které jsou v Hessenbergoveˇ tvaru • pro QR rozklad matice H(k −1) staˇcí eliminovat pouze n − 1 nenulových prvku˚ pod diagonálou

• k tomu staˇcí provést n − 1 Givensových rotací • definujme

Trojúhelníková metoda LR algoritmus QR algoritmus Gramuv-Schmidt ˚ uv ˚ ortonormalizaˇcní proces

Q(k )

• pak je


T

=

G(kn−1) . . . G(k1 )

Q(k ) H(k −1) = R(k ) T


Shrnutí a otázky

• a tudíž

H(k −1) = Q(k )R(k )

• a z definice QR iterací je

H(k ) = R(k )Q(k ) = R(k ) G(k1 )

T

...

G(kn−1)

T 68 / 70

Tomáš Oberhuber



Lemma 41 Výsledkem násobení


R(k ) G(k1 )


T

...

G(kn−1)

T


ˇ matice v Hessenbergovu tvaru. je opet • celou iteraci se nám podaˇrilo provést s O(n2 )

operacemi • na konci mámeˇ opet ˇ matici v Hessenbergovu tvaru

Shrnutí a otázky

69 / 70



Shrnutí a otázky

Shrnutí a otázky • cˇ ásteˇcný problém vlastních cˇ ísel • mocninná metoda • jak v praxi zaruˇcit její podmínky konvergence? • jak napoˇcítat nejmenší vlastní cˇ íslo? • redukˇcní metoda • úplný problém vlastních cˇ ísel • trojúhelníková metoda • LR algoritmus • QR algoritmus • • • •

ˇ proces Grammuv-Schmidt ˚ uv ˚ ortogonalizacní Householderovy transformace Givensovy rotace Hessenbergovy iterace

• na cem ˇ závisí rychlost konvergence? • jak konvergenci urychlit? 70 / 70

Faculty of Nuclear Sciences and Physical Engineering Czech Technical University in Prague

Recommend Documents