PA081: Programování numerických výpočtů
7. Rozklady matic
Aleš Křenek
Přehled a
motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Dekompozice na singulární hodnoty
Vlastní hodnoty a vektory
jaro 2013
Motivace
danou matici A vyjádříme jako součin
PA081: Programování numerických
výpočtů
A. Křenek
Přehled a motivace
A = Ai A2 ... A,
LU
dekompozice
► matice Ai, A2,... An mají nějaké speciální vlastnosti
► vlastnosti rozkladu lze prakticky využít
► řešení systému lineárních rovnic
► existují implementace algoritmů se známými vlastnostmi
► zpravidla numericky stabilní
► optimalizované na konkrétní CPU, maximálně efektivní
zřejmější vlastnosti problému popsaného původní A
► kritické body vektorového pole
► statisticky významné vlastnosti nějakého jevu
► podmíněnost systému lineárních rovnic
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
QR
dekompozice
■O0.O
Přehled
► A
LU, resp. PA = LU
► L a U jsou dolní a horní trojúhelníková
► P je permutace řádků Choleského: A = LLr
► pro symetrickou pozitivně definitní A A = QR
► Q ortogonální, R horní trojúhelníková
► symetrické varianty RQ, QL a LQ singulární hodnoty: A = USV
► U, V ortogonální, 2 diagonální spektrální: A = VAV-1
► A diagonální, vlastní hodnoty
► varianta Jordánova: blokově diagonální A, násobné vlastní hodnoty
Shurova: A = VSVr
► V ortogonální
► S horní trojúhelníková s vlastními hodnotami A na diagonále
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
3/30
LU dekompozice
Princip
předpokládejme, že dokážeme rozložit A = LU tak, že
► L je spodní trojúhelníková matice (prvky nad diagonálou jsou nulové)
► U je horní trojúhelníková matice (prvky pod diagonálou jsou nulové)
potom lze psát
Ax = (LU)x = L(Ux) = b původní systém lze vyřešit postupným řešením
Ly = b     a     Ux = y to je triviální dopřednou a zpětnou substitucí
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
4/30
LU dekompozice
Výpočet rozkladu
prvky rozkladu A = LU lze, s ohledem na nuly psát
kj: aij í = j: aij i> j: aij
unhj + ui2hj unhj + ui2hj unhj + ui2hj
UaXij uidjj
je to systém N2 rovnic v N2 + N neznámých
► diagonála je pokryta u i l
► můžeme tedy volit la = 1
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
5/30
LU dekompozice
Croutův algoritmus
► postupně pro j = 1,2,...,N: ► pro i = 1,2,... j spočítáme na základě uvedených rovnic
i-1
Uij — Clij — ^ lik^-kj fc=l
► pro i = j + l,j + 2,
1
JJ \ k=l
postup vždy využívá dříve spočtené prvky
k numerické stabilitě je třeba navíc dodat pivoting Ijj
Vlastní hodnoty a vektory
■O0.O
LU dekompozice
Shrnutí a použití
řešení lineárních rovnic pro libovolný počet b
► Croutův algoritmus O (N3)
► dopředná a zpětná substituce O (N2)
► všechna b není třeba znát dopředu - hlavní přednost metody
výpočet inverzní matice
► řešení systému pro b = bázové vektory
► potřebujeme-li počítat A_1B, je lepší přímo pro sloupce B než explicitní vyjádření A-1
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Choleského dekompozice
předpokládáme A po prvcích platí
LL
la
a i
i-1 fc=0
ik
podobně jako při LU dekompozici vždy využíváme dříve spočtené prvky
odmocninu lze vždy spočítat pro symetrickou pozitivně definitní A
► omezenější, ale stále významná třída problémů algoritmus je efektivnější a numericky stabilní
► cca. 2x méně operací než LU, menší paměťová náročnost
► má smysl používat speciální funkce z knihoven
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
8/30
QR dekompozice
► rozklad A = QR
► Q ortogonální, R trojúhelníková
► systém Ax = b lze psát
Rx = Qrb
► jednoduché řešení substitucí
► lepší numerické vlastnosti
► metody konstrukce - nulování prvků pod diagonálou
► Gram-Schmidtův proces
► Householderovy transformace
► Givensovy rotace
PA081: Programování numerických
výpočtů
A. Křenek
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Dekompozice na singulární hodnoty
Vlastní hodnoty a vektory
QR dekompozice
Gram-Schmidtův proces
Gram-Schmidtův ortogonalizační proces
Ul
ui = ai vi
ui I
u2 = a2 - projVla2 v2 =
I u21
u3 = a3 - proj a3 - projv a3        v3 = —-
IIU3I
potom
Q = (vi... v„
R
/ vi ■ ai   vi ■ a2 0     v2 ■ a2
numerický problém, jsou-li některé au a j skoro kolmé
QR dekompozice
Householderova transformace
► zrcadlení daného vektoru podle nadroviny
► zarovná s bázovým vektorem, zachová velikost
► pro libovolné x:
U T
u = x + cxei     v = -—-     Q = I - 2w
l|u||
► potom Qx = (cx, 0,..., 0)r
QR dekompozice
Householderova transformace
zrcadlení daného vektoru podle nadroviny
► zarovná s bázovým vektorem, zachová velikost pro libovolné x:
u = x + cxei v
potom Qx = (a, 0,..., 0)r triangulace A
í-2wJ
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
R = an !... QiA     a tedy
Q=OjoJ...Oj_1
Základní dekompozice - shrnutí
LU
► ve variantě PA = LU existuje pro všechny čtvercové matice
► odpovídá Gaussově eliminaci - trpí stejnými numerickými problémy
Choleského
čtvercové, symetrické, pozitivně dennitní matice
► numericky stabilní (to je ale LU pro tyto matice také)
► cca. 2 x méně operací
QR
► obecné m x n matice
► existují numericky stabilní konstrukce (Householderova transformace)
► výpočetně náročnější
použití především k řešení systémů lineárních rovnic
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
12/30
Dekompozice na singulární hodnoty
Základní tvrzení
libovolnou reálnou (i komplexní) matici lze rozložit
UM
vJVxJV
V
NxN>
U je sloupcově ortogonální
► až na nulové sloupce v případě M < N S diagonální a V ortogonální rozklad je unikátní až na
► současnou perumutaci sloupců všech tří matic
► lineární kombinaci sloupců U, V odpovídajících nulovým 07
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Geometrický význam
► A je složení transformací
► rotace/zrcadlení V-1
► zvětšení/zmenšení faktory 07 ve směrech e;, včetně degenerace (07 = 0)
► rotace/zrcadlení a projekce do méně/více dimenzí U
Dekompozice na singulární hodnoty
Geometrický význam
A je složení transformací
► rotace/zrcadlení V-1
► zvětšení/zmenšení faktory 07 ve směrech e;, včetně degenerace (07 = 0)
► rotace/zrcadlení a projekce do méně/více dimenzí U
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
obor hodnot zobrazení A
► sloupce U odpovídající nenulovým 07 jsou jeho generátory
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Geometrický význam
A je složení transformací
► rotace/zrcadlení V-1
► zvětšení/zmenšení faktory 07 ve směrech eu včetně degenerace (07 = 0)
► rotace/zrcadlení a projekce do méně/více dimenzí U
obor hodnot zobrazení A
► sloupce U odpovídající nenulovým 07 jsou jeho generátory
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
nulový prostor N(A) = {x e RN : Ax = 0}
► řádky VT odpovídající nulovým crr jsou jeho generátory
Dekompozice na singulární hodnoty
Numerický význam
„oddělení zrna od plev"
sloupce U a V jsou kolmé a normované
veškeré potenciální degenerace soustředěny do 2
► singularity A odpovídají nulovým 07
► včetně numerických (crr » 0)
numericky velmi stabilní algoritmus dekompozice lze použít na řešení systémů lineárních rovnic
► M < N a M = N singulární: reprezentant řešení + generátor prostoru
► M > N: nejbližší řešení
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
□     s   4 :
15/30
Dekompozice na singulární hodnoty
Použití pro M = N
řešení systému rovnic, resp. výpočet inverzní matice A"1 = V[diag(l/o-r)]Ur
kdy to nejde
► jedno nebo více 07 je nulových - A byla singulární *• Omin/Omax < e (špatně podmíněná matice) - standardní metody řešení selhaly
Dekompozice na singulární hodnoty
Použití pro M = N
► řešení systému rovnic, resp. výpočet inverzní matice
► kdy to nejde
► jedno nebo více 07 je nulových - A byla singulární *• Omin/Omax < e (špatně podmíněná matice) - standardní metody řešení selhaly
► označme
► rovnice Ax = b nemusí mít řešení, přesto zkusíme
x = vrurb
A"1 = V[diag(l/o-r)]Ur
Dekompozice na singulární hodnoty
Použití pro M = N
► hledáme nejbližší řešení, tj. minimalizujeme |Ax - b|
► pro libovolné x' je A(x + x') - b = Ax - b + b', kde b' =
|Ax-b + b'| = |(USVr)(VS'Urb) - b + b'| = |(USS'Ur - J)b + b'| = |U((SS' - J)Urb + Urb')| = |(ZZ' - J)Urb + Urb'|
► 55' - J je diagonální s nenulovými prvky pro 07 = 0
► b' je v oboru hodnot A, tedy Urb' má nenulové prvky právě pro 07 =/= 0
► minimum právě pro b' = 0 a tedy i x' = 0
Ax'
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Použití pro M = N
SVD solution of A x = d
Dekompozice na singulární hodnoty
Použití pro M = N - prakticky
► singularitu A detekujeme podle 07 = 0
► vypočteme nejbližší řešení jako x = VS'Urb
► dosazením ověříme, zda je to přesné řešení
► když ne, víme, že přesné řešení neexistuje
► máme nejbližší aproximaci
Dekompozice na singulární hodnoty
Použití pro M = N - prakticky
singularitu A detekujeme podle 07 = 0 vypočteme nejbližší řešení jako x = VS'Urb dosazením ověříme, zda je to přesné řešení
► když ne, víme, že přesné řešení neexistuje
► máme nejbližší aproximaci
špatně podmíněná matice | crmax | » | crmin |
► lépe v 2' vynulovat i taková 07
► paradoxní - zahazujeme část vstupní informace
► v praxi dává lepší výsledky - právě tento vstup má tendenci škodit
► stanovení prahu 07 ~ 0 není triviální
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
□     s   4 :
■OQ.O
19/30
Dekompozice na singulární hodnoty
Použití pro M * JV
méně rovnic, M < N nekonečně mnoho řešení
rozklad na singulární hodnoty - N - M nulových 07
► nemusí být přesně nulové (numerické nepřesnosti)
► může jich být více díky dalším singularitám
5/ vypočítáme vynulováním problematických 07 přímo vypočteme reprezentativní řešení x
► včetně ověření, zdaje skutečně řešením
sloupce V odpovídající nulovaným 07 generují prostor dalších řešení
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Použití pro M * JV
více rovnic, M > N
neexistuje přesné řešení, hledáme nejbližší aproximaci rozklad na singulární hodnoty
► obecně nemusí dát žádná nulová 07
► získáme nejbližší aproximaci řešení, viz naznačený důkaz
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Použití pro M * JV
více rovnic, M > N
neexistuje přesné řešení, hledáme nejbližší aproximaci rozklad na singulární hodnoty
► obecně nemusí dát žádná nulová 07
► získáme nejbližší aproximaci řešení, viz naznačený důkaz (skoro) nulové singulární hodnoty
► skrytá degenerace systému
► může vést na jedno nebo i více přesných řešení
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Použití pro M * JV
více rovnic, M > N
neexistuje přesné řešení, hledáme nejbližší aproximaci rozklad na singulární hodnoty
► obecně nemusí dát žádná nulová crr
► získáme nejbližší aproximaci řešení, viz naznačený důkaz (skoro) nulové singulární hodnoty
► skrytá degenerace systému
► může vést na jedno nebo i více přesných řešení
velmi malé singulární hodnoty
► ukazují na nízkou citlivost problému
► právě ve směrech odpovídajících sloupců V
► zpravidla lépe vynulovat v 2'
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
□     s   4 :
■OQ.O
21/30
Dekompozice na singulární hodnoty
Aproximace matíc
původní matici lze vyjádřit
A
k
je-li většina 07 skoro nulových
má smysl ukládat jen několik sloupců U a V
► stále dostáváme poměrně přesnou aproximaci A násobení Ax je výrazně efektivnější - K(M + N) operací
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Dekompozice na singulární hodnoty
Algoritmus
první fáze - redukce na bidiagonální formu
► využívá Householderovy transformace
druhá fáze - iterační, varianta výpočtu vlastních hodnot výjimečně stabilní
► zaměření na vytažení problematických vlastností A do 2 použijeme existující implementaci:-)
► už to za nás jednou někdo udělal
► další vylepšující triky dodavatelů knihoven
► nezbavuje to odpovědnosti za interpretaci výsledku
původní algoritmus Golub a Reinsch, Singular value decomposition and least squares solutions, 1970
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Vlastní hodnoty a vektory
PA081: Programování numerických
výpočtů
A. Křenek
Přehled a
► vlastní hodnoty a vektory matice (transformace) motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Dekompozice na singulární hodnoty
Vlastní hodnoty a vektory
► obecná reálná matice nemusí mít reálné vlastní hodnoty
► více viz kursy lineární algebry
Vlastní hodnoty a vektory
Použití
hledání kořenů polynomů ► výpočty vyšších mocnín matíc
An = VAV^VAV"1... VAV"1
VAnV"1
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Dekompozice na singulární hodnoty
Mastili hodnoty a vektory
4 □ ►   4       ►   4 = >   < -Š ►
Vlastní hodnoty a vektory
Použití
► hledání kořenů polynomů
► výpočty vyšších mocnín matíc
A" = VAV^VAV"1... VAV-1 = VA"^1
kritické body funkce více proměnných
► první parciální derivace jsou nulové
► Hessián (matice druhých parciálních derivací) určuje aproximaci kvadrikou v daném bodě
► symetrická matice - reálné hodnoty, ortonormální vektory
► extrémy - všechny A kladné, sedla - některé záporné
► absolutní hodnoty A určují tvar, vektory orientaci
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
□     s   4 :
■OQ.O
25/30
Vlastní hodnoty a vektory
Použití
► kritické body vektorového pole
► velikost vektoru je nulová
► Jakobián (matice prvních parciálních derivací)
Rapelling focua R1 .R2> 11 - -I2 <>
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
26/30
Vlastní hodnoty a vektory
Použití
► analýza hlavních komponent
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
téma příští přednášky
□        4 S ►   4 =  >   < -š ►
■O0.O
Vlastní hodnoty a vektory
Algoritmy
iterační algoritmus
► A0 := A
► dekompozice A^ = Q^R^
► položíme Afc+i := R^Q^
platí Afc+i = RkQt = Qj^RfcQfc = Q^AfcQfc
tedy iterační krok zachovává vlastní hodnoty
lze ukázat, že za jistých okolností konverguje k Shurově formě
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
Vlastní hodnoty a vektory
Algoritmy
numericky dost nepříjemný problém
ještě jasnější případ, kdy sáhnout k hotovým řešením
různé varianty pro různé případy
► reálné a komplexní
► vlastní hodnoty, vektory, obojí
► různé speciální typy matic
vztah k singulárním hodnotám
► sloupce U v SVD jsou vlastní vektory AAT
► nenulové singulární hodnoty jsou odmocniny nenulových vlastních hodnot AAT
Přehled a motivace
LU
dekompozice
Choleského dekompozice
QR
dekompozice
Vlastní hodnoty a vektory
■O0.O
29/30
Shrnutí
► základní rozklady matic
► LU, Cholského, QR
► použití především k řešení systémů lineárních rovnic
► existují i další varianty
► SVD
► náročnější postup, větší stabilita
► špatně podmíněné systémy
► větší náhled do problému - další aplikace
► vlastní hodnoty
► rozklad matice A = VAV-1
► přímé využití pro charakteristiku řady problémů
► (více v příští přednášce)
► existující implementace
► téměř vždy je použijeme
► je nutné vědět, co děláme a co můžeme od dané implementace čekat