Řešení soustav lineárních rovnic Matematické algoritmy (11MAG)
Jan Přikryl 11. přednáška 11MAG pondělí 9. prosince 2013 verze:2013-12-10 15:33
Obsah 1 Soustavy
2
1.1
Formulace úlohy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2
Existence a jednoznačnost řešení . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.3
Podmíněnost soustav lineárních rovnic, odhady chyb . . . . . . . . . . . . . . . .
3
1.3.1
Vektorové normy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.3.2
Maticové normy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.3.3
Podmíněnost matic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.3.4
Odhady chyb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.3.5
Reziduum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2 Numerické metody 2.1
10
Transformace řešené úlohy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.1
Permutace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.2
Řádkové škálování . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.3
Lineární kombinace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2
Řešení soustav rovnic s trojúhelníkovou maticí . . . . . . . . . . . . . . . . . . . 12
2.3
Gaussova eliminační metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.4
LU rozklad matice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5
Algoritmus LU rozkladu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6
Výběr hlavního prvku . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3 Dodatky 3.1
22
Numerické chyby při Gaussově eliminaci . . . . . . . . . . . . . . . . . . . . . . . 22 1
3.2
Metoda LU rozkladu na počítači . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.3
Výpočetní složitost řešení lineárních soustav . . . . . . . . . . . . . . . . . . . . . 23
1
Soustavy lineárních algebraických rovnic
Mnohé matematické a počítačové modely fyzikálních a technických jevů vedou v konečné fázi na řešení soustav lineárních algebraických rovnic.Velmi často jsou tyto soustavy rozsáhlé, mají běžně statisíce nebo miliony rovnic a neznámých.Z tohoto důvodu se jak v teoretické, tak v praktické numerické matematice (tvorbě matematického softwaru) tématice numerické lineární algebry věnovala a věnuje značná pozornost a tato oblast je dnes důkladně prozkoumána. Používají se jak přímé, tak iterační metody řešení soustav lineárních rovnic. Přímé metody dávají přesné řešení soustavy, alespoň teoreticky, v konečném počtu kroků. Z tohoto hlediska je tedy řešení soustavy lineárních rovnic typická numerická úloha. Při realizaci přímé metody na počítači ovšem musíme brát v úvahu také vlastnosti počítačové aritmetiky, především zaokrouhlovací chyby. Protože velká část používaných počítačových modelů vychází z matematických modelů formulovaných typicky jako okrajové úlohy pro diferenciální rovnice a řešené soustavy lineárních rovnic jsou pouze aproximacemi těchto okrajových úloh, je často zbytečné usilovat o jejich přesné vyřešení. Stačí nám přibližné řešení získané vhodnou iterační metodou s přesností srovnatelnou s přesností aproximace teoretické okrajové úlohy. Takové přibližné řešení se dá iteračními metodami zpravidla získat mnohem úsporněji (uvažte velikost řešených soustav!) než řešení vypočítané přímou metodou. Ale i v jednotlivých krocích iteračních metod se může přímá metoda vyskytnout, navíc jsou tyto metody velmi užitečné při řešení soustav s maticemi menších rozměrů. Také teorie přímých metod, alespoň ta, kterou vyložíme zde, je podle našeho názoru snazší než je tomu u metod iteračních a patří víceméně ke všeobecnému vzdělání inženýra. V této části kurzu se proto budeme zabývat přímými metodami a ohledně iteračních metod, na které nám již nezbude čas a prostor, čtenáře musíme odkázat na literaturu. Pro detailnější obeznámení s pojmy, uváděnými níže, doporučuji i zde konzultovat knihu Michaela T. Heathe [4] nebo monografii Nicka Highama [5]. Čtenář může případně použít i některou z českých učebnic či mnoha skript o numerické matematice, která v posledních letech vyšla – například [1, 2, 3] (části skripta [3] jsou dostupné i on-line). Mnohé ze zde použitých materiálů jsme převzali právě z [4].
1.1
Formulace úlohy
Zopakujeme nejprve stručně formulaci úlohy a některé základní informace o řešení soustav a o jejich maticích, které by čtenář měl již znát z kurzu lineární algebry. V lineární algebře jsme se dozvěděli, že pohodlný způsob, jak vyjádřit lineární zobrazení (transformaci) mezi dvěma konečně-rozměrnými vektorovými prostory, spočívá v použití matic.V tomto textu se budeme zabývat pouze soustavami n rovnic o n neznámých, zmíněné vektorové prostory tedy budou mít stejnou dimenzi a příslušné matice budou čtvercové.Poznamenejme, že pod pojmem vektor zde rozumíme vždy sloupcový vektor, pokud výslovně neuvedeme jinak. V maticovém zápisu má soustava lineárních algebraických rovnic tvar Ax = b, 2
(1)
kde A je čtvercová matice řádu n (nebo také typu n × n), a x, b jsou n-složkové vektory. Taková soustava lineárních rovnic vlastně před nás klade otázku „Je možné vektor b vyjádřit jako lineární kombinaci sloupců matice A?“. Nebo, což je totéž, „Leží vektor b v prostoru span(A) = {Ax : x ∈ Rn }?“ Pokud tomu tak je, soustava má zřejmě řešení a nazývá se konzistentní. Řešení soustavy může i nemusí existovat, a pokud existuje, může i nemusí být jediné. Detaily jsou studovány v lineární algebře, ale zde se omezíme pouze na soustavy rovnic s regulárními maticemi, kde je otázka existence a jednoznačnosti řešení jednoduše vymezena. Pro jednoduchost se omezíme také pouze na soustavy lineárních algebraických rovnic s reálnými koeficienty, jejichž řešeními jsou pak nutně opět reálné vektory. Ovšem komplexní soustavy lineárních rovnic mohou být zpracovány zcela obdobně.
1.2
Existence a jednoznačnost řešení
Soustavy lineárních algebraických rovnic se čtvercovými regulárními maticemi mají vždy řešení (pro každou pravou stranu b) a toto řešení je jediné. Definice 1 (Regulární matice). Čtvercová matice A řádu n je regulární, jestliže vyhovuje kterékoli z následujících ekvivalentních podmínek: 1. k matici A existuje inverzní matice A−1 (pro takovou matici platí AA−1 = A−1 A = I, kde I je jednotková matice řádu n); 2. det(A) 6= 0 (tj. determinant matice je nenulový); 3. hodnost matice A je rovna n (hodnost matice je maximální počet lineárně nezávislých sloupců nebo řádků, které matice obsahuje); 4. soustava Az = o, kde o je nulový vektor, má pouze triviální řešení (řešením je opět pouze nulový vektor). Matice, která není regulární, se nazývá singulární. Existence a jednoznačnost řešení soustavy (1) jednoduchým způsobem závisí na regularitě matice A: Jestliže je matice A regulární, pak existuje její inverzní matice A−1 a soustava (1) má vždy právě jedno řešení x = A−1 b nezávisle na hodnotě pravé strany b. Jestliže je naproti tomu matice A singulární, pak je řešitelnost soustavy (1) závislá na pravé straně b; pro některou pravou stranu nemusí mít soustava žádné řešení, ale pokud existuje alespoň jedno řešení, pak jich existuje nekonečně mnoho. Řešení tedy v takovém případě není určeno jednoznačně. Geometricky ve dvou dimenzích každá z rovnic soustavy 2 × 2 představuje rovnici přímky a řešením soustavy je tedy společný bod obou přímek. Jsou-li tyto dvě přímky různoběžné, mají právě jeden průsečík, řešení existuje a je jediné (regulární případ), Jsou-li obě přímky rovnoběžky, pak buď nemají průsečík vůbec nebo splývají a pak je zde nekonečně mnoho společných bodů, řešení (singulární případ).
1.3
Podmíněnost soustav lineárních rovnic, odhady chyb
Poté, co jsme zhruba popsali otázky existence a jednoznačnosti řešení lineárních soustav, se nyní zaměříme na citlivost řešení x na poruchy ve vstupních datech řešené numerické úlohy, což jsou 3
zde matice A a vektor pravých stran b. Abychom takové vektorové a maticové poruchy mohli měřit, budeme potřebovat pro vektory a matice nějak zobecnit pojem velikosti, jako který nám u reálných a komplexních čísel slouží jejich absolutní hodnota. Příslušným zobecněním, které by již čtenář mohl znát z kurzu lineární algebry, je pojem normy. Pro jistotu zde příslušné informace uvedeme. 1.3.1
Vektorové normy
Budeme zde používat vektorové normy, které jsou vesměs speciálními případy tzv. p-norem, které jsou pro kladné celé číslo p > 0 a n-složkový vektor x definovány jako kxkp =
n X
!1/p p
|xi |
.
i=1
Důležité speciální případy představují • jedničková norma kxk1 =
n X
|xi |,
i=1
• euklidovská norma n X
kxk2 =
!1/2 2
|xI |
,
i=1
která odpovídá obvyklému pojmu vzdálenosti v euklidovských prostorech (nebo chcete-li, délce vektoru), a • ∞-norma (čteno nekonečno norma) kxk∞ = max |xi |, 1≤i≤n
která je limitní hodnotou p-normy při p → ∞ a říká se jí také Čebyševova norma. Při práci se všemi těmito normami dostáváme kvalitativně podobné výsledky, ovšem jednotlivé normy se můžou analyticky nebo výpočetně zpracovávat snáze nebo hůře. Při analýze podmíněnosti soustav se většinou pracuje s jedničkovou nebo Čebyševovou normou. Euklidovská norma se v této souvislosti používá méně, ale nachází účelnou aplikaci v jiných oblastech numerické lineární algebry. Dá se ukázat, že při pevně daném n se libovolné dvě z uvedených norem liší nanejvýš nevelkou multiplikativní konstantou, která sice závisí na n, ale nezávisí na konkrétním měřeném vektoru. Tím pádem jsou všechny tři normy ekvivalentní v tom smyslu, že pokud je jedna z nich malá, jsou úměrně malé i ty zbývající. V dané situaci tedy můžeme při našich úvahách volit tu z norem, s níž se nám pohodlně pracuje. Z tohoto důvodu budeme v další části tohoto textu index u norem používat pouze tam, kde to bude hrát nějakou roli, a tam, kde bude jedno, která norma se použije, budeme index vynechávat. Každá vektorová p-norma má pro libovolné dva vektory x, y tyto důležité vlastnosti: 1. kxk ≥ 0 a kxk = 0 právě pro x = o. 4
2. kαxk = |α| · kxk pro každý skalár (číslo) α. 3. kx + yk ≤ kxk + kyk (trojúhelníková nerovnost). Pro výpočet p-norem vektorů je v Matlabu k dispozici funkce norm. 1.3.2
Maticové normy
Budeme potřebovat také nějaký způsob, jak posuzovat „velikost“ matic. Jsou zde sice možné rozličné definice norem, ale ty normy, které budeme používat, budou definovány prostřednictvím již zavedených vektorových norem. Definice 2 (Maticová norma matice A). Při dané vektorové normě budeme maticovou normu matice A definovat jako kAxk kAk = max . x6=o kxk O takové maticové normě se říká, že je indukovaná příslušnou vektorovou normou nebo že je této normě podřízená. Intuitivně se na tuto definici můžeme dívat tak, že norma matice měří to, jak transformace danou maticí může maximálně protáhnout ten který vektor, když délku vektorů měříme použitou vektorovou normou. S některými maticovými normami se pracuje snáze než s jinými. Příklad 3 (Maticová norma odpovídající jedničkové vektorové normě). Tato norma se rovná jednoduše maximální hodnotě sloupcových součtů prvků matice braných v absolutních hodnotách, kAk1 = max j
n X
|aij |.
i=1
Příklad 4 (Maticová norma odpovídající Čebyševově vektorové normě). Maticová norma odpovídající Čebyševově vektorové normě je prostě maximální řádkový součet prvků matice braných v absolutních hodnotách, kAk∞ = max i
n X
|aij |.
j=1
Maticová norma odpovídající euklidovské vektorové normě se dá definovat pomocí vlastních čísel matic, její explicitní výpočet je tudíž složitější a my se jí zda zabývat nebudeme. Takto definované maticové normy mají pro libovolné matice A, B tyto důležité vlastnosti: 1. kAk ≥ 0 a kAk = 0 právě když A = O (nulová matice). 2. kαAk = |α| · kAk pro každý skalár (číslo) α. 3. kA + Bk ≤ kAk + kBk. 4. kABk ≤ kAk · kBk. 5. kAxk ≤ kAk · kxk pro každý vektor x.
5
Poznamenejme, že obecně se za maticovou normu považuje jakákoli funkce matic, která má první tři právě uvedené vlastnosti. Normy, které mají navíc poslední dvě vlastnosti, se pak nazývají submultiplikativní nebo konzistentní. Zdůrazňujeme znovu, že maticové normy indukované vektorovými p-normami jsou s nimi vždy konzistentní, mají tedy všech pět výše uvedených vlastností. Námi popsané maticové normy se v Matlabu dají snadno počítat pomocí funkce norm. 1.3.3
Podmíněnost matic
Číslo podmíněnosti regulární čtvercové matice A se v dané maticové normě definuje jako cond(A) = kAk · kA−1 k. Pro singulární matici A se úmluvou klade cond(A) = ∞. Při hlubším studiu se dá ukázat, že takto zavedené číslo podmíněnosti je konzistentní s definicí podanou v Přednášce 8, a to v tom smyslu, že ohraničuje poměr mezi relativní změnou řešení soustavy lineárních rovnic a danou relativní změnou ve vstupních datech. Z definice je zřejmé, že hodnota čísla podmíněnosti závisí na použité maticové normě, což se někdy označuje tím, že se použije příslušný index; píšeme pak například cond1 (A) nebo cond∞ (A). Z toho, co jsme říkali v Odstavci 1.3.1 ale vyplývá, že hodnoty čísla podmíněnosti se při různých normách mohou lišit maximálně nevelkou multiplikativní konstantou (která závisí na n) a jsou tedy jako měřítka podmíněnosti stejně užitečné. Jako zajímavost uvádíme, že se dá ukázat, že číslo podmíněnosti také měří poměr mezi maximálním relativním protažením a maximálním relativním zkrácením libovolných dvou n-složkových vektorů poté, co na ně aplikujeme matici A. Z definice čísla podmíněnosti matic lze vcelku snadno odvodit následující jeho důležité vlastnosti, které platí při použití jakékoli námi zavedené normy: 1. Pro všechny matice A platí cond(A) ≥ 1. 2. Pro jednotkovou matici je cond(I) = 1. 3. Pro jakoukoli matici A a každé číslo α různé od nuly je cond(αA) = cond(A). 4. Pro každou diagonální matici D = diag(d) platí cond(D) = (maxi |di |)/(minj |dj |). Číslo podmíněnosti je mírou toho, jak blízko je daná matice k tomu, aby byla singulární: matice s velkým číslem podmíněnosti (co to znamená, řekneme v následujícím odstavci) je téměř singulární, zatímco matice, jejíž číslo podmíněnosti je řádu jednotek, má k singularitě daleko. Z definice je také zřejmé, že regulérní matice a její inverze mají stejné číslo podmíněnosti. Všimněme si také toho, že přestože determinant singulární matice je nulový, není hodnota determinantu sama o sobě měřítkem toho, jak blízko nebo daleko má daná matice A k singulární matici. Vezměme si jako příklad matici 0,1 I, jejíž determinant je roven 0,1n , což už u malé matice s n = 30 dává jako hodnotu determinantu 10−30 , tedy velmi malé číslo. Přitom je tato matice velmi dobře podmíněná, její číslo podmíněnosti je 1 v jakékoli normě. Brzy uvidíme, že užitečnost čísla podmíněnosti se ukáže při posuzování přesnosti řešení soustav lineárních rovnic. V definici čísla podmíněnosti ovšem vystupuje inverzní matice, takže počítat 6
jeho hodnotu by byla obvykle netriviální záležitost. Ve skutečnosti by výpočet čísla podmíněnosti z jeho definice vyžadoval podstatně více výpočetní práce než samotné řešení posuzované soustavy lineárních rovnic. V praxi lze naštěstí jako vedlejší produkt při řešení soustav získat alespoň dobrý odhad čísla podmíněnosti, souhlasící s jeho skutečnou hodnotou přinejmenším řádově. V Matlabu je pro tento účel k dispozici funkce cond a některé další funkce. Příklad 5 (Odhad čísla podmíněnosti). Uvažujme matici "
A=
0,913 0,659 0,457 0,330
#
.
Odhadneme její čísla podmíněnosti v Matlabu; použijeme k tomu funkci cond. Dostaneme cond1 (A) = cond∞ (A) ≈ 1,6958 · 104 , cond2 (A) ≈ 1,2485 · 104 . Odhady se řádově shodují. Vzhledem k nevelkému řádu matice a počtu uvedených platných cifer můžeme danou matici považovat za špatně podmíněnou.Jak uvidíme v příštím odstavci, pokud by prvky této matice byly naměřeny s chybou velikosti řádově 10−4 , nemohli bychom při řešení soustav s touto maticí zaručit ve výsledku ani jednu platnou číslici. 1.3.4
Odhady chyb
Kromě toho, že je číslo podmíněnosti spolehlivým indikátorem blízkosti k singulární matici, nám dává také užitečné kvantitativní odhady chyby pro vypočítaná řešení lineárních soustav.Uvedeme zde pouze výsledky, jakkoli matematika za nimi skrytá není nikterak obtížná. Případné zájemce v tomto směru odkazujeme například na [4, 5]. Uvažujme nejprve poruchy v pravé straně soustavy lineárních rovnic. Nechť x je řešení regulární ˆ je řešením soustavy Aˆ soustavy lineárních rovnic Ax = b a nechť x x = b + ∆b s porušenou ˆ − x. Pak se dá odvodit odhad pravou stranou. Označme ∆x = x k∆xk k∆bk ≤ cond(A) . kxk kbk Číslo podmíněnosti matice je tedy jakýsi „amplifikační faktor“, který pomáhá odhadnout maximální relativní změnu řešení způsobenou danou relativní poruchou ve vektoru pravých stran (porovnejte to s obecně zavedeným pojmem čísla podmíněnosti z Přednášky 8). Podobný výsledek se dá odvodit pro relativní změny v prvcích matice A. Jestliže Ax = b a (A + E)ˆ x = b, pak k∆xk kEk ≤ cond(A) . kˆ xk kAk Obecně se pak dá říci, že při současných změnách prvků matice i pravé strany platí (až na veličiny vyšších řádů) k∆xk k∆bk kEk ≤ cond(A) + . kxk kbk kAk Vidíme tedy opět, že relativní změna řešení je ohraničená součinem čísla podmíněnosti a relativní změny v datech úlohy.
7
Obrázek 1: Dobře podmíněná (vlevo) a špatně podmíněná (vpravo) soustava dvou lineárních rovnic o dvou neznámých.
Geometrická interpretace popsaných výsledků ve dvou dimenzích je ta, že jsou-li přímky definované dvěma rovnicemi téměř rovnoběžné, pak pokud mohou být tyto přímky zatíženy zaokrouhlovacími chybami nebo chybou měření, není jejich průsečík definován ostře. Pokud ale tyto přímky vůbec nejsou rovnoběžné, jsou třeba téměř na sebe kolmé, je jejich průsečík relativně ostře vymezen. Oba tyto případy ilustrujeme na Obrázku 1, kde čárkované přímky vyznačují oblast nepřesnosti pro každou plně vytaženou přímku, takže průsečík skutečných přímek (které přesně neznáme) by mohl být kdekoli ve vystínovaném rovnoběžníku. Velké číslo podmíněnosti je tedy spojeno s velkou nepřesností získaného řešení. Máme-li tedy shrnout to, jsme se zatím naučili, pak v případě, že jsou vstupní data přesná na strojovou přesnost (relativní chyba v Matlabu tedy cca 10−16 ), můžeme za rozumný odhad relativní chyby ve vypočítaném řešení považovat kˆ x − xk < ≈ cond(A) mach kxk kde mach je relativní přesnost aritmetiky počítače (angl. unit roundoff, definici a výpočet máte popsány v podkladech pro 8. cvičení a také v [4] dopsat do přednášky ).Jednoduchý způsob interpretace získaných výsledků je ten, že vypočítané řešení soustavy ztrácí během řešení zhruba log10 (cond(A)) desítkových číslic přesnosti vůči přesnosti vstupních dat. Podíváme-li se na matici z Příkladu 5, vidíme, že její číslo podmíněnosti je řádově 104 , takže při řešení soustavy dvou rovnic s touto maticí nemůžeme v výsledku očekávat jedinou správnou číslici, pokud prvky matice nejsou přesné na více než čtyři platné číslice a pokud řešení nepočítáme v aritmetice používající alespoň čtyři platné desítkové číslice. Jako kvantitativní míra citlivosti tak číslo podmíněnosti matice hraje při řešení soustav lineárních rovnic stejnou roli, jako tomu bylo u obecně zavedeného čísla podmíněnosti v Přednášce 8. Důležitý rozdíl je zde ale v tom, že maticové číslo podmíněnosti nemůže nikdy být menší než 1. 1.3.5
Reziduum
Jedním ze způsobů, jak ověřit řešení rovnice, je dosadit je do řešené rovnice a podívat se, jak se ˆ lineární soustavy po dosazení shoduje její levá a pravá strana. Reziduum přibližného řešení x rovnic Ax = b je rozdíl r = b − Aˆ x. Pokud je A regulární, platí pro chybu teoreticky, že k∆xk = kˆ x − xk = 0 tehdy a jen tehdy, je-li krk = 0. V praxi však nemusí být obě tyto veličiny malé současně. Všimněme si především toho, 8
že vynásobíme-li soustavu Ax = b libovolnou nenulovou konstantou, její řešení se nezmění, ale reziduum samé se násobí stejným číslem. Reziduum tedy můžeme udělat libovolně velké nebo malé podle toho, jakým faktorem soustavu vynásobíme (přeškálujeme ji, změníme měřítko), a tím pádem je velikost rezidua nesmyslné kritérium, pokud ji nějak nevztáhneme k velikosti dat úlohy a řešení. Z tohoto důvodu se zavádí pojem relativní reziduum, které se definuje jako R(ˆ x) =
krk . kAk · kˆ xk
Abychom uvedli relativní reziduum do vztahu s chybou získaného přibližného řešení, všimneme si, že k∆xk = kˆ x − xk = kA−1 (Aˆ x − b)k = k − A−1 rk ≤ kA−1 k · krk. Vydělíme-li obě strany získané nerovnosti kˆ xk a použijeme definici čísla podmíněnosti matice A, dostáváme odsud odhad k∆xk krk ≤ cond(A) = cond(A) · R(ˆ x). kˆ xk kAk · kˆ xk Malé relativní reziduum tedy znamená malou relativní chybu vypočítaného řešení tehdy a jen tehdy, když je A dobře podmíněná matice. Pro algoritmy přímých metod řešení soustav lineárních rovnic se dá ukázat, že jimi vypočítané ˆ je přesným řešením nějaké porušené soustavy řešení x (A + E)ˆ x = b. Popravdě řečeno je takových porušených soustav celá řada, stačí si uvědomit, že při daných A ˆ máme ve výše uvedeném vztahu určit n2 prvků matice E z n lineárních rovnic. Ukazuje se ax však, že mezi poruchami E se dá vybrat porucha s minimální maticovou normou, a té se pak říká zpětná chyba vypočítaného řešení (na rozdíl od přímé chyby ∆x). Odhady zpětné chyby jsou pro algoritmy přímých metod k dispozici a můžeme u nich tedy říci, že jimi vypočítané numerické řešení je přesným řešením původní soustavy porušené o poruchu, jejíž velikost (v normě) umíme odhadnout shora. Výše provedené úvahy o relativním reziduu nám snadno ale dávají dolní odhad takové zpětné chyby, Je totiž krk = kb − Aˆ xk = kEˆ xk ≤ kEk · kˆ xk, odkud plyne po vydělení normou A nerovnost krk kEk ≤ . kAk · kˆ xk kAk Velké relativní reziduum tak implikuje velkou zpětnou chybu, což znamená, že algoritmus použitý k výpočtu řešení nebyl stabilní a získané řešení nemůže být kvalitní. Na druhé straně, stabilní algoritmy nutně produkují malá relativní rezidua bez ohledu na podmíněnost úlohy, a tím pádem malé reziduum vypovídá o kvalitě získaného řešení jen málo. Příklad 6 (Malé reziduum). Uvažujme soustavu lineárních rovnic "
Ax =
0,913 0,659 0,457 0,330
#"
9
x1 x2
#
"
=
0,254 0,127
#
= b,
s jejíž maticí jsme se již setkali v Příkladu 5. Řekněme, že máme k dispozici dvě přibližná řešení "
ˆ1 = x
−0,0827 0,5
#
"
0,999 −1,001
#
.
a
ˆ2 = x
a
kr2 k1 = 2,4 · 10−2 .
Normy odpovídajících reziduí jsou kr1 k1 = 2,1 · 10−4
Které z obou řešení tedy máme považovat za přesnější? Díky mnohem menšímu reziduu bychom ˆ 1 . Ale přesné řešení této soustavy je (dá se to mohli mít sklon k tomu, abychom řekli, že je to x snadno ověřit) vektor " # 1 x= , −1 ˆ 2 je ve skutečnosti mnohem přesnější než x ˆ 1 . Důvodem této možná překvapivé skutečnosti takže x je, že matice A je špatně podmíněná, jak jsme ostatně viděli v Příkladu 5. Díky jejímu velkému číslu podmíněnosti zde malé reziduum neznamená, že máme přesné řešení.
2
Numerické metody pro řešení soustav lineárních rovnic
Přímé numerické metody pro řešení soustav lineárních algebraických rovnic, kterými se nyní budeme zabývat, nejsou vlastně ničím jiným, než zobecněním postupů známých ze střední školy. Tam jsme u soustavy dvou rovnic o dvou neznámých tím či oním způsobem soustavu upravili: z jedné rovnice jsme vyloučili (eliminovali) druhou neznámou a převedli ji tak na rovnici pro jedinou neznámou. Tuto rovnici jsme vyřešili, výsledek dosadili do některé z původních rovnic a z ní pak vypočítali zbývající neznámou. Ukazuje se, že všechny takové úpravy se dají přehledně zapsat ve formě maticových transformací. Numerické metody, které zde budeme popisovat, vlastně provádějí na danou soustavu (její matici) určitou posloupnost úprav, které se dají chápat jako maticové transformace a jejichž cílem je danou soustavu převést na soustavu jednodušší (s maticí speciálního tvaru), která je s původní soustavou ekvivalentní, ale řeší se podstatně přímočařeji.
2.1
Transformace řešené úlohy
Prozatím zde pouze popíšeme maticové transformace, které nemění řešení soustavy, a ukážeme, že k nim mimo jiné patří ty transformace, které si lze představit jako elementární operace s rovnicemi soustavy, speciálně jako vynásobení rovnice nenulovou konstantou nebo jako přičtení násobku jedné rovnice k rovnici jiné. Naším cílem bude posléze převést řešenou soustavu (její matici) ne takový tvar, v němž ji půjde už snadno a přímo vyřešit. Abychom nebyli příliš tajemní, poznamenáváme již zde, že naším cílem bude získat ekvivalentní soustavu, jejíž matice bude trojúhelníková. Je především jasné, že pokud obě strany soustavy Ax = b s regulární maticí vynásobíme zleva jakoukoli regulární maticí M, na řešení soustavy to nebude mít vliv. Abychom se o tom přesvědčili, všimněme si, že řešení transformované soustavy MAz = Mb je dáno jako z = (MA)−1 Mb = A−1 M−1 Mb = A−1 b = x. 10
2.1.1
Permutace
Důležitým příkladem regulární transformace je skutečnost, že řádky matice A a odpovídající složky ve vektoru b lze přeházet, aniž bychom změnili řešení x. Intuitivně je to zřejmé: všechny rovnice soustavy musí být splněny současně, takže na pořadí, v jakém je zapíšeme, nezáleží. Formálně se takové přehazování řádků dá v maticovém tvaru zapsat jako násobení obou stran soustavy zleva permutační maticí P, což je čtvercová matice (opět řádu n), mající v každém řádku a sloupci právě jednu jedničku a všechny ostatní prvky rovné nule (je to vlastně jednotková matice s přeházenými řádky a sloupci). Podívejme se, jak permutační matice může například přeházet řádky sloupcového vektoru:
v3 v1 0 0 1 1 0 0 v2 = v1 . v2 v3 0 1 0 Permutační matice je vždy regulární; ve skutečnosti je její inverzí prostě její transpozice, P−1 = PT . Pro pořádek poznamenejme, že transpozicí matice M, kterou značíme MT , rozumíme matici, jejímiž sloupci jsou řádky matice M. Je-li tedy N = MT , platí nij = mji . 2.1.2
Řádkové škálování
Jiným jednoduchým, ale důležitým případem regulární maticové transformace je diagonální škálování. Připomínáme, že matice D je diagonální, jestliže v ní pro všechna i 6= j platí dij = 0, což znamená, že jedinými prvky, které v ní mohou být nenulové, jsou prvky na hlavní diagonále dii , i = 1, . . . , n. Determinantem diagonální matice je součin jejích prvků na hlavní diagonále, takže regulární diagonální matice má na hlavní diagonále pouze nenulové prvky. Vynásobíme-li obě strany soustavy Ax = b zleva regulární diagonální maticí, DAx = Db je to totéž, jako bychom každý z řádků matice soustavy a pravé strany vynásobili odpovídajícím diagonálním prvkem dii , a mluvíme zde proto o řádkovém škálování. Řádky, kterým v použité diagonální matici odpovídají na hlavní diagonále jedničky, se přitom nemění. Řádkové škálování může být v praxi někdy užitečné, protože nemění řešení soustavy a přitom může pozitivně ovlivnit chování použitého algoritmu na počítači.
2.1.3
Lineární kombinace řádků matice
Další ekvivalentní úprava soustavy lineárních rovnic spočívá v tom, že provedeme určitou lineární kombinaci řádků soustavy a přičteme ji k některé jiné rovnici. Speciální případ, který zde popíšeme, protože hraje významnou roli v numerických metodách řešení soustav, je ten, kdy určitý násobek jednoho řádku matice přičteme k jinému řádku (nebo jej od něj odečteme). Ukážeme, že taková operace se dá opět zapsat jako transformace jistou regulární maticí, takže nemění řešení soustavy (koneckonců se to zdá být zřejmé intuitivně). Především je vhodné si zopakovat definici násobení matic a uvědomit si, že prvek na místě ij ve výsledné matici není nic jiného než skalární součin i-tého řádku prvního činitele s j-tým sloupcem činitele druhého. Pak 11
už bychom měli snadno pochopit, že přičtení m-násobku řádku i v dané matici k jejímu j-tému řádku se dá maticově zapsat tak, že danou matici vynásobíme maticí M, která se od jednotkové matice liší pouze tím, že má navíc na místě ji číslo m. Matice M je tedy trojúhelníková matice, která má například pro i < j tvar podobný 1 ··· . . .. .. 0 ··· M= 0 ··· .. . . . .
0 ···
0 ··· .. . . . . 1 0 ··· m 1 ··· .. .. . . . . . 0 0 ··· 0 .. .
0 .. . 0 0 .. .
.
1
Jestliže je v uvedené matici prvek m ve sloupci k a řádku k + 1, pak MA realizuje přičtení m-násobku řádku k v matici A k řádku k + 1 v téže matici. Příklad 7 (Lineární kombinace řádků matice). Mějme transformační matici M a matici A zadané jako 1 0 0 1 2 3 M = 0 1 0 a A = 0 −1 1 . 0 m 1 −1 2 1 Potom
1 2 3 0 −1 1 MA = −1 + 0m 2 − 1m 1 + 1m Matice podobného typu (třeba i s více nenulovými prvky m v některém sloupci) budeme v dalším nazývat elementární eliminační matice. Všimněme si, že elementární eliminační matice mají na diagonále samé jedničky, takže jejich determinant je roven jedné a matice jsou regulární. Popsané kombinování řádků soustavy tedy nemění její řešení.
2.2
Řešení soustav rovnic s trojúhelníkovou maticí
První otázka, kterou si při konstrukci přímých metod řešení soustav lineárních rovnic klademe je, jaké soustavy leze řešit nejsnáze. Na ně se pak budeme obecné soustavy snažit ekvivalentními úpravami převést. Představme si, že máme soustavu Ax = b, v níž se najde rovnice, která obsahuje jenom jednu neznámou (jednu složku hledaného řešení). To znamená, že v daném řádku matice A je jednom jeden nenulový prvek. Z takové rovnice se pak dá snadno (prostým dělením) vypočítat příslušná složka řešení. Teď si představme, že v soustavě je ještě jiná rovnice, v níž vystupují pouze dvě neznámé, a že jedna z nich je ta, kterou jsme právě vypočítali. Jestliže nyní do této druhé rovnice dosadíme hodnotu právě vypočítané složky řešení, můžeme z ní opět snadno vypočítat tu druhou. Pokud je struktura matice i nadále podobná, to jest pokud nám v každé rovnici přibude pouze jedna nová neznámá, dají se tak postupně snadno vypočítat všechny složky řešení. Matici s touto speciální vlastností se říká trojúhelníková, a to z důvodu, který bude brzy zřejmý. Protože soustavy lineárních rovnic s trojúhelníkovými maticemi se snadno řeší, jsou právě ony cílem úprav prováděných v přímých numerických metodách na obecných soustavách rovnic. 12
Pro výpočetní účely se ukazuje jako užitečné zavést trojúhelníkové matice dvou speciálních tvarů. Řekneme, že matice L je dolní trojúhelníková, jestliže všechny její prvky nad hlavní diagonálou jsou nulové (tj. `ij = 0 pro i < j). Podobně řekneme, že matice U je horní trojúhelníková, jestliže všechny její prvky pod hlavní diagonálou jsou nuly (tj. uij = 0 pro i > j. Poznamenáváme, že trojúhelníkovou matici definovanou výše v obecném smyslu můžeme vždy permutacemi řádků nebo sloupců převést na dolní nebo horní trojúhelníkovou matici. Příklad 8 (Horní a dolní trojúhelníková matice).
1 0 0 L = 0 0 0 . a −1 0 1
1 2 3 U = 0 −1 1 0 0 1
Jakkoli je způsob řešení soustav s trojúhelníkovými maticemi nyní již snad zřejmý, považujeme za užitečné algoritmy pro jejich řešení popsat podrobněji. Proces řešení dolní trojúhelníkové soustavy Lx = b postupným dosazováním se nazývá přímá substituce a matematicky se dá zapsat jako x1 = b1 /`11 ,
xi = bi −
i−1 X
`ij xj /`ii ,
i = 2, . . . , n.
j=1
Formální zápis příslušného algoritmu může vypadat například následujícím způsobem. Algoritmus 1 Přímá substituce pro soustavu s dolní trojúhelníkovou maticí for j = 1 to n do {cyklus přes sloupce} if `jj = 0 then {stop při singulární matici} return error end if xj ← bj /`jj {výpočet složky řešení} for i = j + 1 to n do bi ← bi − `ij xj {aktualizace pravé strany} end for end for Podobně se u horní trojúhelníkové soustavy Ux = b proces řešení postupným dosazováním nazývá zpětná substituce a matematicky se dá zapsat jako
xn = bn /un ,
xi = bi −
n X
uij xj /uii ,
i = n − 1, . . . , 1.
j=i+1
Uvádíme jeden ze způsobů algoritmické realizace zpětné substituce. V obou těchto algoritmech jsme zvolili pořadí indexů v cyklech tak, že se k prvkům matic přistupuje po sloupcích (a ne po řádcích) a že operace, kterou provádí vnitřní cyklus, je skalár krát vektor plus vektor, operace zkracovaná saxpy (z angl. scalar times a vector plus a vector) a používaná v základních knihovnách numerické lineární algebry. Ukládání maticí po sloupcích je typické pro některé programovací jazyky, jako je například Matlab nebo Fortran. Naproti tomu třeba v jazyku C se matice ukládají po řádcích a při větších soustavách by výše popsané algoritmy bylo třeba přepsat tak, aby se v nich k maticím přistupovalo po řádcích. To nepředstavuje zásadní problém, zmiňujeme se o tom pouze proto, abychom čtenáře upozornili, že způsob implementace matematického algoritmu může podstatně ovlivnit jeho efektivitu, a to v závislosti na použitém programovacím jazyku a výpočetním systému. Poznamenáváme ještě, že výskyt 13
Algoritmus 2 Zpětná substituce pro soustavu s horní trojúhelníkovou maticí for j = n to 1 do {zpětný cyklus přes sloupce} if ujj = 0 then {stop při singulární matici} return error end if xj ← bj /ujj {výpočet složky řešení} for i = 1 to j − 1 do bi ← bi − uij xj {aktualizace pravé strany} end for end for nulového prvku na diagonále bude znamenat selhání uvedených algoritmů, což ovšem není překvapivé, protože trojúhelníková matice s nulovým prvkem na diagonále má nulový determinant a je tedy singulární. Příklad 9 (Soustava rovnic s trojúhelníkovou maticí). Mějme soustavu rovnic s horní trojúhelníkovou maticí x1 3 1 2 2 0 −1 −6 x2 = −6 . x3 1 0 0 1 Z poslední rovnice je dáno x3 = 1. Tuto hodnotu pak můžeme dosadit do druhé rovnice a z ní vypočítat −x2 −6 = −6, x2 = 0. Nakonec do první rovnice dosadíme jak x3 tak x2 a vypočítáme z ní x1 + 0 + 2 = 3, x1 = 1.
2.3
Gaussova eliminační metoda
Naší strategií nyní je navrhnout regulární lineární transformaci, která bude převádět danou obecnou soustavu lineárních algebraických rovnic na ekvivalentní soustavu s trojúhelníkovou maticí, již pak budeme moci snadno řešit postupnými substitucemi. Budeme tedy potřebovat transformaci, která by nahrazovala vybrané nenulové prvky dané matice nulami (vyeliminovala je). Uvidíme, že se toho dá dosáhnout prováděním vhodných lineárních kombinací řádků dané matice. Základem přímých numerických metod pro řešení obecných soustav lineárních algebraických rovnic je tak postup, který nyní popíšeme na příkladu soustavy čtyř rovnic o čtyřech neznámých a který se nazývá Gaussova eliminační metoda. Uvažujme tedy soustavu lineárních rovnic a11 x1 + a12 x2 + a13 x3 + a14 x4 = b1 a21 x1 + a22 x2 + a23 x3 + a24 x4 = b2 a31 x1 + a32 x2 + a33 x3 + a34 x4 = b3 a41 x1 + a42 x2 + a43 x3 + a44 x4 = b4 Položíme mi1 = ai1 /a11 ,
14
i = 2, 3, 4,
a odečteme mi1 krát první rovnici od i-té rovnice pro i = 2, 3, 4. Dostaneme tak upravenou ekvivalentní soustavu a11 x1 + a12 x2 + a13 x3 + a14 x4 = b1 a022 x2 + a023 x3 + a024 x4 = b02 a032 x2 + a033 x3 + a034 x4 = b03 a042 x2 + a043 x3 + a044 x4 = b04 kde a0ij = aij − mi1 a1j
a
b0i = bi − mi1 b1 .
Vidíme, že neznámá x1 byla z posledních tří rovnic vyeliminována. Protože při eliminaci se čísly mi1 násobí první rovnice říká se jim multiplikátory (z angl. multiply, násobit). Nyní položíme mi2 = a0i2 /a022 ,
i = 3, 4,
a odečteme mi2 krát druhou rovnici od i-té rovnice (4i = 3, 4). Výsledkem je soustava a11 x1 + a12 x2 + a13 x3 + a14 x4 = b1 a022 x2 + a023 x3 + a024 x4 = b02 a0033 x3 + a0034 x4 = b003 a0043 x3 + a0044 x4 = b004 kde a00ij = a0ij − mi2 a02j
a
b00i = b0i − mi2 b02 .
A nakonec položíme mi3 = a00i3 /a0033 ,
i = 4,
a odečteme mi3 krát třetí rovnici od rovnice čtvrté. Výsledkem je soustava s horní trojúhelníkovou maticí a11 x1 + a12 x2 + a13 x3 + a14 x4 = b1 a022 x2 + a023 x3 + a024 x4 = b02 a0033 x3 + a0034 x4 = b003 000 a000 44 x4 = b4
kde 00 00 a000 ij = aij − mi3 a3j
a
00 00 b000 i = bi − mi3 b33 .
Protože má tato soustava horní trojúhelníkovou matici, lze ji snadno vyřešit zpětnou substitucí.Postupu, kterým jsme původní soustavu lineárních rovnic převedli na ekvivalentní soustavu s trojúhelníkovou maticí se také říká přímý chod a řešení výsledné trojúhelníkové soustavy zpětný chod Gaussovy eliminace.Prvkům matic, které postupně objevují na hlavní diagonále, tedy v našem případě a11 , a022 , a0033 , a000 44 , se říká hlavní prvky nebo také pivoty. Popsaný postup je proveditelný pouze, pokud jsou při něm všechny hlavní prvky nenulové. Pokud během eliminace narazíme na nulový hlavní prvek, postupujeme tak, že rovnici s tímto prvkem prohodíme z některou ze zbývajících rovnic, která má v daném sloupci nenulový koeficient. Pokud by ovšem v daném sloupci matice už byly všechny prvky pod diagonálou nulové 15
(včetně hlavního prvku), není těžké ukázat, že by to znamenalo, že matice řešené soustavy je singulární (a tedy by soustava buď neměla řešení nebo by jich měla nekonečně mnoho). Gaussovu eliminační metodu ilustrujeme na jednoduchém konkrétním příkladu řešení soustavy tří rovnic o třech neznámých. Příklad 10 (Gaussova eliminace). Budeme řešit soustavu lineárních rovnic x1 + 2 x2 + 2 x3 = 3, 4 x1 + 4 x2 + 2 x3 = 6, 4 x1 + 6 x2 + 4 x3 = 10, která se dá maticově zapsat jako
3 x1 1 2 2 Ax = 4 4 2 x2 = 6 = b. 10 x3 4 6 4 Pro účely Gaussovy eliminace je vhodné zapsat data úlohy ve formě tak zvané rozšířené matice soustavy 1 2 2 3 6 4 4 2 4 6 4 10 a postup Gaussovy eliminace provádět na této matici; je to tak přehlednější. Po prvním kroku eliminace, během kterého odečteme čtyřnásobek prvního řádku rozšířené matice od druhého a třetího řádku, rozšířená matice přejde na
1 2 2 3 −6 0 −4 −6 . 0 −2 −4 −2 Abychom nyní vynulovali člen pod diagonálou v druhém sloupci získané rozšířené matice, odečteme druhý řádek vynásobený jednou polovinou od třetího řádku a dostaneme tak matici
1 2 2 3 0 −4 −6 −6 . 1 0 0 −1 To je již rozšířená matice trojúhelníkové soustavy rovnic, z níž můžeme zpětnou substitucí dostat hledané řešení jako −1 x = 3 . −1 Algoritmus, který jsme výše obecně popsali pro soustavy čtyř rovnic, se dá snadno zobecnit pro soustavy libovolného řádu. Ukážeme nyní, že výše odvozený postup Gaussovy eliminační metody se dá elegantně zapsat pomocí maticových transformací, v nichž opět figurují trojúhelníkové matice.
16
2.4
LU rozklad matice
Ukážeme opět na příkladu obecné soustavy 4 × 4, že během Gaussovy eliminace vlastně konstruujeme dolní trojúhelníkovou matici L a horní trojúhelníkovou matici U takové, že A = LU. Položíme A1 = A a dále
M1 =
1 −m21 −m31 −m41
0 1 0 0
0 0 1 0
0 0 0 1
,
kde čísla mij jsou multiplikátory zavedené v předchozím odstavci. Pak se dá snadno ověřit (vzpomeňte si na Odstavec 2.1.3), že platí
A2 = M1 A1 =
a11 0 0 0
a12 a022 a032 a042
a13 a023 a033 a043
a14 a024 a034 a044
,
což je matice soustavy, která výše vznikla po prvním kroku Gaussovy eliminační metody. Nyní položíme 1 0 0 0 0 1 0 0 M2 = . 0 −m32 0 0 0 −m42 0 1 Opět nyní není těžké ukázat, že
A3 = M2 A2 =
a11 a12 0 a022 0 0 0 0
a13 a023 a0033 a0043
a14 a024 a0034 a0044
,
což je matice soustavy, kterou jsme výše obdrželi po provedení druhého kroku eliminace. Konečně položíme 1 0 0 0 0 1 0 0 M3 = . 0 0 1 0 0 0 −m43 1 Potom
U = M3 A3 =
a11 a12 a13 0 a022 a023 0 0 a0033 0 0 0
a14 a024 a0034 a000 44
je horní trojúhelníková matice výsledné soustavy získané Gaussovou eliminační metodou. Jak jsme již poznamenali v Odstavci 2.1.3, matice Mk , které se během právě popsaného postupu vyskytly, se v literatuře nazývají elementární eliminační matice. Zapíšeme-li vše maticově, vidíme, že pro horní trojúhelníkovou matici U platí U = M3 M2 M1 A. Položíme-li tedy −1 −1 L = M−1 1 M2 M3 ,
17
pak A = LU. Protože inverze dolní trojúhelníkové matice je opět dolní trojúhelníková a protože součin dolních trojúhelníkových matic je opět dolní trojúhelníková matice, je výše sestrojená matice L dolní trojúhelníková matice, a ukázali jsme tak tedy, že Gaussova eliminační metoda vlastně realizuje LU rozklad matice A na součin dolní a horní trojúhelníkové matice. Této skutečnosti se hojně využívá při realizaci algoritmů Gaussovy eliminace na počítači. Dosud jsme se nijak nezabývali výpočtem prvků matice L. Jejich stanovení je však překvapivě jednoduché, neboť se dá ukázat, že matice L má na hlavní diagonále jedničky a pod diagonálou má v pozici `ij , i > j, multiplikátor mij . Není tak těžké to ověřit. Poznamenáváme opět, že zde popsané úvahy o trojúhelníkovém rozkladu matice nezávisejí na řádu matice a dají se snadno zobecnit pro obecné n.
2.5
Algoritmus LU rozkladu
Pokud již máme k dispozici rozklad A = LU, můžeme snadno převést úlohu řešit soustavu Ax = b na řešení dvou soustav rovnic s trojúhelníkovými maticemi. Je to okamžitě vidět z řešené soustavy, kam za A dosadíme její trojúhelníkový rozklad. Nejprve vyřešíme soustavu Lz = b, jejíž řešení z je vlastně eliminací transformovaná pravá strana původní soustavy, a pak soustavu Ux = z, jejímž řešením je řešení původní soustavy rovnic. Poznamenáváme ještě, že pokud jsme pro danou matici jednou stanovili její trojúhelníkový rozklad, můžeme tímto způsobem řešit s menší pracností i více soustav rovnic, které se liší pouze pravou stranou. To bývá někdy v praxi užitečné. Zbývá nám pro případné zájemce popsat algoritmus výpočtu matic L a U. Na počítači zpravidla prvky těchto matic ukládáme postupně na místa původně obsazená prvky matice A. Co do počtu míst to vychází, protože víme, že diagonální prvky matice L jsou jedničky a nemusíme je tedy ukládat. Příslušný algoritmus má řadu variant, uvádíme zde jednu z možností. Algoritmus 3 LU rozklad Gaussovou eliminací for k = 1 to n − 1 do {cyklus přes sloupce} if akk = 0 then return {stop, je-li hlavní prvek nulový} end if for i = k + 1 to n do {výpočet multiplikátorů pro daný sloupec} aik ← aik /akk end for for j = k + 1 to n do for i = k + 1 to n do {aplikujeme transformaci na zbývající submatici} aij ← aij − aik akj end for end for end for
18
2.6
Výběr hlavního prvku
Gaussova eliminační metoda v sobě nese jeden zřejmý problém, o kterém jsme se již zmínili, ale také jeden další, méně zřetelný problém. Možným zřejmým problémem je to, že celý proces selže, pokud se během eliminace vyskytne nulový hlavní prvek, protože by se při výpočtu multiplikátorů v odpovídajícím sloupci dělilo nulou. Ale řešení tohoto problému je téměř stejně zřejmé: pokud je v kroku k eliminace na diagonále nulový hlavní prvek, prohodíme řádek k v soustavě (tedy v její matici i v její pravé straně, jinak řečeno, řádek k v rozšířené matici soustavy) s nějakým z pod ním ležících řádků, jehož prvek ve sloupci k je nenulový. Z Odstavce 2.1.1 již víme, že tím se řešení soustavy nezmění. S nenulovým diagonálním prvkem jako hlavním prvkem pak celý postup pokračuje jako obvykle. Takové přehazování řádků se pak nazývá výběr hlavního prvku nebo také pivotace. Příklad 11 (Pivotace a singularita). Případná nutnost provádět výběr hlavního prvku a přehazovat řádky nemá nic společného s tím, zda daná matice je singulární či nikoliv. Tak například matice " # 0 1 A= 1 0 je zřejmě regulární (spočítejte si determinant), ale pokud nepřehodíme řádky, nemá žádný LU rozklad. Naproti tomu singulární matice "
A=
1 1 1 1
#
má (bez přehazování řádků) LU rozklad "
A=
1 1 1 1
#
"
=
1 0 1 1
#"
1 1 0 0
#
= LU.
Ale co dělat v případě, že na hlavní diagonále ani pod ní ve sloupci k není žádný nenulový prvek? Pak v daném kroku není třeba provádět nic, protože všechny prvky, které bychom měli eliminovat jsou nulové, a přejdeme tedy prostě k dalšímu sloupci (je tedy Mk = I). Všimněme si toho, že tento krok zanechá na diagonále nulový prvek, takže výsledná matice U bude singulární, nicméně LU rozklad lze stejně dokončit. Znamená to ovšem, že s takovou maticí selže zpětný chod, protože se v něm vyskytuje dělení každým z diagonálních prvků matice U, ale to není až tolik překvapivé, protože v takovém případě je nutně singulární už původní matice A (stačí spočítat si její determinant z determinantů jejích trojúhelníkových faktorů, což jsou součiny jejich diagonálních prvků). Oním méně zřetelným problémem je to, že v počítačové aritmetice nám na diagonále nemusí vyjít přesná nula, ale pouze nějaké velmi malé číslo, což nás vede k následující úvaze. Jako hlavní prvek se při výpočtu multiplikátorů dá v zásadě využít jakékoli nenulové číslo, ale v počítačové aritmetice s konečnou přesností bychom se měli zajímat také o to, jak minimalizovat šíření zaokrouhlovacích chyb během výpočtu. Konkrétně je naším cílem omezit velikosti multiplikátorů, takže se již existující zaokrouhlovací chyby vzniklé v předchozích krocích eliminace nebudou zesilovat při následujících transformacích realizovaných elementárními eliminačními maticemi. Absolutní hodnota multiplikátorů zaručeně nikdy nepřevýší jedničku, jestliže v každém sloupci vybereme za hlavní prvek ten z prvků na diagonále nebo pod ní, který má největší absolutní hodnotu. Takové strategii se říká částečný výběr hlavního prvku a v praxi je tento výběr pivota podstatný pro to, abychom mohli sestrojit numericky stabilní implementaci 19
Gaussovy eliminace pro obecné matice. Poznamenejme pro pořádek, že slovo částečný je zde použito proto, že vhodný hlavní prvek hledáme pouze v prvním sloupci dosud nevyeliminované submatice. Mohli bychom jej vybírat i v celé této submatici, pak bychom ovšem možná museli permutovat i některé sloupce. Navíc tento způsob volby pivota (úplný výběr hlavního prvku) nepřináší žádné podstatné výhody. Ukážeme nyní na extrémním, ale jednoduchém příkladu, jak se může špatná volba pivota při výpočtu projevit. Podobný efekt bychom mohli pozorovat s méně přehnanou volbou pivota u větších soustav či matic. Příklad 12 (Malé hlavní prvky). Pracujeme-li v počítačové aritmetice, která má pouze konečnou přesnost, musíme se vyhýbat nejen nulovým hlavním prvkům, ale také hlavním prvkům, které jsou malé vzhledem k ostatním prvkům daného sloupce, a to proto, abychom zabránili nepřijatelnému růstu chyb. Ukážeme to na velmi zjednodušeném příkladu, kde symbol ≈ znamená výsledek počítačové operace. Označíme zde jakékoli kladné číslo, které je menší než strojové epsilon mach , tedy jakékoli dostatečně malé číslo x, pro které platí 1 + x ≈ 1. Pokud chce čtenář být konkrétní, může vzít v jednoduché aritmetice = 10−8 a v dvojnásobné přesnosti, se kterou normálně pracuje Matlab, vzít = 10−16 . Vezměme nyní matici "
#
1 1 1
A=
.
To je regulární matice, jejíž determinant je zhruba roven −1 a která je dobře podmíněná, má číslo podmíněnosti řádu jednotek. Pokud nepřehodíme řádky, bude hlavním prvkem , takže výsledný multiplikátor je −1/ a " # 1 0 L= . 1/ 1 Jako výsledná horní trojúhelníková matice po eliminaci nám na počítači vyjde "
U=
#
1 0 1 − 1/
"
1 0 −1/
≈
#
.
Vynásobíme-li ale nyní matice L a U získané na počítači, dostaneme "
1 1 0
LU =
#
6= A.
Použití malého prvku a tím pádem velkého multiplikátoru způsobilo velkou ztrátu informace ve výsledné matici U. Přehodíme-li řádky, bude nyní hlavní prvek roven 1 a odpovídající multiplikátor bude −, takže "
1 0 1
L=
#
.
Jako výsledná horní trojúhelníková matice nám na počítači tentokrát vyjde "
U=
#
1 1 0 1−
"
≈
1 1 0 1
Vynásobením takto získaných matic L a U nyní dostaneme "
LU =
1 1 1
#
což je přesně původní matice A s přehozenými řádky. 20
,
#
.
Dá se ukázat, že přehazování řádků během přímého chodu Gaussovy eliminace se dá shrnout jako provedení jedné permutační transformace na původní matici a že tedy na rozdíl od přímého LU rozkladu matice A, který nemusí vždy existovat (viz výše), ke každé čtvercové matici A existuje permutační matice P tak, že k permutované matici lze LU rozklad nalézt, tedy PA = LU. Algoritmus LU rozkladu popsaný v předchozím odstavci se při výběru hlavních prvků příliš nemění, je třeba si pouze uvědomit, že nyní místo soustavy Ax = b řešíme ekvivalentní soustavu PAx = Pb. Po získání LU rozkladu matice PA (matice P vlastně vzniká teprve během tohoto rozkladu) nejprve řešíme soustavu Lz − Pb s dolní trojúhelníkovou maticí přímou substitucí a poté řešíme zpětnou substitucí soustavu Ux = z s horní trojúhelníkovou maticí. V numerické lineární algebře se ukazuje, že existují třídy matic, pro které je výběr hlavního prvku zbytečné provádět a kde tedy výše uvedená permutační matice je jednotková. Patří se například pozitivně definitní symetrické matice, s nimiž se často setkáváme v odborné praxi. Pro nedostatek místa zde explicitně algoritmus LU rozkladu s výběrem hlavního prvku neuvádíme. Příklad 13 (Gaussova eliminace s částečnou pivotací). V Příkladu 10 jsme nepoužívali permutace řádků a v důsledku toho byly některé multiplikátory větší než jedna. Jako ilustraci nyní tento příklad zopakujeme, ale budeme tentokrát používat částečný výběr hlavního prvku. Soustava, kterou jsme řešili v Příkladu 10, je
1 2 2 x1 3 Ax = 4 4 2 x2 = 6 = b 4 6 4 x3 10 a její rozšířená matice je
1 2 2 3 6 . 4 4 2 4 6 4 10 Největší prvek v prvním sloupci je 4, takže prohodíme první dva řádky a permutovaná rozšířená matice bude 4 4 2 6 3 . 1 2 2 4 6 4 10 Abychom vynulovali prvky pod diagonálou v prvním sloupci, odečteme od druhého řádku 0,25 krát první řádek a od třetího řádku první řádek krát jedna. jako transformovanou matici tak dostaneme 4 4 2 6 0 1 1,5 1,5 . 0 2 2 4 Poslední prvek pod diagonálou v druhém sloupci je 2, takže přehodíme poslední dva řádky a dostaneme tak opět permutovanou matici
6 4 4 2 4 . 0 2 2 0 1 1,5 1,5
21
Abychom vynulovali subdiagonální prvek ve druhém sloupci, odečteme od třetího řádku 1/2násobek druhého řádku a dostaneme tak konečný tvar vyeliminované matice jako
6 4 4 2 4 . 0 2 2 0 0 0,5 −0,5 Převedli jsme tak původní soustavu rovnic na ekvivalentní soustavu rovnic s horní trojúhelníkovou maticí, kterou můžeme nyní vyřešit zpětnou substitucí a dostáváme tak stejně jako v Příkladu 10 −1 x = 3 . −1 Pro zájemce uvádíme ještě výsledek ve formě LU rozkladu. Výsledná permutační matice je
0 1 0 P= 0 0 1 1 0 0 a máme
4 4 2 1 0 0 1 0 0 2 2 = LU. PA = 1 0 0 0,5 0,25 0,5 1
3 3.1
Dodatky Numerické chyby při Gaussově eliminaci
Jak jsme již uvedli v Odstavci 1.3.5, platí pro relativní reziduum vypočítaného řešení nerovnost krk kEk ≤ , kAk · kˆ xk kAk kde E je zpětná chyba v matici A. Ale jak velká asi bude kEk při praktických výpočtech? V numerické lineární algebře se ukazuje, že pro LU rozklad Gaussovou eliminací platí odhad tvaru kEk ≤ ρ n mach , kAk kde ρ se nazývá růstový faktor a mach je strojové epsilon charakterizující přesnost počítačové aritmetiky. V Matlabu je standardně strojové epsilon k dispozici jako proměnná eps a je řádové velikosti 10−16 . Pokud jde o růstový faktor, je to poměr největšího prvku matice U k největšímu prvku matice A braný v absolutních hodnotách. Pokud neprovádíme výběr hlavního prvku, může být tento faktor ρ libovolně velký, a tím pádem není Gaussova eliminační metoda bez pivotace obecně numericky stabilní (nemusí poskytovat malou zpětnou chybu). Pokud provádíme částečný výběr hlavního prvku, může být růstový faktor teoreticky velikosti 2n−1 (protože v nejhorším možném případě se velikost prvků může v každém kroku eliminace zdvojnásobit), ale takové chování Gaussovy eliminační metody je v praxi nesmírně vzácné. V prakticky řešených příkladech dochází při částečném výběru hlavních
22
prvků pouze k malému nebo žádnému růstu mezivýsledků, takže se dá přibližně brát jako odhad zpětné chyby kEk < ≈ n mach . kAk Tento vztah znamená, že při řešení soustav lineárních rovnic Gaussovou eliminací s částečným výběrem hlavního prvku téměř vždy dostáváme velmi malá relativní rezidua, a to bez ohledu na to, jak špatně podmíněná řešená soustava je. Malé relativní reziduum tak nutně neznamená, že jsme získali přesný výsledek; to platí pouze u dobře podmíněných soustav. Doporučujeme čtenáři, aby se v této souvislosti znovu podíval na Příklad ??.
3.2
Metoda LU rozkladu na počítači
Existuje několik způsobů jak implementovat metodu LU rozkladu na počítači, které se liší v prvé řadě tím, jak se přistupuje k prvkům matice A (po řádcích nebo po sloupcích). Účinnost těchto implementací se na daném počítači může výrazně lišit a jedna a ta samá implementace se může na různých výpočetních systémech chovat různě. Některé z variant Gaussovy eliminační metody si dokonce vysloužily vlastní pojmenování, takže máme například Croutovu metodu nebo Doolittlovu metodu. Speciální varianta metody LU rozkladu pro řešení soustav rovnic s třídiagonálními maticemi bývá v inženýrské praxi často nazývána Thomasovým algoritmem. Jakkoli se ale různé implementace LU rozkladu mohou významně lišit co do výkonnosti, produkují pro danou matici A v podstatě tentýž výsledný rozklad. Pokud jde o hospodaření s pamětí, je běžné, že vytvářené trojúhelníkové faktory postupně přepisují původní prvky matice A. Prostorově to přesně vychází, protože matice L má na diagonále jedničky a ty není třeba ukládat. Aby se minimalizoval fyzický pohyb dat, neprovádí se přehazování řádků při výběru hlavních prvků explicitně. Místo toho řádky zůstávají na svých původních místech, ale k zaznamenávání jejich pořadí se používá pomocný celočíselný vektor indexů. Je totiž zřejmé, že konečný efekt všech záměn řádků je pouze nějaká permutace čísel 1, . . . , n. V Matlabu je algoritmus metody LU rozkladu s částečným výběrem hlavního prvku (včetně řešení obou trojúhelníkových soustav) zabudován do operátoru \ (zpětné lomítko, angl. backslash). Řešení soustavy rovnic Ax = b zapisujeme v Matlabu jako x = A\b. Matlab má také k dispozici funkci inv pro výpočet inverzní matice. Výslovně ale upozorňujeme, že počítat řešení soustavy jako x = inv(A)*b je zcela nesprávné. Takovýto způsob výpočtu je totiž podstatně pracnější než použití zpětného lomítka a je také horší z hlediska numerické stability, zejména u větších matic. Zápis x = A−1 b je tedy užitečný pouze v teoretických úvahách na papíře.
3.3
Výpočetní složitost řešení lineárních soustav
Výpočet LU rozkladu matice řádu n Gaussovou eliminací vyžaduje zhruba n3 /3 násobení v pohyblivé řádové čárce a zhruba stejný počet sčítání. Řešení vzniklých trojúhelníkových soustav pro jednu pravou stranu si vyžádá asi n2 násobení a podobný počet sčítání. Jestliže tedy řád matice n roste, začne v celkové výpočetní práci stále více převládat konstrukce LU rozkladu. Z hlediska výpočetní složitosti je výpočet inverzní matice samotné asi třikrát náročnější než řešení soustavy s jednou pravou stranou. Dalším postupem, kterému je třeba se při výpočtech řešení soustav vyhnout, je Cramerovo pravidlo, v němž se každá složka řešení počítá jako podíl dvou determinantů. Pro větší matice (už jistě od řádu dvacet) je tento postup přímo astrono23
micky náročný na výpočetní práci a Cramerovo pravidlo se tak dá použít buď u matic velmi malého řádu nebo v teoretické lineární algebře.
Reference [1] MÍKA, Stanislav. MVŠT. Numerické metody algebry. Praha: SNTL, 1982, 169 s. [2] MÍKA, Stanislav. MVŠT. Numerické metody algebry. 2. opravené vyd. Praha: SNTL, 1985, 169 s. [3] MÍKA, Stanislav a Marek BRANDNER. Numerické metody I. Plzeň: FAV ZČU, 2000, 169 s. [4] HEATH, M. T. Scientific Computing: An Introductory Survey. 2nd Edition. New York: McGraw-Hill, 2002, 563 s. [5] HIGHAM, N. J. Accuracy and Stability of Numerical Algorithms. 2nd Edition. Philadelphia: SIAM, 2002, 680 s.
24
