PA081: Programování numerických výpočtů A. Křenek PA081: Programování numerických výpočtů 7. Systémy lineárních rovnic Aleš Křenek jaro 2011 1/31 Přehled metod Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody Motivace ► hledáme řešení systému a\\X\ + aizXz + &2\X\ + 0.22X2 + + 0.2NXN = b2 dM\X\ + CLM2X-2 + ■ ■ ■ + 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(iV) 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(iV) 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 ► 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, LU 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 Přehled metod Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 8/31 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 (v příští přednášce) ► ř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 Přehled metod Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 9/31 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 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 a.21, o-mi odečtením ^ násobku prvního řádku obdobně pokračujeme druhým sloupcem od a.32 výsledkem je matice s vynulovanými prvky pod diagonálou řešení systému získáme zpětnou substitucí ► poslední rovnice (ImmXm = k>m je triviální ► do předposlední dosadíme získanou hodnotu % atd. Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 11/31 Gaussova eliminace Numerické problémy ► diagonální prvek 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 1 1 0 1 zaokrouhlení může vést až k a.22 = -\, to odpovídá původní matici " e 1 1 0 ► metoda v této podobě je numericky nestabilní Gaussova eliminace LU dekompozice Ilerační zpřesnění Řídké matice Iterační metody 12/31 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 k vybíráme z ajk pro j > k, tj. jen z nulovaného sloupce ► plný pivoting ► vybíráme z pro i, j > k, tj. z celé podmatice doprava a dolů ► vyžaduje udržovat vznikající permutaci proměnných Gaussova eliminace Pivoting ► za pivota volíme maximum z kandidátů ► pro všechny faktory v eliminačních krocích platí | ^ | ► parciální i plný pivoting jsou pak numericky stabilní ► přínos plného pivotingu je jen okrajový 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-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 Gaussova eliminace LU dekompozice Ilerační zpřesnění Řídké matice Iterační metody 15/31 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 *b 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í ► viz Gaussova metoda LU dekompozice Výpočet rozkladu ► prvky rozkladu A = LU lze, s ohledem na nuly psát i < j: atj = unhj + ui2hj + ■ ■ ■ + Uiikj í = j: atj = unhj + Ui2hj + ■ ■ ■ + uuljj í> j: atj = unhj + Ui2hj + ■ ■ ■ + i/,,'// ► je to systém N2 rovnic v N2 + N neznámých ► diagonála je pokryta u i i ► můžeme tedy volit la = 1 LU dekompozice Croutův algoritmus ► postupně pro j = 1, 2,... ,N: ► pro í = 1,2,... j spočítáme na základě uvedených rovnic í-i Uij = Clij — ^ líkUkj k=l ► pro í = j + 1, j + 2,... ,N 3-1 - I au - Ui ► postup vždy využívá dříve spočtené prvky ► k numerické stabilitě je třeba navíc dodat pivoting Ijj 1 ( 3 lij = - Clij - X ^ífc^fcj Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 19/31 LU dekompozice Shrnutí a použití ► řešení lineárních rovnic pro libovolný počet b ► Croutův algoritmus O (JV3) ► dopředná a zpětná substituce O (JV2) ► 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-l 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 + ôx; položíme A(x + ôx) = b + 5b odečtením dostaneme Aôx = 5b = A(x + ôx) - b pravou stranu umíme spočítat, řešíme nový systém rovnic pro ôx ► 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 Přehled metod Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 22/31 Řídké matice ► typicky rozsáhlé problémy (miliony rovnic) ► ale počet nenulových koeficientů je malý, zpravidla O(iV) ► 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(iV) 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 Přehled metod Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 23/31 Řídké matice Některé speciální vzory Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 24/31 Ří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 í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 [] Gaussova eliminace LU dekompozice Ilerační zpřesnění Řídké matice Iterační metody 25/31 Ří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 Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 26/31 Ří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 8 v) ► lze ukázat (AT^ÍA+Oiavjr^A-1- (A-^aívA-i) 1 + vA :u 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(N) ► postup lze opakovat pro různá u, v ► existují další zobecnění (Woodbury, viz literatura) Gaussova eliminace LU dekompozice Iterační zpřesnění Řídké matice Iterační metody 27/31 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(iV) ► 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ě ► limk^oo |(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 |A| ► spektrální norma |A| p(ArA) 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. \ j=í,j*i i ► Gauss-Seidel: x = (D - D^Ux + (D - D^Ub 11 V j=i j=í+i / l 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 směrů ► 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) = |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, 1 rozšířit ► při výpočtu je třeba jen násobit Ax ► to lze u řídkých matic velmi efektivně