PA081: Programování numerických výpočtů 5. Systémy lineárních rovnic a optimalizované implementace lineární algebry Aleš Křenek jaro 2015 Motivace aparát vektorových a maticových operací součást metod řešení nelineárních rovnic, optimalizací atd. diferenciální rovnice ► simulace dynamických systémů ► metoda konečných prvků realistické osvětlení scény (radiosita) a mnoho dalších ... Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Motivace aparát vektorových a maticových operací součást metod řešení nelineárních rovnic, optimalizací atd. diferenciální rovnice ► simulace dynamických systémů ► metoda konečných prvků realistické osvětlení scény (radiosita) a mnoho dalších ... relativně snadná implementace, podpora v HW Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Motivace Simulace pohybu tělesa ► založena na Newtonově zákoně F = ma ► vede na systém diferenciálních rovnic (13 proměnných) dy(t) dt ( dx/dt \ dq/dt dP/dt y dLjdt i í -P \ m \wqq F T pro simulaci potřebujeme opakovaně řešit Yn ~ Yn—1 dt triková aproximace - zanedbání vyšších členů Taylorova rozvoje dy/dt jako funkce y konečný krok od yn_i k yn znamená řešení systému 13 lineárních rovnic Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Motivace Modely měkkých tkání Interakce s měkkými tkáněmi ► tréning chirurgů - vytvořit simulaci chování tkání ► potřebujeme vytvořit matematický model tkáně ► s touto tkání pak můžeme interagovat (např. hapticky) ► je třeba simulovat chování tkáně s dostatečnou přesností použité matematické modely ► deformovatelná tělesa ► chování je popsáno parciálními diferenciálními rovnicemi, zpravidla neznáme analytické řešení ► řešíme numericky pomocí diskretizace metodou konečných prvků (FEM) Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sinmjlár 4/76 Motivace Metoda konečných prvků ► simulovaná tkáň je geometricky aproximována meshem deformace popisuje teorie elasticity geometricky lineární model - systém lineárních rovnic geometricky nelineární model - systém nelineárních rovnic ► v každém kroku iterační metody systém lineárních rovnic systémy rovnic jsou řídké (ovlivňují se jen přímo spojené uzly) Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Motivace Rekonstrukce signálu z ultrazvuku pro vývoj filtrů rekonstruujících data z ultrazvuku je výhodné mít počítačový model přístroje a vyšetřovaného objektu snáze pak ladíme metody rekonstrukce, nevnášíme nepřesnost danou nedokonalým modelem vyšetřovaného objektu či nedokonalým přístrojem pro výpočet šíření vln v objektu použijeme FEM ► na vlnovou délku ultrazvuku potřebujeme několik lineárních elementů ► velikost elementu je pak v desetinách mm ► vyšetřovaný objekt má velikost v jednotkách až desítkách centimetrech systémy o jednotkách až desítkách miliónů rovnic ► vysoké nároky na výpočetní kapacitu ► vysoké paměťové nároky Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Výkon současných CPU ► taktování typicky 2-4 GHz ► 4-16 jader v socketu ► v jednom taktu na jednom jádru až 8 aritmetických operací ve f 1 oat ► podmíněno dostatečným přísunem instrukcí a dat ► potřebný datový tok 4 GHz x 8 operací x 4 B x 10 jader = l,28TB/s Vektorové instrukce ► historie - např. Cray v předsálí budovy D ► současné vektorové systémy (např. NEC SX) ► řešení velmi specializovaných problémů ► nízký výkon na skalárních operacích ► vektorová rozšíření v architektuře x86 ► MMX: pouze celá čísla, kolize s float registry ► SSE: 8 (16) registrů po 128 bitech, postupně přidávané instrukce (rsqrt) ► AVX: 256 bitové registry, 3-operandové instrukce ► snazší využití plného výkonu procesoru Vektorové instrukce triviální příklad float a[4],b[4]; i nt i; for (i=0; i<4; a[i] += b[i]; movaps addps movaps 32(%rsp), %xmmO (%rsp), %xmmO %xmmO, 32(%rsp) vektorizaci kódu provede chytřejší kompilátor (icc) ► musíme hlídat, abychom tomu nebránili ► zarovnání dat, recyklace proměnných, vedlejší efekty in-line funkcí, ... ► vyplatí se kontrola vygenerovaného assembleru Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sinmjlár 9/76 Vektorové instrukce zarovnání dat #define SIZ 200 #define DIM 3 f 1 oat x [SIZ] [DIM] , y [SIZ] [DIM], z [SIZ] [DIM] ; f or (i=0; i Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Asymptotická složitost ► časová složitost násobení matic je 0(n3) ► Strassen (1969): 0(n281) ► rozdělení matic na 2 x 2 bloky, potom Mi = (Au +A22)(Bii +B22) M2 = (A21 + A22)Bii M3 = An(Bi2 +B22) M4 = A22(B2i - B11) M5 = (Au + Ai2)B22 M6 = (A2i -Aii)(Bn +B12) M7 = (A12 -A22MB21 +B22) C11 = Mi + M4 - M5 + M7 C12 = M3 + M5 C2i = M2 + M4 C22 = Mi - M2 + M3 + M6 rekurzivní aplikace, počet operací T(n) T{n) = 7T(n/2) + 18(n/2)2 = 0(nlo&7) Coppersmith-Winograd (1987): 0(n2376) 0(n 2.81 > Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Systémy lineárních rovnic hledáme řešení systému maticové vyjádření Ax = b Klasifikace problémů *■ M = N - dobře podmíněné systémy ► naděje na jednoznačné řešení (je-li A regulární) ► M < N - nedostatečně podmíněné systémy ► žádné řešení ► nekonečně mnoho řešení (N - M-rozměrný prostor) ► M > N - příliš podmíněné systémy ► přesné řešení zpravidla neexistuje ► hledáme „nejlepší kompromis" Klasifikace problémů *■ M = N - dobře podmíněné systémy ► naděje na jednoznačné řešení (je-li A regulární) ► M < N - nedostatečně podmíněné systémy ► žádné řešení ► nekonečně mnoho řešení (N - M-rozměrný prostor) ► M > N - příliš podmíněné systémy ► přesné řešení zpravidla neexistuje ► hledáme „nejlepší kompromis" husté vs. řídké ► 0(N2) vs. O(N) nenulových prvků Klasifikace problémů M = N - dobře podmíněné systémy ► naděje na jednoznačné řešení (je-li A regulární) M < N - nedostatečně podmíněné systémy ► žádné řešení ► nekonečně mnoho řešení (N - M-rozměrný prostor) M > N - příliš podmíněné systémy ► přesné řešení zpravidla neexistuje ► hledáme „nejlepší kompromis" husté vs. řídké ► 0(N2) vs. O(N) nenulových prvků velikost - dopady kumulace chyby ► N ~ 50 - zpravidla stačí f 1 oat ► N -200-500 - double ► větší - vyžadují speciální metody Přehled metod prime ► daný algoritmus, vždy stejný počet operací ► nemusí fungovat dobře pro „skoro singulární" nebo příliš velké systémy ► preferované pro „normální" problémy ► Gaussova eliminace, různé dekompozice, ... iterační ► postupně vylepšované řešení dosažení kritéria konvergence ► různá náročnost pro různé vstupy ► Jacobi, Gauss-Seidel, sdružené směry ... Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Přehled metod ► specializované metody pro řídké matice ► pro různé vzory řídkých matic ► výrazně efektivnější, jediná možnost řešení velkých systémů ► přímé i iterační ► metody pro singulární systémy ► analyticky i numericky („skoro") singulární ► nalezení celého prostoru řešení ► standardní metody zhavarují ► dekompozice na singulární hodnoty řešení příliš podmíněných systémů ► specializované metody nejmenších čtverců PA081: Programování numerických výpočtů A. Křenek Motivace Výkon současných CPU Blokové algoritmy BLAS Asymptotická složitost Systémy lineárních rovnic Přehled metod Gaussova eliminace Rozklady matic LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineu 24/76 Dostupné knihovny zpravidla sáhneme k hotové knihovně ► nejsme první, kdo tento problém řeší ► implementace optimalizované pro různé případy ► k volbě vhodné metody je třeba rámcově rozumět jejich principům Numerical Recipes in C ► hlavní zdroj pro tuto přednášku ► jednoduché přehledné implementace LAPACK http://www.netli b.org/1apack/ ► plnohodnotná volně dostupná knihovna ► využívá nízkoúrovňové knihovny BLAS, zpravidla implementované strojově závisle ► obecné a specializované funkce pro různé typické případy NAG http://www.nag.co.uk/ ► rozsáhlá komerční numerická knihovna a další... Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 25/76 Gaussova eliminace Základní postup standardní „školní" metoda ► řešení systému se nezmění, nahradíme-li libovolný řádek A a odpovídající prvek b lineární kombinací tohoto řádku s libovolným dalším vynulujeme 0,21, cím\ odečtením ^ násobku prvního řádku obdobně pokračujeme druhým sloupcem od a^2 výsledkem je matice s vynulovanými prvky pod diagonálou řešení systému získáme zpětnou substitucí ► poslední rovnice aMMxM = bM je triviální ► do předposlední dosadíme získanou hodnotu Xm atd. Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Gaussova eliminace Numerické problémy ► diagonální prvek akk je 0 nebo příliš malý ► v k-té iteraci se jím dělí ► dojde k přetečení ► při odečítání dochází k přílišné ztrátě platných cifer " e 1 " ' e 1 1 1 -N 0 1 - zaokrouhlení může vést až k 0,22 = -7, to odpovídá původní matici I" e 1 " 1 0 ► metoda v této podobě je numericky nestabilní Gaussova eliminace Pivotíng ► další tvrzení o systému ► řešení systému se nezmění výměnou libovolné dvojice řádků ► řešení systému se nezmění výměnou libovolné dvojice sloupců, zaměníme-li zároveň příslušné komponenty vektoru x ► pivot je prvek, který po aplikaci těchto kroků použijeme k dělení při eliminaci fe-tého sloupce ► částečný pivoting (parciální) ► v iteraci k vybíráme z ojt pro j > k, tj. jen z nulovaného sloupce ► plný pivoting ► vybíráme z ay pro i, j > k, tj. z celé podmatice doprava a dolů ► vyžaduje udržovat vznikající permutaci proměnných Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Gaussova eliminace Pivotíng ► za pivota volíme maximum z kandidátů ► pro všechny faktory v eliminačních krocích platí | ^ | < 1 parciální i plný pivoting jsou pak numericky stabilní přínos plného pivotingu je jen okrajový Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Gaussova eliminace Pivotíng ► za pivota volíme maximum z kandidátů ► pro všechny faktory v ehminačních krocích platí | ^ | < 1 parciální i plný pivoting jsou pak numericky stabilní přínos plného pivotingu je jen okrajový volba pivota závisí na řádu koeficientů původního systému ► vynásobení jedné rovnice faktorem 1010 ... ► za pivota lze volit koeficient, který by byl největší, kdyby byl celý původní systém normalizovaný, tj. největší koeficient všech rovnic roven 1 ► vyžaduje dodatečnou údržbu faktorů škálování, může přinést větší robustnost Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 29/76 Gaussova eliminace Gauss-Jordánova eliminace Gaussovu metodu lze použít pro více pravých stran současně speciálně pro N bázových vektorů dostaneme výpočet matice A-1 rozšíření - Gauss-Jordanova eliminace ► v každé iteraci eliminujeme prvky nad i pod diagonálou A ► výsledkem je jednotková matice ► řešení systému je přímo modifikovaný vektor b ► resp. dostáváme A-1 srovnatelné se základní implementací ► reálně provádíme tytéž operace Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na Sinei] 30/76 Gaussova eliminace Výhody a nevýhody ► jednoduchá, dobře pochopitelná a stabilní metoda ►•je třeba znát pravé strany dopředu ► jsou součástí celého výpočtu ► výpočet A-1 je zatížen poměrně velkou numerickou chybou ► přímý výpočet řešení systému Ax = b je přesnější než výpočet A-1 a následné x = A_1b Rozklady matic danou matici A vyjádříme jako součin A = AiA2 ...An matice Ai, A2,... An mají nějaké speciální vlastnosti 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 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í Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Rozklady matic ► 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ě deflnitní 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álni, 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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í Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice LU dekompozice Výpočet rozkladu ► prvky rozkladu A = LU lze, s ohledem na nuly psát i < j: ciij = unhj + ui2hj + ■ ■ ■ + uuUj i = j: a.ij = unhj + ui2hj + ■ ■ ■ + uuljj i> j: ciij = unhj + ui2hj + ■ ■ ■ + »,,/,, ► je to systém N2 rovnic v N2 + N neznámých ► diagonála je pokryta uil ► můžeme tedy volit la = 1 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 + 1, 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 QR dekompozice 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 0(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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 37/76 Choleského dekompozice ► předpokládáme A = LLr ► po prvcích platí la a i i-1 lile fc=0 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. 2 x méně operací než LU, menší paměťová náročnost ► má smysl používat speciální funkce z knihoven Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 QR dekompozice Gram-Schmidtův proces Gram-Schmidtův ortogonalizační proces ui = ai u2 = a2 - projVla2 u3 = a3 - projVla3 - projV2a3 Vl v2 v3 lluil U2 u2 U3 l|u3| Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic potom Q = (vi ... V* R / vi ■ ai vi ■ a2 0 V2 ■ 3iZ \ numerický problém, jsou-li některé a;, a^ skoro kolmé Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 U T — Q=/-2wr u Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 = —^— Q = J - 2w l|u|| potom Qx = (a, 0,..., 0)r triangulace A R = Qn !... QiA a tedy Q = Q[qJ ... 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ě definitní 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 nulový prostor JV(A) = {x G RN : Ax = 0} ► řádky Vr odpovídající nulovým 07 jsou jeho generátory Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 Z ► singularity A odpovídají nulovým 07 ► včetně numerických (07 ^ 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í Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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-;)]Ur kdy to nejde ► jedno nebo více 07 je nulových - A byla singulární * ^min/^max < £ (špatně podmíněná matice) - standardní metody řešení selhaly Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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' Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na singulární hodnoty Použití pro M = N null space of A, solutions of solutions of / sh. ■ X = C' A-x = d^W~ \ \ SVD "solution" \ V-^^ of A x = c .. \ range oi A \ \ ^-^^ \ \ ^^^^ * \ \ \ / \ •c SVD solution of A x = d Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineii 48/76 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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í Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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í 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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í Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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' Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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í Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Vlastní hodnoty a vektory PA081: Programování numerických výpočtů A. Křenek Motivace ► vlastní hodnoty a vektory matice (transformace) I výkon současných CPU Blokové algoritmy B LAS Asymptotická složitost Systémy lineárních rovnic Přehled metod Gaussova eliminace ► obecná reálná matice nemusí mít reálné vlastní hodnoty Rozklady matic ► více viz kursy lineární algebry dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 54/76 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 > 0 11 - -I2 o 0 Výkon současných CPU Blokové algoritmy Systémy lineárních rovníc Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 56/76 Vlastní hodnoty a vektory Použití ► analýza hlavních komponent téma příští přednášky 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ě Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 AAr ► nenulové singulární hodnoty jsou odmocniny nenulových vlastních hodnot AAr Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice 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 Iterační zpřesnění ► i přímé metody řešení lineárních rovnic ztrácí přesnost ► v závislosti na velikosti systému a poměru koeficientů ► kumulace chyb pocházejících ze sčítání/odčítání ► v optimistickém případě 2-3 platná místa ► u „skoro singulárních" matic podstatně horší ► lze vylepšit iteračním procesem ► x je přesné řešení Ax = b, známe ale jen x + 5x; položíme A(x + 5x) = b + 5b odečtením dostaneme Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic A5x = 5b = A(x + 5x) - b pravou stranu umíme spočítat, řešíme nový systém rovnic pro 5x ► pro výpočet pravé strany je nutná vyšší přesnost odčítání ► s výhodou využijeme recyklaci LU dekompozice výsledným 5x korigujeme původní x Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Iterační zpřesnění ► schematické znázornění -»-f A postup lze opakovat, zpravidla stačí 1-2 krát Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 62/76 Řídké matice typicky rozsáhlé problémy (miliony rovnic) ► ale počet nenulových koeficientů je malý, zpravidla O(N) v extrémním případě neřešitelné standardními metodami ► příliš velká paměťová a časová náročnost ► de-facto nepřekonatelné problémy s přesností výpočtu specializované metody ► využívají nulových koeficientů ► vyžadují jen O (N) operací i paměti ► závislé na konkrétním vzoru nenulových prvků ► např. LU dekompozice tridiagonální matice - jen N prvkový vektor navíc a 2 cykly o N - 1 iteracích v některých případech aproximační ► nenulové prvky se během řešení „množí" nad míru ► je třeba některé zanedbávat Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Řídké matice Některé speciální vzory Řídké matice Uložení v paměti problematická není jen výpočetní náročnost, ale zejména nároky na paměť ► N = 106 by vyžadovalo ve floatech 4 TB existují různá schémata, např. ► pole hodnot f 1 oat val [] a indexů i nt idx[] ► prvních N prvků val [] jsou všechny diagonální prvky, zpravidla bývají nenulové ► prvních N prvků i dx [] jsou indexy ve val [], kde jsou uloženy první nenulové nediagonální prvky jednotlivých řádků matice ► i dx [N + 1] ukazuje za poslední prvek val [] ► další prvky val [] jsou nenulové nediagonální prvky matice uspořádané po řádcích a sloupcích ► další prvky i dx [] jsou indexy sloupců odpovídajících prvků val [] Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sinmj 65/76 Řídké matice Uložení v paměti ► původní matice je reprezentována 3 0 1 0 0 0 4 0 0 0 0 7 5 9 0 0 0 0 0 2 0 0 0 6 5 1 2 3 4 5 6 7 8 9 10 11 idx[] 7 8 8 10 11 12 3 2 4 5 4 val[] 3 4 5 0 5 n/a 1 7 9 2 6 Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Dekompozice na sineij 66/76 Řídké matice Sherman-Morrisonův postup řídké matice, které se „trochu liší" od známého vzoru ► navíc řádek, sloupec apod. ► obecně A' = A + (u ® v) lze ukázat (A')"1 = (A+ (U0V))-1 = A"1 (A-1!!) 0 (vA'1) 1 + vA xu některou z hotových metod vypočteme A-1 aplikací vzorce získáme (A')-1 ► je-li A-1 řídká v řádu O(N), náročnost celé metody je O postup lze opakovat pro různá u, v existují další zobecnění (Woodbury, viz literatura) Iterační metody ► obecně méně přesné než přímé metody ► řešení řídkých systémů ► neexistuje přímá metoda pro konkrétní vzor ► výpočet iteračního kroju je zpravidla O(N) ► nevyžaduje další paměť ► téměř singulární systémy ► přímé metody ztrácí přesnost kumulací chyby Iterační metody Prostá iterace rovnici Ax = b převedeme na tvar x = A'x + b' opakovaně počítáme iterační krok kritérium konvergence ► zobrazení x <- A'x + b' musí být kontrakce ► stejné jako u nelineárních rovnic naplnění předpokladů věty o pevném bodě ► lmifc^K. l(A')fc| = 0 pro nějakou maticovou normu ► Frobeniova norma _ |A| spektrální norma |A| =p{A1A) kde pO je maximální absolutní hotnota vlastních hodnot matice maximální součet sloupce nebo řádku Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Iterační metody Prostá iterace - příklady konkrétních metod Jacobi: z rozkladu A = D - L - U dostaneme x = D X(L + U)x + D xb, tj. Gauss-Seidel: x = (D - L)_1Ux + (D - L) ^Ub , / i—l n au \ j=i j=í+i lze recyklovat uložení x konvergence není zaručena automaticky ► pro konkrétní A některé metody konvergují, jiné ne ► specifická kritéria viz literatura Iterační metody Metoda konjugovaných gradientů ► lze využít i metody řešení nelineárních rovnic ► smysl to má jen ve velmi specifických případech ► výjimkou je metoda konjugovaných gradientů ► minimalizujeme funkci /(x) 1 xAx - bx v minimu je její gradient nulový, tj. 0 = V/ = Ax - b tedy minimum / přesně odpovídá řešení Ax = b metoda tedy najde minimum po N iteracích ► / je kvadratická forma, viz minulá přednáška takto jednoduše funguje jen pro pozitivně definitní A, rozšířit při výpočtu je třeba jen násobit Ax ► to lze u řídkých matic velmi efektivně lze Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Bloková LU dekompozice ► rekurzivní formulace algoritmu / A, skalární a bloková verze a2i »12 = a2i/ofn zůstává A22 = A22 - a2iai2 podstatná část operací je násobení matic takto implementuje knihovna LAPACK LU dekompozice Choleského dekompozice QR dekompozice Bloková LU dekompozice ► rekurzivní formulace algoritmu / A, skalární a bloková verze a2i »12 = a2i/ofn zůstává A22 = A22 - a2iai2 podstatná část operací je násobení matic takto implementuje knihovna LAPACK problém algoritmu - „dlouhé" matice A21 a A12 Quintana et. al (2008): algoritmus s bloky 2 p x p LU dekompozice Choleského dekompozice QR dekompozice Paralelní blokové algoritmy ► paralelismus na úrovni BLAS ► existují optimalizované implementace ► implicitní synchronizace na konci každého volání ► nemusí být dosažitelný optimální výkon ► blokové operace jsou efektivní provedené vcelku ► všechna data se dostanou naráz do cache ► na úrovni L1/L2 se vyplatí na jednom jádru CPU ► vzájemně nezávislé blokové operace na více jádrech Paralelní blokové algoritmy Quintana et. al (2008) - prototypová implementace SuperMatrix ► konkrétni algoritmus proběhne „abstraktně" ► blokové operace s deklarovanými závislostmi se pouze zařadí do fronty ► následně se vyhodnotí pořadí zpracování ► operace se provedou potenciálně paralelně čitelná formulace algoritmu bez explicitního paralelismu efektivní provedení ► matice n > 5000 ► 16 jader CPU ► LU dekompozice na více než 50% teoretického výkonu Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Využití GPU téměř vše, co platí o cache, platí také o GPU ► rychlá paměť omezené velikosti ► netriviální režie kopírování z/do hlavní paměti paralelní kód ► daleko masivnější (desítky až stovky) viz PV197: GPU Programming Výkon současných CPU Blokové algoritmy Systémy lineárních rovnic Gaussova eliminace LU dekompozice Choleského dekompozice QR dekompozice Shrnutí ► řešení lineárních rovnic je součást řady problémů v operacích LA se tráví velká část výpočetního času ► má smysl hledat optimalizované a numericky stabilní algoritmy ► neexistuje univerzální řešení ► k dispozici řada optimalizovaných implementací