x=[88.00001].Tato úloha má přesné řešení . Pokud ji ale perturbujeme nabude mít najednou přesné řešení . Špatná podmíněnost je způsobena tím, že matice je hodně blízko singulární matici. Geometricky si to můžeme představit tak, že když máme dvě téměř rovnoběžné přímky, poloha jejich průsečíku se bude rychle měnit i při malé změně náklonu.Pokud máme špatně podmíněnou úlohu, nepomůže nám ani zcela přesný numerický algoritmus. V takovém případě musíme opravit chybu v některém z předchozích kroků, tedy ve fyzikální a matematické formulaci.
Aritmetika s pohyblivou řádovou čárkou
Definice Nechť . Množina čísel s pohyblivou řádovou čárkou je konečná množina , kdeČíslo je mantisa a je exponent.Příklad Pro mámePříklad Standardní aritmetika IEEE, která se běžně používá v počítačích, má dvě nejčastější varianty:- Jednoduchá přesnost (32 bitů): .
- Dvojitá přesnost (64 bitů): .
Pro ušetření bitů se používá normalizace: jelikož pro nenulové číslo je první bit mantisy vždy , vynechává se. Největší a nejmenší hodnota exponentu se používají pro speciální účely, viz níže.Lemma Nechť jsou po sobě jdoucí (tedy ) se stejným exponentem . Potom .Definice Strojová přesnost čísel s pohyblivou desetinnou čárkou, značená , je vzdálenost od nejbližšího souseda v . Zaokrouhlovací jednotka je .Příklad Pro IEEE jednoduchou přesnost je .V praxi se množina rozšiřuje o subnormální čísla, k jejichž reprezentaci se používá nulový exponent:
Navíc jsou-li exponent i mantisa nulové, máme reprezentaci nuly. Dokonce máme kladnou i zápornou nulu, které se rovnají, ale nejsou shodné. Nejvyšší exponent se potom používá k reprezentaci nekonečna a neplatných hodnot.
Přesnost zobrazení reálného čísla v
Definice Nechť . Potom je nejbližší číslo z k .Věta Nechť . Potom , kde .Vliv zaokrouhlovacích chyb v aritmetice s konečnou přesností
Příklad Mějme úlohu: Spočtěte hodnotu funkce pro .Zaokrouhlíme-li na 10 míst, dostaneme . Budeme-li dále již přesně počítat s touto hodnotou, dostaneme . Pomocí metod z analýzy ale snadno ověříme, že ve skutečnosti .
Cvičení Mějme rekurenciSpočtěte s pomocí počítače . Porovnejte s přesným řešením .Řešení Spočtěme to v jazyce Janet, který používá IEEE aritmetiku s dvojitou přesností.(def c1 (/ 34 11))
(def c2 (/ 3 11))
(def a @[1 (/ 11)])
(for n 2 21
(array/push a (- (* c1 (a (- n 1))) (* c2 (a (- n 2))))))
(string/format "Výsledek: %g, přesná hodnota: %g." (a 20) (math/pow 11 -20))
Výsledek: 1.15985e-08, přesná hodnota: 1.48644e-21.
Zkusíme ještě :
Výsledek: 40.4416, přesná hodnota: 2.20949e-42.
To, co v těchto příkladech způsobuje absurdně velké chyby, se nazývá katastrofické krácení (cancellation effect), k němuž dochází, když odčítáme dvě téměř stejná čísla zatížená chybou. Mějme nějaká čísla a jejich aproximace . Potom
Tedy relativní chyba rozdílu je , což může být velké číslo, i když jsou malé.
Z toho plyne poučení, že bychom se měli vyhýbat odčítání podobně velkých čísel. Ale to se lehko řekne. Například při aproximaci derivace počítáme
Kdybychom tuto aproximaci počítali přesně, výsledek by pro menší byl lepší. Ale při počítání se zaokrouhlováním to díky zmíněnému efektu nefunguje a pod určitou mezí už se přesnost bude zhoršovat. Ukazuje se, že optimální hodnota je přibližně .
Často je ale možné se nebezpečnému odčítání nějak vyhnout. Například pokud si u příkladu s funkci rozvineme do Taylorova polynomu, dokážeme ji spočíst mnohem přesněji.
Přímá a zpětná stabilita algoritmu
Mějme například úlohu řešení soustavy lineárních algebraických rovnic . Použijeme nějaký algoritmus, například Gaussovu eliminaci, který by nám v teorii vrátil nějaký výsledek v závislosti na vstupních parametrech . Při počítání v aritmetice s konečnou přesností ale dostaneme výsledek . (Pozor, tím se nemyslí skutečná hodnota zaokrouhlená na nejbližší pohyblivé číslo, ale výsledek algoritmu, který celou dobu počítá s konečnou přesností.) Chceme najít absolutní chybu .
Mohli bychom to zkusit analyzovat přímo, ale to je pro naprostou většinu algoritmů neproveditelné, jenom u hodně jednoduchých operací jako výpočet skalárního součinu. Místo toho provedeme zpětnou analýzu chyby. Budeme hledat vstupní data taková, že řešení původní úlohy algoritmem v konečné aritmetice bylo stejné jako přesné řešení pro . Tedy má platit
Definice Algoritmus je zpětně stabilní, jestliže se chyby vzniklé v průběhu výpočtu promítnou do malé změny vstupních dat. Jinými slovy, algoritmus dává malou zpětnou chybu.Přímou chybu můžeme potom spočítat jako
Stačí už nám tedy jen odhad na citlivost úlohy.
Citlivost vlastních čísel
Mějme matice , kde vyjadřuje chybu. Zajímá nás vztah mezi a . Budeme značit
Definice je číslo podmíněnosti matice .Věta Pro každou matici platí .Důkaz Věta Schurova Pro libovolnou matici existuje unitární matice taková, že matice je horní trojúhelníková. Přitom je možné zvolit tak, že bude mít na diagonále vlastní čísla v předepsaném pořadí.Důkaz Indukcí podle . Pro volíme . Předpokládejme, že věta platí pro . Nechť pro nějaké dané vlastní číslo . Vytvoříme unitární matici tak, že doplníme na ortonormální bázi: .PotomPodle indukčního předpokladu existuje unitární matice taková, že je horní trojúhelníková. Definujme matici , u níž snadno ověříme, že je také unitární a je horní trojúhelníková.Definice Matice je normální, pokud .Věta Schurova pro normální matice Je-li normální, potom v jejím Schurově rozkladu je diagonální.Důkaz Indukcí podle . Pro volíme . Předpokládejme, že věta platí pro . Nechť je Schurův tvar . Máme je tedy také normální. To znamenáSpeciálně pro levý horní roh máme , tedy . Pro pravou dolní část máme , tedy i je normální a můžeme použít indukční předpoklad.Poznámka Snadno ověříme, že dá-li se matice zapsat jako , kde je unitární a diagonální, potom je normální. Platí tedy ekvivalence.Věta Normalizované vlastní vektory normální matice tvoří ortonormální bázi .Důkaz Podle Schurovy věty existuje rozklad , kde je diagonální. Potom sloupce jsou vlastní vektory .Definice Matice je unitární, pokud .Definice Matice je hermitovská, pokud .Věta Matice je unitární, právě když je normální a její vlastní čísla leží na jednotkové kružnici.Věta Matice je hermitovská, právě když je normální a její vlastní čísla jsou reálná.Definice Matice je diagonalizovatelná, pokud existuje regulární matice taková, že je diagonální.Věta Množina diagonalizovatelných matic je hustá v .Důkaz Nechť je Schurův rozklad nediagonalizovatelné matice. Potom na diagonále se alespoň jedno vlastní číslo opakuje. Pro libovolné zvolíme matici takovou, že a má na diagonále všechna čísla různá. Označíme-li , potom je diagonalizovatelná a platíPoznámka Díky této větě by se mohlo zdát, že při analýze citlivosti se stačí omezit na diagonalizovatelné matice a ty ostatní jimi aproximujeme. Problém je ale v tom, že malá změna prvků matice může vyvolat velkou změnu spektra.Příklad Mějme maticePotom , ale . To znamená, že při poruše řádu se mohou prvky spektra posunout o . A u velkých matic by to bylo mnohem horší.Věta o spojitosti vlastních čísel Nechť a s násobností . Potom pro každé existuje takové, že pro každou obsahuje kruh právě vlastních čísel .Poznámka Tato věta neodporuje předchozímu zjištění. To, že je funkce spojitá, nám samo o sobě nedává žádný odhad na její citlivost. (K tomu by byla potřeba lipschitzovská spojitost.)Měření vzdálenosti spekter
Budeme značit vlastní čísla matice a vlastní čísla matice .
Definice Spektrální variace matice vzhledem k matici jePoznámka Spektrální variace není metrika.Definice Hausdorffova vzdálenost matic jePoznámka Hausdorffova vzdálenost už je metrika.Definice Párová vzdálenost matic jeVěta Lemma Hadamardovo přičemž rovnost nastává, právě když má nulový sloupec nebo má vzájemně ortogonální sloupce.Důkaz Nechť , kde je unitární a je horní trojúhelníková. Jsou-li prvky na diagonále , potomCo se týče rovnosti:- (⇒)
- Platí-li rovnost, potom speciálněJe-li tento součin roven nule, potom má nulový sloupec. Jinak musí být , takže je diagonální, z čehož plyne, že má ortogonální sloupce.
- (⇐)
- Má-li nulový sloupec, potom se determinant rovná nule. Má-li ortogonální sloupce, potom je diagonální, takže v nerovnosti máme rovnost.
Věta o citlivosti spekter obecných matic Důkaz Jelikož pravá strana je symetrická vůči záměně , stačí nerovnost dokázat pro . Nechť je vlastní číslo, pro které se v definici nabývá maxima, je příslušný znormovaný vlastní vektor a jsou libovolné vektory, které ho doplní na ortonormální bázi. Označme . PotomZ Hadamardovy nerovnosti dostávámeBudeme odhadovatDosazením dostaneme kýženou nerovnost.Věta Elsner, Ostrowski Nechť . PotomVěta Bauer-Fike Nechť je regulární a . PotomDůkaz PlatíZ předpokladu levá strana rovnosti je singulární a levá část součinu na pravé straně je regulární, takže pravá část musí být singulární. Tedy existuje vektor takový, žeZ toho plynePodělením dostaneme kýženou nerovnost.Poznámka Věta platí pro libovolnou maticovou normu.Věta Geršgorinova OznačmePotomPřitom je-li sjednocení z těchto kruhů disjunktní od sjednocení ostatních, potom v něm leží přesně vlastních čísel .Důkaz Použijeme maticovou normuPoužijeme Bauerovu-Fikeovu větu (nebo lépe řečeno část jejího důkazu) s a maximovou normou. Tím dostávámePříklad Mějme maticePlatí . Z Elsnerovy věty dostávámeTo není moc dobrý odhad, jelikož spektrum se ve skutečnosti posunulo jen o . Použijeme-li Geršgorina, dostaneme odhad , což je výrazně lepší, ale pořád daleko. Můžeme to vylepšit tím, že použijeme podobnostní transformaciTato matice má stejné spektrum jako a pro různé volby můžeme z Geršgorina dostat dobrý odhad pro obě vlastní čísla. Za domácí úkol rozmyslet.Citlivost vlastních čísel diagonalizovatelných a normálních matic
Věta Nechť je diagonalizovatelná jako a je konzistentní norma splňujícíPotomDůkaz první nerovnosti Nechť . Z Bauer-Fikea mámePřitomtakžeDůsledek Nechť je normální. PotomDůkaz V takovém případě je unitární, takže .Zpětná analýza citlivosti vlastních čísel
Věta Nechť . Potom .Důkaz Poznámka K čemu je taková banální věta dobrá? Nechť jsou aproximace vlastního čísla a normovaného vlastního vektoru . Ve větě zvolíme . To nám říká, že je vlastní číslo matice , kde je reziduum.Věta Nechť je diagonalizovatelná jako , je aproximace vlastního čísla a je normovaná aproximace vlastního vektoru. Potom existuje takové, že , kde .Důkaz Označme . Podle předchozí úvahy je . Z věty o citlivosti vlastních čísel diagonalizovatelných matic mámePřitomDůsledek Speciálně pro normální matici je .Citlivost řešení soustav lineárních algebraických rovnic
Mějme soustavu , kde je regulární a . Chceme porovnat její řešení s řešením perturbované soustavy , kde , a . Rozlišíme tři případy:
- ,
- ,
- .
Věta odhad chyby pro přesnou matici Nechť . PotomDůkaz Napíšeme perturbovanou soustavu:Roznásobíme:Jelikož , máme . Z tohoPoznámka Pokud , potom není zaručeno, že při malém bude malé .Lemma Platí-lipotom je regulární.Důkaz Předpoklad je ekvivalentní tomu, že . Předpokládejme pro spor, že je singulární, tedy existuje takové, že neboli neboli . Z tohoTo je spor s předpokladem.Věta odhad chyby pro přesnou pravou stranu Nechť . Platí-lipotomDůkaz Napíšeme si perturbovanou rovnici:Roznásobením a poškrtáním dostanemeZ tohoVydělením číslem v závorce, které je podle předpokladu kladné, a přepsáním z definice dostaneme tvrzení.Poznámka Pokud , věta dává dobrý odhad chyby. Pokud je rozdíl jen malý, věta nedává smysluplný odhad. Pokud je znaménko nerovnosti opačné, věta nedává vůbec žádný odhad.Věta obecný odhad chyby Platí-lipotomDůkaz Opět si rozepíšemeÚpravami analogickými jako v důkazech předchozích vět dostaneme cílový odhad.Poznámka Je-li a , potom .Zpětná analýza citlivosti řešení soustav lineárních algebraických rovnic
Věta Rigal, Gaches Nechť je aproximace řešení soustavy . Potom existují perturbace takové, že , a platíDůkaz Označme reziduum . ZvolmePotomZ toho mámeDokážeme, že při této konkrétní volbě nastane rovnost ve znění věty:kde jsme využili rovnostiZbývá dokázat, že je to minimum. Predpokládejme pro spor, že existují nějaká s aPlatí . Odhadujme:To je podle přédpokladu menší než , což je spor.QR rozklady matic
Velká část této sekce bude opakování z Numerické matematiky 1.
Definice Flop (floating-point operation) je jedna operace s čísly s pohyblivou řádovou čárkou, přičemž jedno sčítání/odčítání může být provedeno dohromady s jedním násobením/dělením.Věta o QR rozkladu Pro libovolnou existuje takové, že je unitární, je horní trojúhelníková a .Důkaz Ke konstrukci rozkladu použijeme Givensovy rotace. Pro daný úhel můžeme sestrojit unitární matici rotace:Pro daný vektor vždy existuje takové, že . Konkrétně potřebujeme, aby platiloStačí zvolit , potomTo můžeme zobecnit na unitární matici, která pro libovolné a vynuluje -tou složku a jinak změní jen -tou. Vezmemekde volíme tak, aby . (Ve skutečnosti ji nemusíme počítat, stačí spočíst ). Nyní uvažujme celou matici . Najdeme takovou, aby . Poté najdeme takovou, aby . To nám přitom nerozhodí předtím vytvořenou nulu. Tímto způsobem vynulujeme celý první sloupec kromě první složky. Poté zvolíme takovou, aby . Iterací přes všechna vynulujeme všechno pod hlavní diagonálou. Tím dostaneme a stačí zvolitPoznámka Při praktickém využití tohoto algoritmu by mohlo dojít k podtečení. Tomu zabráníme tím, že zvolíme chytře. Je-li , položíme . Jinak:- pokud , potom
- pokud , potom
Poznámka Jaká je složitost algoritmu? Pokud to implementujeme efektivně, při nulování každého sloupce budeme mít řádově sčítání, násobení a odmocnin. Z toho celkově dostáváme flopů.Poznámka Jaká je numerická stabilita algoritmu? Ve skutečnosti spočteme nějaká , kde . Pro násobení jednou Givensovou maticí platíkde , přičemž je řádu jednotek.Definice Nechť . Potom Householderova transformace je matice .Poznámka Je-li , potomDoplníme-li na ortonormální bázi , z těchto rovností plyne, že Householderova transformace vyjadřuje zrcadlení podle nadroviny dané normálou .Věta Nechť . Potom existuje právě jedna Householderova transformace taková, že .Důkaz ZvolmeDokážeme, že splňuje tvrzení věty. Můžeme psát , kde . PotomZbývá dokázat jednoznačnost. TBDAlgoritmus QR rozklad pomocí Householderových transformací Máme matici . Najdeme Hsusholderovu matici takové, že . Potom má v prvním sloupci kromě první složky samé nuly. Dále najdeme Householderovu transformaci , jejíž normálový vektor bude mít v první složce nulu a zbytek bude udělaný tak, aby se vektor s odstraněnou první složkou zobrazil na s odstraněnou první složkou. A tak dále. Potom zvolíme .Poznámka Jaká je numerická stabilita tohoto algoritmu v porovnání s Givensovými rotacemi? Při naivním výpočtu by mohlo dojít ke katastrofickému krácení. Jeden způsob, jak se tomu vyhnout, je využít volnosti ve volbě znaménka a zvolit ho podle znaménka odpovídající složky sloupce. Například v prvním kroku budeme volitTo povede k tomu, že na diagonále budou střídavá znaménka, což může a nemusí vadit. Pokud to vadí, existuje ještě druhá možnost. Zvolíme si znaménko, jaké chceme. Pokud by mělo dojít ke katastrofickému krácení, zlomek rozšíříme opačným znaménkem.Co se týče aritmetiky s konečnou přesností, dá se ukázat, že . Je-li , potom . Z toho celkově plyne, že algoritmus je zpětně stabilní.
Poznámka Jaká je výpočetní náročnost algoritmu? Při výpočtu máme sčítání, násobení a jednu odmocninu, takže řádově flopů. Na výpočet je potřeba řádově flopů. To je celkově flopů, což je lepší než u Givensových rotací.Věta QR rozklad obdélníkové matice Nechť . Potom existuje rozklad , kde je ortogonální a splňuje , kde je horní trojúhelníková.Důkaz Doplníme na čtvercovou matici přidáním sloupců vpravo. U té provedeme QR rozklad a z něj si vezmeme jen potřebný kus.Algoritmus Gramův–Schmidtův ortogonalizační proces Nechť . Potom má lineárně nezávislých sloupců . Chceme najít ortonormální vektory takové, že .Klasický Gramův–Schmidtův algoritmus funguje tak, že položíme . Následně pro každé položíme
Algoritmus můžeme modifikovat tak, že hned na začátku každé iterace nastavíme a poté budeme pro postupně přiřazovat . Matematicky by měl vyjít stejný výsledek, ale numericky se to ukáže jako podstatná změna.Co jsme z algoritmu vlastně dostali? Platí
Pokračováním stejnou logikou dostanemeTedy označíme-li , dostaneme . Pro obdélníkovou matici to není úplně QR rozklad, protože nemusí být čtvercová, ale je to něco hodně podobného.Příklad Jaký je numerický rozdíl mezi klasickým a modifikovaným Gram-Schmidtem? Mějme maticikde , tedy . Předpokládejme pro jednoduchost, že k žádné jiné zaokrouhlovací chybě nedojde. Budeme klasickým Gramem–Schmidtem počítatVšimněme si, že , takže vektory nejsou ani přibližně ortogonální. Pokud ale použijeme modifikovaný algoritmus, ve třetím kroku dostanemePotom již bude .Poznámka Co se týče numerické stability, u klasického Grama–Schmidta se ortogonalita ztrácí rychle. Máme akorát odhadPro modifikovaného Grama–Schmidta je mnohem lepší odhad (Bjôrk, 1967):kde je polynom stupně maximálně .Poznámka Výpočetní náročnost Grama–Schmidta je flopů, což je horší než u předchozích dvou algoritmů. Klasický Gram–Schmidt se také dá paralelizovat, zatímco u ostatních to moc nejde.Metoda nejmenších čtverců
Mějme . Řešíme soustavu , tedy soustavu, která má víc rovnic než proměnných. Jelikož řešení nemusí existovat, chceme se alespoň dostat co nejblíž, tedy minimalizovat .
K řešení použijeme QR rozklad. Nechť , kde s čtvercovou. Potom
Označíme-li prvních řádků a zbylých řádků, dostáváme podle Pythagorovy rovnosti
Tento výraz nabývá minima, právě když . Tato soustava má jednoznačné řešení, takže tím dostáváme jednoznačné řešení problému nejmenších čtverců.
Metody Krylovových podprostorů
Definice Pro danou a vektor . Potom -tý Krylovův podprostor jeAlgoritmus Arnoldiho Pomocí tohoto algoritmu můžeme vygenerovat ortonormální bázi . Položíme . Následně pro každé spočtemeJe-li , hned zastavíme, jinak položíme .Poznámka Vlastně jde o specializovanou verzi modifikovaného Grama–Schmidta.Poznámka Pro Arnoldiho algoritmus lze vytvořit kompaktní schéma. Pro Arnoldiho vektory definujme matici . Dále označíme matici . Tato matice je v Hessenbergově tvaru: pod hlavní diagonálou je jedna vrstva nenulových prvků a všechno pod ní už je nulové. Platía tak dále. To můžeme kompaktně napsat jakoJe-li , což nastane nejpozději v -tém kroku, potomAlgoritmus metoda zobecněných minimálních reziduí [GMRES] Jde o iterační metodu pro řešení soustavy s . Zvolíme nějakou počáteční aproximaci . Spočteme si reziduum . Můžeme předpokládat, že , jinak máme hotovo. Pomocí Arnoldiho algoritmu najdeme ortonormální bázi . Budeme hledat řešení ve tvaru , kde . TBDAlgoritmus Lanczosův Tento algoritmus slouží pro konstrukci ortonormální báze pro symetrickou matici . Funguje to stejně jako Arnoldiho algoritmus, akorát vzniklé vektory budeme značit , aby se to nepletlo.Poznámka Kompaktní schéma je stejné jako u Arnoldiho algoritmu:Opět v posledním kroku máme , takže můžeme psátNyní využijeme předpoklad, že je symetrická> Z konstrukce plyne, že i je symetrická. Jelikož zároveň je v Hessenbergově tvaru, musí být tridiagonální. Označíme siZ kompatního schématu dostávámea tak dále. Na rozdíl od Arnoldiho algoritmu tedy stačí tříčlenné rekurence.Jak zkonstruovat koeficienty ?Koeficienty chceme zvolit tak, aby vzniklé vektory byly kolmé, což nás přinutí k volbě . Ještě musíme ověřít, že při této volbě budou kolmé i :TBDPoznámka Jelikož , algoritmus můžeme použít k výpočtu vlastních čísel. Spočítat celou matici by ale bylo náročné. Nestačila by pro nějaké ?TBDMetody Krylovových podprostorů pro řešení soustav rovnic (přehled)
Řešíme soustavu , kde je regulární. Zvolíme počáteční aproximaci a spočeteme reziduum . Poté budeme postupně nějakým způsobem hledat . K tomu jsou různé přístupy podle toho, čeho chceme dosáhnout:
- Ritzova–Galerkinova projekce
- Hledáme takové, že . Druhy: FOM (Full Orthogonalization Method), Lanczosova metoda, CG (metoda sdružených gradientů)
- Minimalizace rezidua
- Hledáme takové , pro které je minimální. Druhy: GMRES, MINRES (pro symetrické matice), ORTHODIR
- Petrovova–Galerkinova projekce
- Hledáme takové , že , kde je nějaký podprostor dimenze . Druhy: BiCG (metoda bikonjugovaných gradientů), QMR (kvazi-minimální reziduum)
- Minimalizace chyby
- Hledáme takové , aby bylo minimální na . Druhy: SYMMLQ, GMERR
Co se týče rychlosti konvergence, všechny metody konvergují rychle, pokud . Pro dosažení tohoto můžeme řešit ekvivalentní (předpodmíněnou) soustavu:
- levá
- pravá
Jak najít matici tak, aby se úloha řešila co nejlépe a zároveň se dala snadno najít?Singulární rozklady matic
Definice Nulový prostor matice je množinaa obraz jeVěta o dimenzi Lemma .Důkaz - (⊃)
- (⊂)
Lemma .Důkaz Dokážeme , zbytek známe z LAL.