AdamátorZápiskyHlášky

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

Ú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ů
    • Garoldiho 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 ‖E‖≪‖A‖,‖e‖≪‖b‖ 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 A−1 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=[10−2]. Š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,t≥1,emin≤emax. Množina čísel s pohyblivou řádovou čárkou je konečná množina ℱ⊂ℝ, kde
ℱ≔{±m⋅βe−t|m,e∈ℤ,0≤m<βt,emin≤e≤emax}.
Čí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 e2−e1=2e−t.
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=2−23.

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

𝒮≔{±m⋅βemin−t|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)≔1−cosxx2 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
a0≔1,a1≔111,an≔3411an−1−311an−2.
Spočtěte s pomocí počítače a20. Porovnejte s přesným řešením an=11−n.
Ř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⋅δx−y−y⋅δy=(x−y)⋅(1+x⋅δx−y⋅δyx−y).

Tedy relativní chyba rozdílu je x⋅δx−y⋅δyx−y, 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)=limh→0f(x+h)−f(x)h≈f(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.

Také se to může objevit při řešení kvadratické rovnice. Počítáme-li diskriminant b2−4ac, může se při nevhodném vstupu stát, že hodnoty b2 a 4ac budou podobné (konkrétně když kořeny rovnice jsou blízko sobě). V takovém případě může být lepší použít jiný vzoreček.

Často je ale možné se nebezpečnému odčítání nějak vyhnout. Například pokud si u příkladu s 1−cosx 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,E∈ℂN×N, kde E vyjadřuje chybu. Zajímá nás vztah mezi σ(A) a σ(A+E). Budeme značit

A*≔A𝖳¯,
⟨x,y⟩≔y*x=∑i=1Nxiyi¯,
‖x‖=‖x‖2≔⟨x,x⟩=∑i=1N|xi|2,
σ(A)≔{λ∈ℂ|∃x≠0:Ax=λx},
ρ(A)≔maxλ∈σ(A)|λ|,
‖A‖=‖A‖2≔max‖x‖=1‖Ax‖=ρ(A*A),
ϰ(A)≔‖A‖⋅‖A−1‖.
Definice ϰ(A) je číslo podmíněnosti matice A.
Věta Pro každou matici platí ϰ(A)≥1.
Důkaz 1=‖𝐈‖=‖A⋅A−1‖≤‖A‖⋅‖A−1‖=ϰ(A).
Věta Schurova Pro libovolnou matici A∈ℂN×N existuje unitární matice U∈ℂN×N taková, že matice R≔U*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 A∈ℂN×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 A∈ℂN×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 A∈ℂN×N je unitární, pokud A⋅A*=𝐈.
Definice Matice A∈ℂN×N je hermitovská, pokud A=A*.
Věta Matice A∈ℂN×N je unitární, právě když je normální a její vlastní čísla leží na jednotkové kružnici.
Věta Matice A∈ℂN×N je hermitovská, právě když je normální a její vlastní čísla jsou reálná.
Definice Matice A∈ℂN×N je diagonalizovatelná, pokud existuje regulární matice X taková, že X−1AX 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í
‖A−Aε‖≤‖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 10−8 se mohou prvky spektra posunout o 10−2. A u velkých matic by to bylo mnohem horší.
Věta o spojitosti vlastních čísel Nechť A∈ℂN×N a λ∈σ(A) s násobností m. Potom pro každé ε>0 existuje δ>0 takové, že pro každou E∈ℂN×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 A∈ℂN×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=1N‖A⋅,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=1N‖R⋅,j‖=∏j=1N‖Q*aj‖=∏j=1N‖aj‖.
Co se týče rovnosti:
(⇒)
Platí-li rovnost, potom speciálně
∏j=1N|ρj|=∏j=1N‖R⋅,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~‖)1−1N⋅‖E‖1N.
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≔(x1⋯xN). Potom
svA(A~)N≤∏j=1N|λ~−λj|=|det(A−λ~𝐈)|=|det((A−λ~𝐈)⋅X)|.
Z Hadamardovy nerovnosti dostáváme
svA(A~)N≤∏j=1N‖(A−λ~𝐈)⋅xj‖=‖(A−λ~𝐈)x1‖⋅∏j=2N‖(A−λ~𝐈)⋅xj‖.
Budeme odhadovat
‖(A−λ~𝐈)x1‖=‖Ex1‖≤‖E‖,
‖(A−λ~𝐈)xj‖≤‖Axj‖+|λ~|⋅≤‖A‖+ρ(A~)≤‖A‖+‖A~‖.
Dosazením dostaneme kýženou nerovnost.
Věta Elsner, Ostrowski Nechť A,E∈ℂN×N,A~≔A+E. Potom
md(A,A~)≤(2N−1)⋅(‖A‖+‖A~‖)1−1N⋅‖E‖1N.
Důkaz Ne.
Věta Bauer-Fike Nechť Q∈ℂN×N je regulární a λ~∈σ(A~)∖σ(A). Potom
‖Q−1(A−λ~𝐈)−1Q‖−1≤Q−1EQ.
Důkaz Platí
Q−1(A~−λ~𝐈)Q=Q−1(A−λ~𝐈)Q+Q−1EQ=Q−1(A−λ~𝐈)⋅(I+Q−1(A−λ~𝐈)QQ−1EQ).
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 x∈ℂN∖{0} takový, že
x+Q−1(A−λ~𝐈)−1QQ−1EQx=0.
Z toho plyne
‖x‖=‖Q−1(A−λ~𝐈)−1QQ−1EQx‖≤‖Q−1(A−λ~𝐈)−1Q‖⋅‖Q−1EQ‖⋅‖x‖.
Podělením ‖Q−1(A−λ~𝐈)Q‖⋅‖x‖ dostaneme kýženou nerovnost.
Poznámka Věta platí pro libovolnou maticovou normu.
Věta Geršgorinova Označme
αi≔∑j=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
‖B‖∞≔maxi∈N^∑j=1NBi,j.
Použijeme Bauerovu-Fikeovu větu (nebo lépe řečeno část jejího důkazu) s Q≔𝐈,A≔diag(A1,1,…,AN,N),A~≔A a maximovou normou. Tím dostáváme
1≤‖(diag(A1,1,…,AN,N)−λ~𝐈)−1(A−diag(A1,1,…,AN,N))‖∞=maxi∈N^∑j=1j≠iN|Ai,j||Ai,i−λ~|=α|Ai,i−λ~|.
Příklad Mějme matice
A≔(1002),E≔(010−410−40).
Platí σ(A~)={2+10−8,1−10−8}. Z Elsnerovy věty dostáváme
hd(A,A~)≤3⋅4+10−8⋅10−4≈2⋅10−2.
To není moc dobrý odhad, jelikož spektrum se ve skutečnosti posunulo jen o 10−8. Použijeme-li Geršgorina, dostaneme odhad λ1∈⟨1−10−4,1+10−4⟩,λ2∈⟨2−10−4,2+10−4⟩, což je výrazně lepší, ale pořád daleko. Můžeme to vylepšit tím, že použijeme podobnostní transformaci
(α001)⋅(110−410−42)⋅(1α001)=(110−4α10−4α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 X−1AX=diag(λ1,…,λN)≕Λ a ‖⋅‖α je konzistentní norma splňující
‖diag(d1,…,dN)‖α=maxi∈N^|di|.
Potom
svA(A~)≤‖X−1EX‖α≤ϰα(x)⋅‖E‖α,
md(A,A~)≤(2N−1)⋅‖X−1EX‖α≤(2N−1)⋅ϰα(x)⋅‖E‖α.
Důkaz první nerovnosti Nechť λ~∈σ(A~)∖σ(A). Z Bauer-Fikea máme
‖X−1(A−λ~𝐈)−1X‖α−1≤‖X−1EX‖α≤ϰα(X)⋅‖E‖α.
Přitom
X−1(A−λ~𝐈)−1X=(X−1(A−λ~𝐈)X)−1=(Λ−λ~𝐈)−1=diag(1λ1−λ~,…,1λN−λ~),
takže
‖X−1(A−λ~𝐈)−1X‖α−1=(maxi∈N^1|λN−λ~|)=mini∈N~|λN−λ~|=svA(A~).
Důsledek Nechť A je normální. Potom
sv(A)≤‖E‖,
md(A,A~)≤(2N−1)⋅‖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,A∈ℂN×N,M∈ℂp×p,X∈ℂN×p,R≔AX−XM,Y*∈ℂp×N,Y*X=𝐈,A~≔A−RY*. Potom A~X=XM.
Důkaz
A~X=(A−RY*)X=AX−RY*X=AX−R=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 X≔x,Y*≔x*,M≔λ. To nám říká, že λ je vlastní číslo matice A~=A−rx*, kde r≔Ax−λx je reziduum.
Věta Nechť A∈CN×N je diagonalizovatelná jako X−1AX=diag(λ1,…,λN), μ je aproximace vlastního čísla a y je normovaná aproximace vlastního vektoru. Potom existuje λ∈σ(A) takové, že |λ−μ|≤ϰ(X)⋅r, kde r≔Ax−λx.
Důkaz Označme A~≔A−ry*. 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*‖=max‖x‖=1‖ry*x‖=max‖x‖=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 A∈ℂN×N je regulární a b∈ℂN,b≠0. 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
‖δx‖‖x‖≤ϰ(A)⋅‖δb‖‖b‖.
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‖=‖A−1δb‖≤‖A−1‖⋅‖δb‖,
‖δx‖‖x‖≤‖A−1‖‖x‖⋅‖δb‖=ϰ(A)‖A‖⋅‖x‖⋅‖δb‖<ϰ(A)⋅‖δb‖‖Ax‖=ϰ(A)⋅‖δb‖‖b‖.
Poznámka Pokud ϰ(A)≫1, potom není zaručeno, že při malém δb bude malé δx.
Lemma Platí-li
‖δA‖‖A‖<1ϰ(A),
potom je A~ regulární.
Důkaz Předpoklad je ekvivalentní tomu, že ‖δA‖⋅‖A−1‖<1. Předpokládejme pro spor, že A~ je singulární, tedy existuje x≠0 takové, že A~x=0 neboli Ax=−δAx neboli x=−A−1δAx. Z toho
‖x‖=‖A−1δAx‖≤‖A−1‖⋅‖A‖⋅‖x‖.
To je spor s předpokladem.
Věta odhad chyby pro přesnou pravou stranu Nechť b~=b. Platí-li
‖δA‖‖A‖<1ϰ(A),
potom
‖δx‖‖x‖≤ϰ(A)⋅‖δA‖‖A‖1−ϰ(A)⋅‖δA‖‖A‖.
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‖=‖A−1δAx+A−1δAδx‖≤‖A−1‖⋅‖δA‖⋅‖x‖+‖A−1‖⋅‖δA‖⋅‖δx‖,
‖δx‖⋅(1−‖A−1‖⋅‖δA‖)≤‖A−1‖⋅‖δA‖⋅‖x‖
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 ‖δA‖‖A‖≪1ϰ(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
‖δA‖‖A‖<1ϰ(A),
potom
‖δx‖‖x‖≤ϰ(A)⋅(‖δA‖‖A‖+‖δb‖‖b‖)1−ϰ(A)⋅‖δA‖‖A‖.
Důkaz Opět si rozepíšeme
(A+δA)(x+δx)=b+δb,
Aδx+δAx+δAδx=δb,
‖δx‖≤‖A−1‖⋅‖δb‖+‖A−1‖⋅‖δA‖⋅‖x‖+‖A−1‖⋅‖δA‖⋅‖δx‖.
Úpravami analogickými jako v důkazech předchozích vět dostaneme cílový odhad.
Poznámka Je-li ‖δA‖‖A‖,‖δb‖‖b‖∼10−α a ϰ(A)∼10β, potom ‖δx‖‖x‖∼10β−α.

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{‖δA‖‖A‖,‖δb‖‖b‖}=‖b−Ax~‖‖A‖⋅‖x~‖+‖b‖.
Důkaz Označme reziduum ι≔b−Ax~. Zvolme
δA≔‖A‖⋅‖x~‖‖A‖⋅‖x~‖+‖b‖⋅ι⋅x~*‖x~‖2,b~≔A~x~.
Potom
A~x=Ax~+δAx~=Ax+ι−‖b‖⋅ι‖A‖⋅‖x~‖+‖b‖=b−‖b‖⋅ι‖A‖⋅‖x~‖+‖b‖.
Z toho máme
δb=−‖b‖⋅ι‖A‖⋅‖x~‖+‖b‖.
Dokážeme, že při této konkrétní volbě δA,δb nastane rovnost ve znění věty:
‖δb‖‖b‖=‖ι‖‖A‖⋅‖x~‖+‖b‖=‖b−Ax~‖‖A‖⋅‖x~‖+‖b‖,
‖δA‖‖A‖=‖A‖⋅‖x~‖‖A‖⋅‖x~‖+‖b‖⋅‖ι⋅x~*‖‖x~2‖=‖b−Ax~‖‖A‖⋅‖x~‖+‖b‖,
kde jsme využili rovnosti
‖ιx~*‖=max‖y‖=1‖ιx~*y‖=‖ι‖⋅max‖y‖=1‖x~*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~)=‖b−Ax~‖‖A‖⋅‖x~‖+‖b‖=‖δA^x~−δb^‖‖A‖⋅‖x~‖+‖b‖≤‖δA^‖⋅‖x‖+‖δb^‖‖A‖⋅‖x~‖+‖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 A∈ℝN×N existuje Q,R∈ℝN×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
c≔cosθ=x1x12+x22,s≔sinθ=x2x12+x22.
To můžeme zobecnit na unitární matici, která pro libovolné x∈ℝN a i≠j vynuluje j-tou složku a jinak změní jen j-tou. Vezmeme
𝒢i,j≔(𝐈i−100000c0−s000𝐈j−i−1000s0c00000𝐈N−j),
kde volíme θ tak, aby s⋅xi+c⋅xj=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
Q≔∏i=1N−1∏j=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 c≔1,s≔0. Jinak:
  • pokud |xj|≥|xi|, potom
    t≔xixj,s≔11+t2,c≔s⋅t,
  • pokud |xj|<|xi|, potom
    t≔xjxi,c≔11+t2,s≔c⋅t.
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 43N2 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|≤u. Pro násobení jednou Givensovou maticí platí
fl(𝒢^i,j𝖳⋅A)=𝒢i,j(A+E),
kde ‖E‖‖A‖=K⋅u, přičemž K je řádu jednotek.
Definice Nechť u∈ℝN,‖u‖=1. Potom Hausholderova transformace je matice H≔𝐈−2uu𝖳.
Poznámka Je-li v∈ℝN,v⟂u, potom
Hu=(𝐈−2uu𝖳)u=u−2u(u𝖳u)=−u,
Hv=(𝐈−2uu𝖳)v=v−2u(u𝖳v)=v.
Doplníme-li u na ortonormální bázi ℝN, z těchto rovností plyne, že Hausholderova transformace vyjadřuje zrcadlení podle nadroviny dané normálou u.