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 2019 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 ► model vyšetření ultrazvukem ► zdroj ultrazvuku, šíření a odraz vln v přístroji a tkáních ► 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 centimetrů ► 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-30 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 20 jader = 2, 56TB/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-512 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; i++) a[i] += b[i]; movaps 32(%rsp), %xmmO addps (%rsp), %xmmO movaps %xmmO, 32(%rsp) ► ruční řešení (intrinsic funkce) je neportabilní, rychle zastarává ► vektorizaci kódu provede chytřejší kompilátor ► 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 singulární/^ Vektorové instrukce ► zarovnání dat #define SIZ 200 #define DIM 3 f 1 oat x [SIZ] [DIM] , y [SIZ] [DIM] , z [SIZ] [DIM] ; foř (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 při naivním přístupu ► N ~ 10 - zpravidla stačí float ► N -50 - double ► větší - vyžadují chytřejší metody (pivoting) 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 ► přímé ► 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 ... 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 siriffu a 23/77 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 si neum 24/77 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 ► numpy.linalg a scipy.linarg https: //sci py. org/ ► Python API, interně využívá zpravidla BLAS/LAPACK ► 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 = 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 sineu a 26/77 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í | akk ► parciální i plný pivoting jsou pak numericky stabilní ► přínos plného pivotingu je jen okrajový < 1 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 siriffu a 29/77 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.^.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 si neum 33/77 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/77 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 si neum 36/77 LU dekompozice Shrnutí a použití PA081: Programování numerických výpočtů A. Křenek ► ř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 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 rW77 LU dekompozice Shrnutí a použití ► tzv. driver routine v knihovně LAPACK ► vlastní rozklad LU ► řešení systému ► iterační zpřesnění (viz dále) LDA, ldb = LDB, info; int n = N, nrhs = NRHS, Ida = i nt i piv[N]; float a[LDA*N] = { ... }; float b[LDB*NRHS] = { ... }; sgesv( &n, &nrhs, a, &lda, i piv, b, &ldb, &info ); ► kompilace icc -mkl source.c ► existují i volné verze 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árnV^ 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/77 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 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/^ 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 si neum řW77 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/77 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/77 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 siriffu a řW77 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 singtt a rW77 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 sineu a rW77 Dekompozice na singulární hodnoty Použiti pro M = N PA081: Programování numerických výpočtů A. Křenek ► řešení 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í ^ crmin/crmax < € (š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 řW77 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í ► crmin/crmax < € (špatně podmíněná matice) - standardní metody řešení selhaly ► označme A -1 V[diag(l/Oi)]U T Oi = 0 ► rovnice Ax = b nemusí mít řešení, přesto zkusíme x = vrurb 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 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 siriffu a rW77 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 sineu a rW77 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 k ► je-li většina di 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í 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 rW77 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 S ► 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 Vlastni hodnoty a vektory Programování numerických výpočtů A. Křenek ► vlastní hodnoty a vektory matice (transformace) Ax = Ax WM" IP' 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 rW77 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 sineua rW77 Vlastní hodnoty a vektory Algoritmy PA081: Programování numerických výpočtů A. Křenek ► iterační algoritmus ► A0 := A ► dekompozice A^ = Q^R^ ► položíme Afc+i := RfcQk ► platí Ak+1 = RfcQfc = QlChRkQk = QkľAk(h ► tedy iterační krok zachovává vlastní hodnoty ► lze ukázat, že za jistých okolností konverguje k Shurově formě 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 smmna 59/77 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 sineu a 63/77 Ří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/77 Ří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 singtt a 65/77 Ří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 i dx [] 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 [] 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 rW77 Ří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 rW77 Ří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 sineu a 68/77 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 krojů je 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 fW77 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 siriffu a 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 sineu a 72/77 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 si neum 73/77 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/77 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/77 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í 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 /nV77