PA081: Programování numerických výpočtů
7. Systémy lineárních rovnic
Aleš Křenek
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
jaro 2012
Motivace
► hledáme řešení systému
az\X\
CI22.X2.
0-2NXN
b2
&M\X\ + CLMlX-2
cimnXn = b
'M
maticové vyjádření
Ax
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 ...
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
■O0.O
2/31
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
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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
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)
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
□       4       ►   4 ■= *   < ►
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
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, 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 ...
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
■O0.O
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ů
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ší...
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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.
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ý pivotíng (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
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
13/31
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ý
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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ý pivotíng 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
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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
► 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 A-1 a následné x = A_1b
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
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
17/31
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 + l,j + 2,
1
JJ \ k=l
postup vždy využívá dříve spočtené prvky
k numerické stabilitě je třeba navíc dodat pivoting Ijj
Řídké matice
Iterační metody
■O0.O
LU dekompozice
Shrnutí a použití
řešení lineárních rovnic pro libovolný počet b
► Croutův algoritmus O (N3)
► dopředná a zpětná substituce O (N2)
► všechna b není třeba znát dopředu - hlavní přednost metody
výpočet inverzní matice
► řešení systému pro b = bázové vektory
► potřebujeme-li počítat A_1B, je lepší přímo pro sloupce B než explicitní vyjádření A-l
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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
Ří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
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
23/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 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 []
Řídké matice
Uložení v paměti
► 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
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 ® 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ě
► limfc^oo l(A')fc| = 0 pro nějakou maticovou normu
► Frobeniova norma
iai= If7~
spektrální norma
|a| = p(ara)
kde pO je maximální absolutní hotnota vlastních hodnot matice
maximální součet sloupce nebo řádku
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
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 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)
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, lze rozšířit
při výpočtu je třeba jen násobit Ax
► to lze u řídkých matic velmi efektivně
Gaussova eliminace
LU
dekompozice
Iterační zpřesnění
Řídké matice
Iterační metody
■O0.O