AdamátorZápiskyHlášky

Pokročilé partie numerické lineární algebry

Zápisky z přednášek prof. Ing. Jiřího Mikyšky, Ph.D.

  1. Úvod
    1. Obsah přednášky
      1. Citlivost úlohy
        1. Aritmetika s pohyblivou řádovou čárkou
          1. Přesnost zobrazení reálného čísla v
            1. Vliv zaokrouhlovacích chyb v aritmetice s konečnou přesností
            2. Přímá a zpětná stabilita algoritmu
              1. Citlivost vlastních čísel
                1. Měření vzdálenosti spekter
                  1. Citlivost vlastních čísel diagonalizovatelných a normálních matic
                    1. Zpětná analýza citlivosti vlastních čísel
                    2. Citlivost řešení soustav lineárních algebraických rovnic
                      1. Zpětná analýza citlivosti řešení soustav lineárních algebraických rovnic
                      2. QR rozklady matic
                        1. Metoda nejmenších čtverců
                          1. Metody Krylovových podprostorů
                            1. Metody Krylovových podprostorů pro řešení soustav rovnic (přehled)

                            Úvod

                            1. reálný problém
                            2. matematický popis (chyba modelu)
                            3. matematický model
                            4. diskretizace (chyba diskretizace)
                            5. numerický model
                            6. linearizace (chyba linearizace)
                            7. lineární algebraický problém
                            8. výpočet na počítači (zaokrouhlovací chyba)
                            9. aproximace řešení

                            Matematický model může být třeba nějaká diferenciální rovnice, která se často nedá exaktně řešit a navíc její řešení je prvek nekonečněrozměrného prostoru, takže se ani nedá popsat. Diskretizací – například pomocí metody sítí – to převedeme do konečněrozměrného prostoru. Tím ale často vznikne nelineární problém, takže ho ještě musíme linearizovat, například pomocí Newtonovy metody.

                            Předmětem tohoto kurzu bude:

                            Při výběru vhodné metody je nutné přihlížet k:

                            Obsah přednášky

                            1. Citlivost a numerická stabilita v numerické lineární algebře
                              • základní pojmy
                              • aritmetika s pohyblivou řádovou čárkou
                              • zaokrouhlovací chyby
                              • citlivost vlastních čísel matic, určení zpětné chyby
                              • citlivost řešení soustav, určení zpětné chyby
                            2. Konkrétní numerické metody
                              • QR rozklady matic
                              • metoda nejmenších čtverců
                              • Arnoldiho algoritmus
                              • Lanczosův algoritmus
                              • SVD rozklady matic

                            Webová stránka předmětu

                            Citlivost úlohy

                            Definice Citlivost úlohy popisuje vliv perturbací vstupních dat na změnu řešení úlohy.
                            Příklad Při řešení soustavy rovnic Ax=b můžeme trochu změnit matici a pravou stranu: A~A+E,b~b+e, kde EA,eb v nějaké rozumné normě. Chyby E,e jsou kombinace chyb vzniklých při měření, modelu, diskretizaci a linearizaci.
                            Definice Úloha je dobře podmíněná, pokud není příliš citlivá na malé změny vstupních parametrů.
                            Příklad U naší soustavy rovnic by to znamenalo, že pro malé E,e je x~x také malé. V praxi říkáme, že soustava je špatně podmíněná, pokud A1 má prvky řádově alespoň 105.
                            Příklad Vezměme úlohu
                            [2626.00001]x=[88.00001].
                            Tato úloha má přesné řešení x=[11]. Pokud ji ale perturbujeme na
                            [2625.99999]x=[88.00002],
                            bude mít najednou přesné řešení x=[102]. Š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ť β,t,emin,emax,β2,t1,eminemax. Množina čísel s pohyblivou řádovou čárkou je konečná množina , kde
                            {±mβet|m,e,0m<βt,emineemax}.
                            Číslo m je mantisa a e je exponent.
                            Příklad Pro β=2,t=3,emin=1,emax=3 máme
                            =±{0,1,2,3,4,5,6,7,12,32,52,72,14,34,54,74,18,38,58,78,116,316,516,716}.
                            Pří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ů): β=2,t=24,emin=125,emax=128.
                            • Dvojitá přesnost (64 bitů): β=2,t=54,emin=???,emax=???.
                            Pro ušetření bitů se používá normalizace: jelikož pro nenulové číslo je první bit mantisy vždy 1, vynechává se. Největší a nejmenší hodnota exponentu se používají pro speciální účely, viz níže.
                            Lemma Nechť y1,y2 jsou po sobě jdoucí (tedy y:y1<y<y2) se stejným exponentem e. Potom e2e1=2et.
                            Důkaz Banální.
                            Definice Strojová přesnost čísel s pohyblivou desetinnou čárkou, značená εM, je vzdálenost 1 od nejbližšího souseda v . Zaokrouhlovací jednotka je 𝓊εM2.
                            Příklad Pro IEEE jednoduchou přesnost je εM=223.

                            V praxi se množina rozšiřuje o subnormální čísla, k jejichž reprezentaci se používá nulový exponent:

                            𝒮{±mβemint|m,0<m<βt}.

                            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ť x. Potom fl(x) je nejbližší číslo z k x.
                            Věta Nechť x[min,max]. Potom fl(x)=x(1+δ), kde |δ|𝓊.

                            Vliv zaokrouhlovacích chyb v aritmetice s konečnou přesností

                            Příklad Mějme úlohu: Spočtěte hodnotu funkce f(x)1cosxx2 pro x=0.000012.

                            Zaokrouhlíme-li cosx na 10 míst, dostaneme 0.9999999999. Budeme-li dále již přesně počítat s touto hodnotou, dostaneme f(x)0.6944. Pomocí metod z analýzy ale snadno ověříme, že ve skutečnosti f(x)<0.5.

                            Cvičení Mějme rekurenci
                            a01,a1111,an3411an1311an2.
                            Spočtěte s pomocí počítače a20. Porovnejte s přesným řešením an=11n.
                            Ř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ě a40:

                            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 x,y a jejich aproximace x~=x(1+δx),y~=y(1+δy). Potom

                            x~y~=x+xδxyyδy=(xy)(1+xδxyδyxy).

                            Tedy relativní chyba rozdílu je xδxyδyxy, což může být velké číslo, i když δx,δy 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

                            f(x)=limh0f(x+h)f(x)hf(x+h)f(x)h.

                            Kdybychom tuto aproximaci počítali přesně, výsledek by pro menší h 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ě hεM.

                            Často je ale možné se nebezpečnému odčítání nějak vyhnout. Například pokud si u příkladu s 1cosx 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 Ax=b. Použijeme nějaký algoritmus, například Gaussovu eliminaci, který by nám v teorii vrátil nějaký výsledek C(A,b) v závislosti na vstupních parametrech A,b. Při počítání v aritmetice s konečnou přesností ale dostaneme výsledek fl(C(A,b)). (Pozor, tím se nemyslí skutečná hodnota C(A,b) zaokrouhlená na nejbližší pohyblivé číslo, ale výsledek algoritmu, který celou dobu počítá s konečnou přesností.) Chceme najít absolutní chybu fl(C(a,b))C(a,b).

                            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 A~,b~ taková, že řešení původní úlohy algoritmem v konečné aritmetice bylo stejné jako přesné řešení pro A~,b~. Tedy má platit

                            fl(C(A,b))=C(A~,b~).
                            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

                            fl(C(A,b))C(A,b)=C(A~,b~)C(A,b).

                            Stačí už nám tedy jen odhad na citlivost úlohy.

                            Citlivost vlastních čísel

                            Mějme matice A,EN×N, kde E vyjadřuje chybu. Zajímá nás vztah mezi σ(A) a σ(A+E). Budeme značit

                            A*A𝖳¯,
                            x,yy*x=i=1Nxiyi¯,
                            x=x2x,x=i=1N|xi|2,
                            σ(A){λ|x0:Ax=λx},
                            ρ(A)maxλσ(A)|λ|,
                            A=A2maxx=1Ax=ρ(A*A),
                            ϰ(A)AA1.
                            Definice ϰ(A) je číslo podmíněnosti matice A.
                            Věta Pro každou matici platí ϰ(A)1.
                            Důkaz 1=𝐈=AA1AA1=ϰ(A).
                            Věta Schurova Pro libovolnou matici AN×N existuje unitární matice UN×N taková, že matice RU*AU je horní trojúhelníková. Přitom je možné U zvolit tak, že R bude mít na diagonále vlastní čísla A v předepsaném pořadí.
                            Důkaz Indukcí podle N. Pro N=1 volíme U(1). Předpokládejme, že věta platí pro N. Nechť Ax=λx,x=1 pro nějaké dané vlastní číslo λ. Vytvoříme unitární matici H tak, že doplníme x na ortonormální bázi: H=(xX).Potom
                            H*AH=(x*X*)A(xX)=(x*Axx*AXX*AxX*AX)=(λx*xx*AXλX*xX*AX)(λb*0M).
                            Podle indukčního předpokladu existuje unitární matice V taková, že V*MV je horní trojúhelníková. Definujme matici U(xXV), u níž snadno ověříme, že je také unitární a U*AU je horní trojúhelníková.
                            Definice Matice AN×N je normální, pokud A*A=AA*.
                            Věta Schurova pro normální matice Je-li A normální, potom R v jejím Schurově rozkladu je diagonální.
                            Důkaz Indukcí podle N. Pro N=1 volíme U(1). Předpokládejme, že věta platí pro N. Nechť R=U*AU=(λv*0R1) je Schurův tvar A. Máme
                            R*R=U*A*UU*AU=U*A*AU=U*AA*U=U*AUU*A*U=RR*.
                            R je tedy také normální. To znamená
                            (λ0vR1*)(λv*0R1)=(λv*0R1)(λ0vR1*).
                            Speciálně pro levý horní roh máme |λ|2=|λ|+v*v, tedy v=0. Pro pravou dolní část máme vv*+R1*R1=R1R1*, tedy i R1 je normální a můžeme použít indukční předpoklad.
                            Poznámka Snadno ověříme, že dá-li se matice zapsat jako U*DU, kde U je unitární a D diagonální, potom je normální. Platí tedy ekvivalence.
                            Věta Normalizované vlastní vektory normální matice AN×N tvoří ortonormální bázi N.
                            Důkaz Podle Schurovy věty existuje rozklad A=U*ΛU, kde Λ je diagonální. Potom sloupce U* jsou vlastní vektory A.
                            Definice Matice AN×N je unitární, pokud AA*=𝐈.
                            Definice Matice AN×N je hermitovská, pokud A=A*.
                            Věta Matice AN×N je unitární, právě když je normální a její vlastní čísla leží na jednotkové kružnici.
                            Věta Matice AN×N je hermitovská, právě když je normální a její vlastní čísla jsou reálná.
                            Definice Matice AN×N je diagonalizovatelná, pokud existuje regulární matice X taková, že X1AX je diagonální.
                            Věta Množina diagonalizovatelných matic je hustá v N×N.
                            Důkaz Nechť A=U*RU je Schurův rozklad nediagonalizovatelné matice. Potom na diagonále R se alespoň jedno vlastní číslo opakuje. Pro libovolné ε>0 zvolíme matici DεN×N takovou, že Dε<ε a R+Dε má na diagonále všechna čísla různá. Označíme-li AεU*(R+Dε)U, potom Aε je diagonalizovatelná a platí
                            AAεU*DεU<ε.
                            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 matice
                            A0(0000000000000000),A1(0100001000010000),E(000000000000ε000).
                            Potom σ(A0)=σ(A1)=σ(A0+E)={0}, ale σ(A1+E)={±ε4,±𝕚ε4}. To znamená, že při poruše řádu 108 se mohou prvky spektra posunout o 102. A u velkých matic by to bylo mnohem horší.
                            Věta o spojitosti vlastních čísel Nechť AN×N a λσ(A) s násobností m. Potom pro každé ε>0 existuje δ>0 takové, že pro každou EN×N,E<δ obsahuje kruh B(λ,ε) právě m vlastních čísel A+E.
                            Důkaz Ne.
                            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 {λi} vlastní čísla matice AN×N a {λ~i} vlastní čísla matice A~A+E.

                            Definice Spektrální variace matice A~ vzhledem k matici A je
                            svA(A~)maximinj|λ~iλj|.
                            Poznámka Spektrální variace není metrika.
                            Definice Hausdorffova vzdálenost matic A,A~ je
                            hd(A,A~)max{svA(A~),svA~(A)}.
                            Poznámka Hausdorffova vzdálenost už je metrika.
                            Definice Párová vzdálenost matic A,A~ je
                            md(A,A~)minπmaxi|λ~π(i)λi|.
                            Věta
                            svA(A~)hd(A,A~)md(A,A~).
                            Lemma Hadamardovo
                            |detA|j=1NA,j,
                            přičemž rovnost nastává, právě když A má nulový sloupec nebo má vzájemně ortogonální sloupce.
                            Důkaz Nechť A=QR, kde Q je unitární a R je horní trojúhelníková. Jsou-li ρ1,,ρN prvky na diagonále R, potom
                            |detA|=|detR|=j=1N|ρj|j=1NR,j=j=1NQ*aj=j=1Naj.
                            Co se týče rovnosti:
                            (⇒)
                            Platí-li rovnost, potom speciálně
                            j=1N|ρj|=j=1NR,j.
                            Je-li tento součin roven nule, potom A má nulový sloupec. Jinak musí být R,j=ρj𝐈,j, takže R je diagonální, z čehož plyne, že A má ortogonální sloupce.
                            (⇐)
                            Má-li A nulový sloupec, potom se determinant rovná nule. Má-li ortogonální sloupce, potom R je diagonální, takže v nerovnosti máme rovnost.
                            Věta o citlivosti spekter obecných matic
                            hd(A,A~)(A+A~)11NE1N.
                            Důkaz Jelikož pravá strana je symetrická vůči záměně A,A~, stačí nerovnost dokázat pro svA(A~). Nechť λ~ je vlastní číslo, pro které se v definici nabývá maxima, x1 je příslušný znormovaný vlastní vektor a x2,,xN jsou libovolné vektory, které ho doplní na ortonormální bázi. Označme X(x1xN). Potom
                            svAN(A~)j=1N|λ~λj|=|det(Aλ~𝐈)|=|det((Aλ~𝐈)X)|.
                            Z Hadamardovy nerovnosti dostáváme
                            svAN(A~)j=1N(Aλ~𝐈)xj=(Aλ~𝐈)x1j=2N(Aλ~𝐈)xj.
                            Budeme odhadovat
                            (Aλ~𝐈)x1=Ex1E,
                            (Aλ~𝐈)xjAxj+|λ~|xjA+ρ(A~)A+A~.
                            Dosazením dostaneme kýženou nerovnost.
                            Věta Elsner, Ostrowski Nechť A,EN×N,A~A+E. Potom
                            md(A,A~)(2N1)(A+A~)11NE1N.
                            Důkaz Ne.
                            Věta Bauer-Fike Nechť QN×N je regulární a λ~σ(A~)σ(A). Potom
                            Q1(Aλ~𝐈)1Q1Q1EQ.
                            Důkaz Platí
                            Q1(A~λ~𝐈)Q=Q1(Aλ~𝐈)Q+Q1EQ=Q1(Aλ~𝐈)(I+Q1(Aλ~𝐈)1QQ1EQ).
                            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 xN{0} takový, že
                            x+Q1(Aλ~𝐈)1QQ1EQx=0.
                            Z toho plyne
                            x=Q1(Aλ~𝐈)1QQ1EQxQ1(Aλ~𝐈)1QQ1EQx.
                            Podělením Q1(Aλ~𝐈)Qx dostaneme kýženou nerovnost.
                            Poznámka Věta platí pro libovolnou maticovou normu.
                            Věta Geršgorinova Označme
                            αij=1N|Ai,j|,σi(A){ξ||ξAi,i|<αi}.
                            Potom
                            σ(A)i=1Nσi(A).
                            Přitom je-li sjednocení m z těchto kruhů disjunktní od sjednocení ostatních, potom v něm leží přesně m vlastních čísel A.
                            Důkaz Použijeme maticovou normu
                            BmaxiN^j=1NBi,j.
                            Použijeme Bauerovu-Fikeovu větu (nebo lépe řečeno část jejího důkazu) s Q𝐈,Adiag(A1,1,,AN,N),A~A a maximovou normou. Tím dostáváme
                            1(diag(A1,1,,AN,N)λ~𝐈)1(Adiag(A1,1,,AN,N))=maxiN^j=1jiN|Ai,j||Ai,iλ~|=α|Ai,iλ~|.
                            Příklad Mějme matice
                            A(1002),E(01041040).
                            Platí σ(A~)={2+108,1108}. Z Elsnerovy věty dostáváme
                            hd(A,A~)34+1081042102.
                            To není moc dobrý odhad, jelikož spektrum se ve skutečnosti posunulo jen o 108. Použijeme-li Geršgorina, dostaneme odhad λ11104,1+104,λ22104,2+104, což je výrazně lepší, ale pořád daleko. Můžeme to vylepšit tím, že použijeme podobnostní transformaci
                            (α001)(11041042)(1α001)=(1104α104α2).
                            Tato matice má stejné spektrum jako A~ 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ť A je diagonalizovatelná jako X1AX=diag(λ1,,λN)Λ a α je konzistentní norma splňující
                            diag(d1,,dN)α=maxiN^|di|.
                            Potom
                            svA(A~)X1EXαϰα(X)Eα,
                            md(A,A~)(2N1)X1EXα(2N1)ϰα(x)Eα.
                            Důkaz první nerovnosti Nechť λ~σ(A~)σ(A). Z Bauer-Fikea máme
                            X1(Aλ~𝐈)1Xα1X1EXαϰα(X)Eα.
                            Přitom
                            X1(Aλ~𝐈)1X=(X1(Aλ~𝐈)X)1=(Λλ~𝐈)1=diag(1λ1λ~,,1λNλ~),
                            takže
                            X1(Aλ~𝐈)1Xα1=(maxiN^1|λNλ~|)1=miniN~|λNλ~|=svA(A~).
                            Důsledek Nechť A je normální. Potom
                            sv(A)E,
                            md(A,A~)(2N1)E.
                            Důkaz V takovém případě je X unitární, takže ϰ(X)=1.

                            Zpětná analýza citlivosti vlastních čísel

                            Věta Nechť p<N,AN×N,Mp×p,XN×p,RAXXM,Y*p×N,Y*X=𝐈,A~ARY*. Potom A~X=XM.
                            Důkaz
                            A~X=(ARY*)X=AXRY*X=AXR=XM.
                            Poznámka K čemu je taková banální věta dobrá? Nechť λ,x jsou aproximace vlastního čísla a normovaného vlastního vektoru A. Ve větě zvolíme Xx,Y*x*,Mλ. To nám říká, že λ je vlastní číslo matice A~=Arx*, kde rAxλx je reziduum.
                            Věta Nechť ACN×N je diagonalizovatelná jako X1AX=diag(λ1,,λN), μ je aproximace vlastního čísla a y je normovaná aproximace vlastního vektoru. Potom existuje λσ(A) takové, že |λμ|ϰ(X)r, kde rAyμy.
                            Důkaz Označme A~Ary*. Podle předchozí úvahy je A~y=μy. Z věty o citlivosti vlastních čísel diagonalizovatelných matic máme
                            svA(A~)ϰ(X)ry*.
                            Přitom
                            ry*=maxx=1ry*xmaxx=1|y*x|r=r.
                            Důsledek Speciálně pro normální matici je |λμ|<r.

                            Citlivost řešení soustav lineárních algebraických rovnic

                            Mějme soustavu Ax=b, kde AN×N je regulární a bN,b0. Chceme porovnat její řešení s řešením perturbované soustavy A~x~=b~, kde A~=A+δA, x~=x+δx a b~=b+δb. Rozlišíme tři případy:

                            1. A~=A,b~b,
                            2. A~A,b~=b,
                            3. A~A,b~b.
                            Věta odhad chyby pro přesnou matici Nechť A~=A. Potom
                            δxxϰ(A)δbb.
                            Důkaz Napíšeme perturbovanou soustavu:
                            A(x+δx)=b+δb.
                            Roznásobíme:
                            Ax+Aδx=b+δb.
                            Jelikož Ax=b, máme Aδx=δb. Z toho
                            δx=A1δbA1δb,
                            δxxA1xδb=ϰ(A)Axδb<ϰ(A)δbAx=ϰ(A)δbb.
                            Poznámka Pokud ϰ(A)1, potom není zaručeno, že při malém δb bude malé δx.
                            Lemma Platí-li
                            δAA<1ϰ(A),
                            potom je A~ regulární.
                            Důkaz Předpoklad je ekvivalentní tomu, že δAA1<1. Předpokládejme pro spor, že A~ je singulární, tedy existuje x0 takové, že A~x=0 neboli Ax=δAx neboli x=A1δAx. Z toho
                            x=A1δAxA1Ax.
                            To je spor s předpokladem.
                            Věta odhad chyby pro přesnou pravou stranu Nechť b~=b. Platí-li
                            δAA<1ϰ(A),
                            potom
                            δxxϰ(A)δAA1ϰ(A)δAA.
                            Důkaz Napíšeme si perturbovanou rovnici:
                            (A+δA)(x+δx)=b.
                            Roznásobením a poškrtáním dostaneme
                            Aδx=δAxδAδx.
                            Z toho
                            δx=A1δAx+A1δAδxA1δAx+A1δAδx,
                            δx(1A1δA)A1δAx
                            Vydělením číslem v závorce, které je podle předpokladu kladné, a přepsáním z definice ϰ(A) dostaneme tvrzení.
                            Poznámka Pokud δAA1ϰ(A), 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í-li
                            δAA<1ϰ(A),
                            potom
                            δxxϰ(A)(δAA+δbb)1ϰ(A)δAA.
                            Důkaz Opět si rozepíšeme
                            (A+δA)(x+δx)=b+δb,
                            Aδx+δAx+δAδx=δb,
                            δxA1δb+A1δAx+A1δAδx.
                            Úpravami analogickými jako v důkazech předchozích vět dostaneme cílový odhad.
                            Poznámka Je-li δAA,δbb10α a ϰ(A)10β, potom δxx10βα.

                            Zpětná analýza citlivosti řešení soustav lineárních algebraických rovnic

                            Věta Rigal, Gaches Nechť x~ je aproximace řešení soustavy Ax=b. Potom existují perturbace δA,δb takové, že A~x~=b~, a platí
                            ν(x~)minδA,δbmax{δAA,δbb}=bAx~Ax~+b.
                            Důkaz Označme reziduum ιbAx~. Zvolme
                            δAAx~Ax~+bιx~*x~2,b~A~x~.
                            Potom
                            A~x=Ax~+δAx~=Ax+ιbιAx~+b=bbιAx~+b.
                            Z toho máme
                            δb=bιAx~+b.
                            Dokážeme, že při této konkrétní volbě δA,δb nastane rovnost ve znění věty:
                            δbb=ιAx~+b=bAx~Ax~+b,
                            δAA=Ax~Ax~+bιx~*x~2=bAx~Ax~+b,
                            kde jsme využili rovnosti
                            ιx~*=maxy=1ιx~*y=ιmaxy=1x~*y=ιx~.
                            Zbývá dokázat, že je to minimum. Predpokládejme pro spor, že existují nějaká δA^,δb^ s A~^x~=b~^ a
                            max{δA^A,δb^b}<ν(x).
                            Platí Ax~b=δb^δA^x~. Odhadujme:
                            ν(x~)=bAx~Ax~+b=δA^x~δb^Ax~+bδA^x+δb^Ax~+b.
                            To je podle přédpokladu menší než ν(x), 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 AN×N existuje Q,RN×N takové, že Q je unitární, R je horní trojúhelníková a A=QR.
                            Důkaz Ke konstrukci rozkladu použijeme Givensovy rotace. Pro daný úhel θ můžeme sestrojit unitární matici rotace:
                            G(cosθsinθsinθcosθ).
                            Pro daný vektor (x1x2) vždy existuje θ takové, že G𝖳(x1x2)=(y10). Konkrétně potřebujeme, aby platilo
                            sinθx1+cosθx2=0.
                            Stačí zvolit θarctan2(x2,x1), potom
                            ccosθ=x1x12+x22,ssinθ=x2x12+x22.
                            To můžeme zobecnit na unitární matici, která pro libovolné xN a ij vynuluje j-tou složku a jinak změní jen j-tou. Vezmeme
                            𝒢i,j(𝐈i100000c0s000𝐈ji1000s0c00000𝐈Nj),
                            kde volíme θ tak, aby sxi+cxj=0. (Ve skutečnosti ji nemusíme počítat, stačí spočíst c,s). Nyní uvažujme celou matici A. Najdeme 𝒢1,2 takovou, aby (𝒢1,2𝖳A)2,1=0. Poté najdeme 𝒢1,3 takovou, aby (𝒢1,3𝖳𝒢1,2𝖳A)1,3=0. 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 𝒢2,3 takovou, aby (𝒢2,3𝖳𝒢1,N𝖳𝒢1,2𝖳)3,2=0. Iterací přes všechna i<j vynulujeme všechno pod hlavní diagonálou. Tím dostaneme R a stačí zvolit
                            Qi=1N1j=i+1N𝒢i,j.
                            Poznámka Při praktickém využití tohoto algoritmu by mohlo dojít k podtečení. Tomu zabráníme tím, že c,s zvolíme chytře. Je-li xj=0, položíme c1,s0. Jinak:
                            • pokud |xj||xi|, potom
                              txixj,s11+t2,cst,
                            • pokud |xj|<|xi|, potom
                              txjxi,c11+t2,sct.
                            Poznámka Jaká je složitost algoritmu? Pokud to implementujeme efektivně, při nulování každého sloupce budeme mít řádově 2N2 sčítání, 4N2 násobení a N odmocnin. Z toho celkově dostáváme 43N3 flopů.
                            Poznámka Jaká je numerická stabilita algoritmu? Ve skutečnosti spočteme nějaká c^=c^(1+εc),s^=s(1+εs), kde |εc|,|εs|𝓊. Pro násobení jednou Givensovou maticí platí
                            fl(𝒢^i,j𝖳A)=𝒢i,j(A+E),
                            kde EA=Ku, přičemž K je řádu jednotek.
                            Definice Nechť uN,u=1. Potom Householderova transformace je matice H𝐈2uu𝖳.
                            Poznámka Je-li vN,vu, potom
                            Hu=(𝐈2uu𝖳)u=u2u(u𝖳u)=u,
                            Hv=(𝐈2uu𝖳)v=v2u(u𝖳v)=v.
                            Doplníme-li u na ortonormální bázi N, z těchto rovností plyne, že Householderova transformace vyjadřuje zrcadlení podle nadroviny dané normálou u.
                            Věta Nechť x,yN,xy,x=y. Potom existuje právě jedna Householderova transformace H taková, že Hx=y.
                            Důkaz Zvolme
                            uxyxy.
                            Dokážeme, že H𝐈2uu𝖳 splňuje tvrzení věty. Můžeme psát x=12(x+y)+12(xy), kde (x+y)(xy). Potom
                            Hx=12H(x+y)+12H(xy)=12(x+y)12(xy)=y.
                            Zbývá dokázat jednoznačnost. TBD
                            Algoritmus QR rozklad pomocí Householderových transformací Máme matici AN×N. Najdeme Hsusholderovu matici H1 takové, že H1(A,1)=±A,1𝐞1. Potom H1A má v prvním sloupci kromě první složky samé nuly. Dále najdeme Householderovu transformaci H2, jejíž normálový vektor bude mít v první složce nulu a zbytek bude udělaný tak, aby se vektor A,1 s odstraněnou první složkou zobrazil na 𝐞2 s odstraněnou první složkou. A tak dále. Potom zvolíme Qj=1NHj.
                            Poznámka Jaká je numerická stabilita tohoto algoritmu v porovnání s Givensovými rotacemi? Při naivním výpočtu u1 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 volit
                            u1A,1+A,1𝐞1sgn(A1,1)A,1+A,1𝐞1sgn(A1,1).
                            To povede k tomu, že na diagonále R 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 fl(H)H<10𝓊. Je-li fl(HA)=H(A+E), potom Ec(N)𝓊A. 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 u1 máme N sčítání, 2N násobení a jednu odmocninu, takže řádově 2N flopů. Na výpočet H1A je potřeba řádově 2N2 flopů. To je celkově 23N3 flopů, což je lepší než u Givensových rotací.
                            Věta QR rozklad obdélníkové matice Nechť AN×M,N>M. Potom existuje rozklad A=QR, kde QN×N je ortogonální a RN×M splňuje R=(RM0), kde RMM×M je horní trojúhelníková.
                            Důkaz Doplníme A na čtvercovou matici přidáním NM 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ť AN×M,NM,h(A)=M. Potom AM lineárně nezávislých sloupců a1,,aMn. Chceme najít ortonormální vektory q1,,qM takové, že [q1,,qM]λ=[a1,,aM]λ.

                            Klasický Gramův–Schmidtův algoritmus funguje tak, že položíme R1,1a1,q1a1r1,1. Následně pro každé k=2,,M položíme

                            Ri,kqi𝖳ak,waki=1k1Ri,kqi,Rk,k(w),qkwrk,k.
                            Algoritmus můžeme modifikovat tak, že hned na začátku každé iterace nastavíme wak a poté budeme pro i=1,,k1 postupně přiřazovat Ri,kqi𝖳w,wwRi,kqi. 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í

                            q1=a1R1,1a1=q1R1,1,
                            q2=a2(q1𝖳a2)q1R2,2a2=q1R1,2+q2R2,2,
                            Pokračováním stejnou logikou dostaneme
                            ak=i=1kqiRi,k.
                            Tedy označíme-li Q(q1qM), dostaneme A=QR. Pro obdélníkovou matici to není úplně QR rozklad, protože Q 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 matici
                            A(111δ000δ000δ),
                            kde 0<δ<εM, tedy fl(1+δ2)=1. Předpokládejme pro jednoduchost, že k žádné jiné zaokrouhlovací chybě nedojde. Budeme klasickým Gramem–Schmidtem počítat
                            R1,1=a1=12+δ2=˙1,
                            q1=a1=(1δ00)𝖳,
                            w2=a2R1,2q1=a2q1=(0δδ0)𝖳,
                            R2,2w2=2δ,
                            q2=w2R2,2=(012120)𝖳,
                            w3=a3R1,3q1R2,3q2=a3q1=(0δ0δ)𝖳,
                            q3=w3R3,3=(012012)𝖳.
                            Všimněme si, že q2𝖳q3=12, takže vektory nejsou ani přibližně ortogonální. Pokud ale použijeme modifikovaný algoritmus, ve třetím kroku dostaneme
                            w3=(0δ2δ2δ)𝖳.
                            Potom již bude q2𝖳q3=0.
                            Poznámka Co se týče numerické stability, u klasického Grama–Schmidta se ortogonalita ztrácí rychle. Máme akorát odhad
                            IQ^𝖳Q^1+N.
                            Pro modifikovaného Grama–Schmidta je mnohem lepší odhad (Bjôrk, 1967):
                            IQ^𝖳Q^φ(M,N)𝓊ϰ(A),
                            kde φ(M,N) je polynom stupně maximálně 2.
                            Poznámka Výpočetní náročnost Grama–Schmidta je (2N)3 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 AN×M,NM,h(A)=M,bN. Řešíme soustavu Ax=b, tedy soustavu, která má víc rovnic než proměnných. Jelikož řešení nemusí existovat, chceme se alespoň dostat co nejblíž, tedy minimalizovat |bAx|2.

                            K řešení použijeme QR rozklad. Nechť A=QR, kde R=(R^0) s R^ čtvercovou. Potom

                            a2bAx2=Q𝖳(bAx)2=Q𝖳bRx2.

                            Označíme-li c prvních M řádků Q𝖳b a d zbylých NM řádků, dostáváme podle Pythagorovy rovnosti

                            a2=cR^x2+d2.

                            Tento výraz nabývá minima, právě když R^x=c. 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 AN×N a vektor vN{0}. Potom n-tý Krylovův podprostor je
                            Kn(A,v)[v,Av,,An1v]λ.
                            Algoritmus Arnoldiho Pomocí tohoto algoritmu můžeme vygenerovat ortonormální bázi Kn(A,v). Položíme v1vv. Následně pro každé kn1^ spočteme
                            wAvk,
                            ik^:hi,kw𝖳vi,wwhi,kvi,
                            hk+1,k|w|.
                            Je-li hk+1,k=0, hned zastavíme, jinak položíme vk+1whk+1,k.
                            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 Vn(v1vn). Dále označíme matici (Hn+1,n)i,jhi,j. 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í
                            v2=Av1h1,1v1h2,1Av1=h1,1v1+h2,1v2,
                            v3=Av2h1,2v1h2,2v2h3,2Av2=h1,2v1+h2,2v2+h3,2v3
                            a tak dále. To můžeme kompaktně napsat jako
                            AVn=Vn+1Hn+1,n=VnHn,n+hn+1,nvn+1en𝖳.
                            Je-li hn+1,n=0, což nastane nejpozději v N-tém kroku, potom
                            AVn=VnHn,n.
                            Algoritmus metoda zobecněných minimálních reziduí [GMRES] Jde o iterační metodu pro řešení soustavy Ax=b s AN×N,bN. Zvolíme nějakou počáteční aproximaci x0. Spočteme si reziduum r0bAx0. Můžeme předpokládat, že r00, jinak máme hotovo. Pomocí Arnoldiho algoritmu najdeme (v1,,vn) ortonormální bázi Kn(A,r0)=[r0,Ar0,,An1r0]λ. Budeme hledat řešení ve tvaru xnx0+zn, kde znKn(A,r0). TBD
                            Algoritmus Lanczosův Tento algoritmus slouží pro konstrukci ortonormální báze Kn(A,v) pro symetrickou matici A. Funguje to stejně jako Arnoldiho algoritmus, akorát vzniklé vektory budeme značit q1,,qn+1, aby se to nepletlo.
                            Poznámka Kompaktní schéma je stejné jako u Arnoldiho algoritmu:
                            AQn=QnHn,n+hn+1,nqn+1en𝖳.
                            Opět v posledním kroku máme hn+1,n=0, takže můžeme psát
                            Hn,n=Qn𝖳AQn.
                            Nyní využijeme předpoklad, že A je symetrická> Z konstrukce plyne, že i Hn,n je symetrická. Jelikož zároveň Hn,n je v Hessenbergově tvaru, musí být tridiagonální. Označíme si
                            Hn,n(α1β200β2α20000αn1βn00βnαn).
                            Z kompatního schématu dostáváme
                            Aq1=α1q1+β2q2,Aq2=β2q1+α2q2+β3q3,Aq3=β3q2+α3q3+β4q4,
                            a tak dále. Na rozdíl od Arnoldiho algoritmu tedy stačí tříčlenné rekurence.
                            Jak zkonstruovat koeficienty αj,βj?
                            q1vv,
                            q~2Aq1α1q1,
                            q2q2q2,
                            Koeficienty αj,βj chceme zvolit tak, aby vzniklé vektory byly kolmé, což nás přinutí k volbě αj(Aqj)𝖳qj,βj|q~j+1|. Ještě musíme ověřít, že při této volbě budou kolmé i q~j+1,qj:TBD
                            Poznámka Jelikož HN,NA, algoritmus můžeme použít k výpočtu vlastních čísel. Spočítat celou matici HN,N by ale bylo náročné. Nestačila by Hn,n pro nějaké nN?TBD

                            Metody Krylovových podprostorů pro řešení soustav rovnic (přehled)

                            Řešíme soustavu Ax=b, kde AN×N je regulární. Zvolíme počáteční aproximaci x0N a spočeteme reziduum r0bAx0. Poté budeme postupně nějakým způsobem hledat xnx0+Kn(A,r0). K tomu jsou různé přístupy podle toho, čeho chceme dosáhnout:

                            Ritzova–Galerkinova projekce
                            Hledáme xn takové, že bAxnKn(A,r0). Druhy: FOM (Full Orthogonalization Method), Lanczosova metoda, CG (metoda sdružených gradientů)
                            Minimalizace rezidua
                            Hledáme takové xn, pro které je bAxn minimální. Druhy: GMRES, MINRES (pro symetrické matice), ORTHODIR
                            Petrovova–Galerkinova projekce
                            Hledáme takové xn, že bAxnSn, kde Sn je nějaký podprostor dimenze n. Druhy: BiCG (metoda bikonjugovaných gradientů), QMR (kvazi-minimální reziduum)
                            Minimalizace chyby
                            Hledáme takové xn, aby xxn bylo minimální na A𝖳Kn(A𝖳,r0). Druhy: SYMMLQ, GMERR

                            Co se týče rychlosti konvergence, všechny metody konvergují rychle, pokud A𝐈. Pro dosažení tohoto můžeme řešit ekvivalentní (předpodmíněnou) soustavu:

                            levá
                            M1Ax=M1b
                            pravá
                            AM1y=b,x=M1y
                            Jak najít matici M 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 AN,M je množina
                            N(A){xM|Ax=𝟎}.
                            a obraz je
                            R(A){Ax|xM}.
                            Věta o dimenzi
                            M=dimN(A)+dimR(A).
                            Důkaz Známe z LAL.
                            Lemma N(A𝖳A)=N(A).
                            Důkaz
                            (⊃)
                            xN(A)Ax=0A𝖳Ax=0xN(AA𝖳).
                            (⊂)
                            xN(AA𝖳)A𝖳Ax=0x𝖳A𝖳Ax=0x=0x=0xN(A).
                            Lemma h(A)=h(A𝖳)=h(AA𝖳)=h(A𝖳A).
                            Důkaz Dokážeme h(A)=h(A𝖳A), zbytek známe z LAL.