A. Křenek 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 2016 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 singulárnV^ 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 ... ► ► 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 singulárrfí/^ 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 ► ► 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 singulárrfí/^ 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 \ dL/dt J ( —P \ m F V pro simulaci potřebujeme opakovaně řešit y n = Yn-1 + dy y n ► triková aproximace - zanedbání vyšších členů Taylorova rozvoje dy/dt jako funkce y ► konečný krok od yn-\ k yn znamená řešení systému 13 lineárních rovnic 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 singulární/^ Motivace Modely měkkých tkání PA081: Programování numerických výpočtů A. Křenek ► ► ► ► ► 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) 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 singulární/^ 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) 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 ► taktování typicky 2-4 GHz ► 4-16 jader v socketu ► v jednom taktu na jednom jádru až 8 aritmetických operací ve float, tj. až 32 Gflop/s ► 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 = 1, 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 ; f oř (i =0; i <4; i ++) a [i ] += b [i ] ; movaps 32(%rsp), %xmm0 addps (%rsp), %xmm0 movaps %xmm0, 32(%rsp) ► ruční řešení (intrinsic funkce) je neportabilní, rychle zastarává ► 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 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 singnlárrft/^ Vektorové instrukce ► zarovnání dat #define SIZ 200 #define DIM 3 float x[SIZ][DIM],y[SIZ][DIM],z[SIZ][DIM]; for (i=0; i N - příliš podmíněné systémy ► přesné řešení zpravidla neexistuje ► hledáme „nejlepší kompromis" 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 a 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é ► O (N2) vs. O (N) nenulových prvků 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 a 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é ► O (N2) vs. O (N) nenulových prvků ► velikost - dopady kumulace chyby ► N ~ 50 - zpravidla stačí float ► N -200-500 - double ► větší - vyžadují speciální metody 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 a 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 ► ► ► ► ► ruzna náročnost pro ruzne vstupy ► Jacobi, Gauss-Seidel, sdružené směry ... 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 a 23/76 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 a 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.netlib.org/lapack/ ► 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ší... Gaussova eliminace Základní postup PA081: Programování numerických výpočtů A. Křenek ► 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, ■ ■ ■, &m\ odečtením ^ násobku prvního řádku ► obdobně pokračujeme druhým sloupcem od ► výsledkem je matice s vynulovanými prvky pod diagonálou ► řešení systému získáme zpětnou substitucí ► poslední rovnice UmmXm = b m je triviální ► do předposlední dosadíme získanou hodnotu Xm atd. 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 singtt a 26/76 Gaussova eliminace Numerické problémy ► diagonální prvek a\± je 0 nebo příliš malý ► v fe-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 zaokrouhlení může vést až k 0122 = -7, to odpovídá původní matici T 6 1 1 0 ► metoda v této podobě je numericky nestabilní e 1 0 1-i _ Gaussova eliminace Pivoting ► 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 fc-tého sloupce ► částečný pivoting (parciální) ► v iteraci fc vybíráme z Ujk pro j > fc, tj. jen z nulovaného sloupce ► plný pivoting ► vybíráme z aij pro í,j > fc, tj. z celé podmatice doprava a dolů ► vyžaduje udržovat vznikající permutaci proměnných Gaussova eliminace Pivoting PA081: Programování numerických výpočtů A. Křenek za pivota volíme maximum z kandidátů ► pro všechny faktory v eliminačních krocích platí aík akk < 1 ► parciální i plný pivoting jsou pak numericky stabilní ► přínos plného pivotingu je jen okrajový 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 a 29/76 Gaussova eliminace Pivoting ► 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ý ► 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 Gaussova eliminace Gauss-Jordanova 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 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 A1 a následné x = A_1b Rozklady matic ► danou matici A vyjádříme jako součin .A. — .A.1.A.2 ... -A.^ ► 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í 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ě deftnitní A ► A = QR ► Q ortogonální, R horní trojúhelníková ► symetrické varianty RQ, QL a LQ ► singulární hodnoty: A = UZV ► U, V ortogonální, S diagonálni ► 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 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 a 33/76 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é) ► Uje 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í LU dekompozice Výpočet rozkladu prvky rozkladu A = LU lze, s ohledem na nuly psát i < j- CLij í = j: Uíj í > j: aij Uiihj + utzhj + ■ ■ ■ + uuUj unhj + utzhj + ■ ■ ■ + uuljj unhj + Ut2hj + ■ ■ ■ + Uijljj ► 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 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 singula 35/76 LU dekompozice Croutův algoritmus PA081: Programování numerických výpočtů A. Křenek ► postupně pro j = 1, 2,..., N: ► pro í = 1,2,... j spočítáme na základě uvedených rovnic í-1 Uij — Cli j ^ hkUkj k=l pro í = j + 1, j + 2,...,N ► postup vždy využívá dříve spočtené prvky ► k numerické stabilitě je třeba navíc dodat pivoting Ijj 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 singtt a 36/76 LU dekompozice Shrnutí a použití ► řešení lineárních rovnic pro libovolný počet b ► Croutův algoritmus O (AT3) ► dopředná a zpětná substituce O (AT2) ► 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í A1 LU dekompozice Shrnutí a použití PA081: Programování numerických výpočtů A. Křenek int n = N, nrhs = NRHS, Ida = i nt i piv[N]; float a[LDA*N] = { ... }; float b[LDB*NRHS] = { ... }; LDA, ldb = LDB, info; sgesv( &n, &nrhs, a, &lda, i piv, b, &ldb, &info ); 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 singula 38/76 Choleského dekompozice ► předpokládáme A = LLr ► po prvcích platí lu = i-1 au -1 ]2 lji = k=0 i ( Tu i-1 k=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. 2x méně operací než LU, menší paměťová náročnost má smysl používat speciální funkce z knihoven ► ► 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 a 39/76 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 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 Dekompozicí na singulártr)/^ QR dekompozice Gram-Schmidtův proces ► Gram-Schmidtův ortogonalizační proces ui = ai u2 = a2 - projVla2 u3 = a3 - projVla3 - projV2a3 vi = v2 = v3 = Ul Ui u2 u2 u3 u3 potom / vi ■ ai vi ■ a2 ... Q= (vi ...vn) R = 0 v2 ■ a2 .. numerický problém, jsou-li některé a^, aj skoro kolmé 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 a QR dekompozice Householderova transformace PA081: Programování numerických výpočtů A. Křenek ► zrcadlení daného vektoru podle nadroviny ► zarovná s bázovým vektorem, zachová velikost ► pro libovolné x: u = x + ŕxei u v = u Q = I - 2vv t ► potom Qx = (a, 0,..., 0) t 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 a 42/76 QR dekompozice Householderova transformace PA081: Programování numerických výpočtů A. Křenek ► zrcadlení daného vektoru podle nadroviny ► zarovná s bázovým vektorem, zachová velikost ► pro libovolné x: u = x + ŕxei u v = u ► potom Qx = (a, 0,..., 0)r ► triangulace A R = Qn-i ■ ■ ■ QiA a tedy Q = / - 2vv t 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 a 42/76 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ě deftnitní 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 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 a 43/76 Dekompozice na singulární hodnoty Základní tvrzení ► libovolnou reálnou (i komplexní) matici lze rozložit AmxN = UmxN " ^NxN " ^NxN ( = ^MxM ' ^MxN ' ^ NxN^ ► Uje sloupcově ortogonální ► až na nulové sloupce v případě M < N ► Z 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 Dekompozice na singulární hodnoty Geometrický význam PA081: Programování numerických výpočtů A. Křenek ► A je složení transformací ► rotace/zrcadlení V1 ► zvětšení/zmenšení faktory 07 ve směrech včetně degenerace (07 = 0) ► rotace/zrcadlení a projekce do méně/více dimenzí U 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 sin gum 45/76 Dekompozice na singulární hodnoty Geometrický význam ► A je složení transformací ► rotace/zrcadlení V1 ► zvětšení/zmenšení faktory 07 ve směrech 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 Dekompozice na singulární hodnoty Geometrický význam ► A je složení transformací ► rotace/zrcadlení V1 ► zvětšení/zmenšení faktory 07 ve směrech 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 JSf (A) = {x e Rn : Ax = 0} ► řádky Vr odpovídající nulovým 07 jsou jeho generátory Dekompozice na singulární hodnoty Numerický význam PA081: Programování numerických výpočtů A. Křenek ► „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í 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 si neum 46/76 Dekompozice na singulární hodnoty Použiti pro M = N PA081: Programování numerických výpočtů A. Křenek reseni systému rovnic, resp. výpočet inverzní matice A"1 = VtdiagU/cr^U7 ► kdy to nejde ► jedno nebo více 07 je nulových - A byla singulární * ^min/cTmax < € (špatně podmíněná matice) - standardní metody řešení selhaly 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 a Dekompozice na singulární hodnoty Použiti pro M = N PA081: Programování numerických výpočtů A. Křenek reseni systému rovnic, resp. výpočet inverzní matice A"1 = V[diag(l/o-i)]Ur ► kdy to nejde ► jedno nebo více 07 je nulových - A byla singulární * ^min/cTmax < € (špatně podmíněná matice) - standardní metody řešení selhaly ► označme diag 1/CTí (Ji ^ 0 0 Oi = 0 ► rovnice Ax = b nemusí mít řešení, přesto zkusíme x = vrurb 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 a Dekompozice na singulární hodnoty Použiti 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' | Ax - b + b' I = I (UZVr) (VZ'Urb) - b + b' | = |(UZZ'Ur-J)b + b'| = |U((ZZ' - J)Urb + Urb')| = |(ZZ' - J)Urb + Urb'| ► ES' - I je diagonální s nenulovými prvky pro ai = 0 ► b' je v oboru hodnot A, tedy Urb' má nenulové prvky právě pro 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 Dekompozice na singulární hodnoty Použiti pro M * N PA081: Programování numerických výpočtů A. Křenek ► 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í 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 si neum Dekompozice na singulární hodnoty Použití pro m =é N ► 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í ► 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 X' Dekompozice na singulární hodnoty Aproximace matic PA081: Programování numerických výpočtů A. Křenek původní matici lze vyjádřit Aij = ŽVkUikVjk k ► je-li většina 0 11 =12 = 0 Repelling nodo R1 ,R2 > 0 n = 12 = 0 Attracting focus R1 = R2 < 0 11 =-I2oO Repelling focus R1 =R2>0 H = -I2 <> 0 Center R1 = R2 = 0 11 = -11 o 0 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 si neum Vlastní hodnoty a vektory Použití analýza hlavních komponent 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 Gaus sova eliminace Rozklady matic téma příští přednášky Dekompozice na sineu a Vlastní hodnoty a vektory Algoritmy ► iterační algoritmus ► A0 := A ► dekompozice A^ = Q^R^ ► položíme Afc+i := R^Q^ ► platí Ak+1 = Rk(h = Q[(hRkQk = Q^A^ ► tedy iterační krok zachovává vlastní hodnoty ► lze ukázat, že za jistých okolností konverguje k Shurově formě 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 sinmila 59/76 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 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 = VAV1 ► 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 Aôx = 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 ôx korigujeme původní x Iterační zpřesnění ► schematické znázornění postup lze opakovat, zpravidla stačí 1-2 krát 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 si neum 63/76 Řídké matice ► ► ► typicky rozsáhlé problémy (miliony rovnic) ► ale počet nenulových koeficientů je malý, zpravidla O (AT) ► 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 (AT) 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 ► ► ► ► ► ► 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 a 64/76 Řídké matice Některé speciální vzory 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 a 65/76 Řídké matice Uložení v paměti PA081: Programování numerických výpočtů A. Křenek ► 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 IN + 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 [] 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 a Řídké matice Uložení v paměti PA081: Programování numerických výpočtů A. Křenek původní matice 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 je reprezentována 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 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 a Ří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 0 v) ► lze ukázat (A ) = (A + (u 0 v)) -A----—r- 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 (AT), náročnost celé metody je O (AT) ► postup lze opakovat pro různá u, v ► existují další zobecnění (Woodbury, viz literatura) -1 PA08i: 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 si neum 68/76 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 krojuje zpravidla O (AT) nevyžaduje další paměť téměř singulární systémy ► přímé metody ztrácí přesnost kumulací chyby ► ► ► 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 si neum £9/76 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ě ► lim^^oo I (A')k| = 0 pro nějakou maticovou normu ► Frobeniova norma kde p{) je maximální absolutní hotnota vlastních hodnot matice ► maximální součet sloupce nebo řádku ► spektrální norma A P (A'A) Iterační metody Prostá iterace - příklady konkrétních metod ► Jacobi: z rozkladu A = D - L - U dostaneme x = D _1(L + U)x + D_1b, tj. x L ( Mi \ 11 bi- (k) X aíjxj j ► Gauss-Seidel: x = (D - D^Ux + (D - L)"1^ L ( Mi \ i-l n \ J = l 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 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 si neum 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 1 f(x) = -xAx-bx ► v minimu je její gradient nulový, tj. 0 = V f = 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, lze rozšířit ► při výpočtu je třeba jen násobit Ax ► to lze u řídkých matic velmi efektivně 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 singn a 72/76 Bloková LU dekompozice rekurzivní formulace algoritmu Aoo aoi A02 d10 «11 d12 A20 a2i A22 ) ► skalární a bloková verze Au a2i = t a12 A22 = a2i/rxn zůstává A22 - a2iai2 11 A21 A12 A22 A21 LllAl2 A22 - L2lAl2 ► podstatná část operací je násobení matic ► takto implementuje knihovna LAPACK PA08i: 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 a 73/76 Bloková LU dekompozice rekurzivní formulace algoritmu Aoo aoi A02 d10 «11 d12 A20 a2i A22 ) ► skalární a bloková verze Au a2i = t a12 A22 = a2i/rxn zůstává A22 - a2iai2 11 A21 A12 A22 A21 LllAl2 A22 - L2lAl2 ► 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 2p x p PA08i: 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 a 73/76 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 Super Ma trix ► 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 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 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 a 76/76 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í