Poznámky k předmětu Numerická lineární algebra I.

Poznámky k předmětu Numerická lineární algebra I.

Michal Merta

∗Katedra aplikované matematiky, VŠB-Technická univerzita Ostrava, e-mail:

michal.merta@vsb.cz

1 Iterační metody pro řešení soustav lineárních rovnic

Přímé metody řešení soustav lineárních rovnic (LU, LDLT, Choleského faktor- izace atd.) vyžadují O(n³)operací a velikost soustavy, kterou jimi dokážeme vyřešit, je značně omezená. V případě husté matice A ∈ R^n×n se přibližná velikost řešitelné soustavy historicky vyvíjela takto:

• 1950: n= 20(Wilkinson)

• 1965: n= 200(Forsythe a Moler)

• 1980: n= 2000(LINPACK)

• 1995: n= 20000(LAPACK)

• 2010: n= 200000(HDSS)

Pracujeme-li s řídkými maticemi, jsme schopni řešit i systémy s mnohem větší dimenzí (miliony, desítky milionů neznámých) – zejména, pokud je řešič schopen pracovat paralelně. V případě použítí přímého řešiče na řídkou matici však může dojít k jejímu zaplnění. Např. využitím prvního řádku následující matice k vynulování prvního sloupce dojde k zaplnění všech ostatních prvků v matici:













× × × × ×

× × 0 0 0

× 0 × 0 0

× 0 0 × 0

× 0 0 0 ×













→













× × × × ×

0 × × × ×

0 × × × ×

0 × × × ×

0 × × × ×













Tento problém lze částečně řešit vhodnou pivotizací.

Mezi další nevýhody přímých řešičů patří:

• Známe-li již přibližné řešení soustavy, nedokážeme tuto znalost využít ke snížení celkového počtu operací a zkrácení doby výpočtu.

• Naopak, pokud nám postačuje znalost pouze přibližného řešení, nemůžeme výpočet pomocí přímého řešiče ukončit předčasně.

Alternativou k přímým řešičům jsou iterační řešiče, které generují posloup- nost přibližných řešení {x^k} a pracují téměř výhradně s násobením matice- vektor, které má náročnostO(n²). Důležitou vlastností každé iterační metody je rychlost konvergence posloupnosti{x^k} k řešení. Může se totiž stát, že pro některé maticeAiterační metoda konverguje velmi pomalu nebo vůbec.

1.1 Lineární iterační metody

Prvním typem iteračních metod, kterým se budeme zabývat, jsou tzv. lineární iterační metody. Ty hledají posloupnost řešení soustavy Ax = b s regulární maticí ve tvaru

x^k+1:=Mx^k+Nb, (1.1)

kdeM aNjsou nějaké matice odpovídajících rozměrů¹.

Definice Lineární iterační metodu nazveme konzistentní, řeší-li rovniciMx+ Nb=xprávě jeden vektor x=A⁻¹b. Je možné ukázat, že metoda je konzis- tentní právě tehdy, je-li splněnoM=I−NA.

Definice Iterační metodu nazveme konvergentní platí-li∀b,∀x⁰ :x^k →x= A⁻¹bprok→ ∞. Je možné ukázat, že metoda je konvergentní, právě tehdy, je-li splněnokMk<1, kdekMk= max_v∈Rⁿ_,v6=o^kMvk_kvk^v

v je maticová norma indukovaná vektorovou normou.

Při odvozování následujících iteračních metod budeme využívat rozkladu maticeAna součet dolní trojúhelníkové, diagonální a horní trojúhelníhové mat- ice, tedyA=L+D+U (pozor, nepleťte si maticeL,D,U se stejně nazvanými maticemi, které se vyskytovaly u přímých řešičů). Např. matici

A=





2 1 0 1 2 1 0 1 2





rozložíme na

L=





0 0 0 1 0 0 0 1 0



,D=





2 0 0 0 2 0 0 0 2



,U=





0 1 0 0 0 1 0 0 0



.

1.1.1 Jacobiho metoda

Vyjděme z rovniceAx=b. DosazenímA=L+D+Udostaneme (L+D+U)x=b.

Roznásobme a přezávorkujme výraz na levé straně Dx+ (L+U)x=b

Jacobiho metodu odvodíme tak, že přidáme indexy k+ 1 a k k příslušným vektorůmx

Dx^k+1+ (L+U)x^k=b.

Osamostatněnímx^k+1 dostaneme předpis prok+ 1 aproximaci vektorux x^k+1:=D⁻¹(b−(L+U)x^k) =−D⁻¹(L+U)

| {z }

=M

x^k+D⁻¹

|{z}

=N

b.

1Indexkoznačuje číslo aktuální iterace.

Jednotlivé složky vektorux^k+1 můžeme vyjádřit jako

(x^k+1)i:= 1 a_i,i



bi−

n

X

j=1,j6=i

ai,j(x^k)j



= (1.2)

= 1 ai,i



bi−

i−1

X

j=1

ai,j(x^k)j−

n

X

j=i+1

ai,j(x^k)j



 (1.3)

Pro snadnější a přehlednější zápis budemei-tý prvek vektorux^k+1 také značit jako x^k+1_i (horní index tedy označuje číslo iterace, dolní index značí pořadí prvku ve vektoru).

Konzistence metody vyplývá z jejího odvození, můžeme však ještě ověřit, že M=I−NA. V případě Jacobiho metody jeM=−D⁻¹(L+U)a N=D⁻¹ (viz výše). Platí tedy

I−NA=I−D⁻¹A=I−D⁻¹(L+D+U) =

=I−D⁻¹(L+U)−D⁻¹D=−D⁻¹(L+U) =M.

Metoda je tedy konzistentní.

Lze dokázat, že metoda je konvergentní právě tehdy, kdyžkMk=kD⁻¹(L+ U)k<1. To splňují např. striktně diagonálně dominantní matice (tedy matice, pro které platí∀i:|ai,i|>Pn

j=1,j6=i|ai,j|). Konvergenci metody pro diagonálně dominantní matice se dá poměrně snadno dokázat. Vyjděme ze vztahu proi-tý prvek aproximovaného vektoru v iteracik+ 1:

x^k+1_i = 1 ai,i



b_i−

n

X

j=1,j6=i

a_i,jx^k_j



.

Jelikož je metoda konzistentní, musí tuto rovnost splňovat i prvky vektoru přes- ného řešeníx:

xi = 1 ai,i



bi−

n

X

j=1,j6=i

ai,jxj



. Odečteme-li od první rovnosti druhou, dostaneme

x^k+1_i −xi

| {z }

=e^k+1_i

=− 1 ai,i

n

X

j=1,j6=i

ai,j(x^k_j−xj

| {z }

=e^k_i

),

kde e^k je vektor chyby v k-tém kroku. Chybu v kroku k+ 1 můžeme tedy odhadnout pomocí vlastností striktně diagonálně dominantní matice:

|e^k+1_i | ≤ 1

|ai,i|

n

X

j=1,j6=i

|ai,j||e^k_j| ≤ 1

|ai,i|

n

X

j=1,j6=i

|ai,j| max

j=1,...,n,j6=i|e^k_j|=

= max

j=1,...,n,j6=i|e^k_j| 1

|ai,i|

n

X

j=1,j6=i

|ai,j|

| {z }

<1

< max

j=1,...,n,j6=i|e^k_j|

Každý prvek vektoru chyby v krokuk+ 1 je tedy v absolutní hodnotě menší než maximální prvek vektoru chyby v předchozím kroku. Vektor chyby tedy konverguje k nulovému vektoru.

1.1.2 Gaussova-Seidelova metoda

Všimněme si, že při výpočtux^k+1_i využíváme v suměPi−1

j=1a_i,jx^k_jve výrazu (1.3) pouze prvkyx^k₁, . . . , x^k_i−1. Tyto prvky tedy můžeme nahradit již vypočtenými prvky aktuálními iterace x^k+1₁ , . . . , x^k+1_i−1. Dostaneme tak předpis Gaussovy- Seidelovy metody:

x^k+1_i := 1 ai,i



bi−

i−1

X

j=1

ai,jx^k+1_j −

n

X

j=i+1

ai,jx^k_j



 (1.4)

Podobně jako v předchozím případě můžeme metodu odvodit, nahradíme-li v soustavě Ax = b matici A součtem L+D+U. Tentokrát ovšem vektor s indexemk+ 1ponecháme u součtuL+D

(L+D)x^k+1=b−Ux^k, (1.5)

tedy

x^k+1=−(L+D)⁻¹U

| {z }

=M

x^k+ (L+D)⁻¹

| {z }

=N

b.

Vztah mezi maticovým zápisem a zápisem po prvcích (1.4) je nejlépe vidět na rovnosti (1.5). Jedná se o soustavu rovnic s dolní trojúhelníkovou maticíL+D, vektorem pravé strany b−Ux^k a neznámým vektorem x^k+1. Všimněte si, že výraz (1.4) pak přesně odpovídá algoritmu pro dopřednou substituci pro řešení takovéto soustavy.

Podobně jako v případě Jacobiho metody můžeme ověřit, zda je metoda konzistentní porovnánímI−NAaM.

I−NA=I−(L+D)⁻¹A=I−(L+D)⁻¹(L+D+U) =

=I−(L+D)⁻¹(L+D)−(L+D)⁻¹U=−(L+D)⁻¹U=M.

Metoda je tedy konzistentní.

Metoda je konvergentní, právě když k(L+D)⁻¹Uk <1, což opět platí pro diagonálně dominantní matice.

1.1.3 Richardsonova metoda

Iterace Richardsonovy metody je dána předpisem x^k+1:=x^k+ωr^k,

kde ω ∈ R+ a r^k = b−Ax^k je reziduum, které určuje, jak dobře je splněna původní rovnice. Vztah mezi reziduem a chybou e^k = x^k −x lze odvodit přenásobením definice chyby maticíA

Ae^k =Ax^k−Ax=Ax^k−b=−r^k. (1.6) Studujme konvergenci metody pro symetrickou pozitivně definitní maticiA. V takovém případě vlastní číslaλ_i a vlastní vektoryv_i (tedy skaláry a vektory, pro které platíAv_i =λ_iv_i,kvk= 1) splňují

0< λ1≤λ2≤. . .≤λn

a z vlastních vektorů lze utvořit ortonormální báziRⁿ.

Díky předchozímu poznatku můžeme reziduum vyjádřit jako lineární kombi- naci prvků báze tvořené vlastními vektory, tedyr^k+1 =Pn

i=1α^k+1_i vi. Studujme nyní, jak se chová reziduum (a tedy i chyba) v jednotlivých iteracích. Na zák- ladě toho se později pokusíme odvodit optimální hodnotuωpro co nejrychlejší konvergenci.

n

X

i=1

α^k+1_i vi=r^k+1 =b−Ax^k+1=b−A(x^k

| {z }

=r^k

+ωr^k) =

=r^k−Aωr^k= (I−ωA) r^k

|{z}

=Pn i=1α^k_iv_i

= (I−ωA)

n

X

i=1

α^k_ivi=

=

n

X

i=1

α^k_iv_i−

n

X

i=1

α^k_iω Av_i

|{z}

=λ_iv_i

=

n

X

i=1

α^k_iv_i−

n

X

i=1

α_i^kωλ_iv_i=

=

n

X

i=1

(1−ωλi)α^k_ivi.

Všimněme si, že jsme vyjádřili koeficienty rozvoje reziduar^k+1 v bázi {vi}ⁿ_i=1 pomocí násobků koeficientů v předchozím kroku (viz podtržené části před- chozího výrazu). Koeficienty se tedy budou zmenšovat (a jednotlivé složky vektoru rezidua budou konvergovat k nule) právě tehdy, když |1−ωλi| < 1 pro všechnai= 1,2, . . . , n. Rychlost konvergence bude záviset na největší hod- notě|1−ωλi|. Shrněme tento poznatek do následující věty.

Věta Richardsonova metoda konverguje, právě když∀i∈ {1,2, . . . , n} :|1− ωλ_i| < 1. Konvergenční faktor ρ = maxi∈{1,2,...,n}{|1−ωλ_i|} určuje rychlost konvergence: kr^k+1k ≤ρkr^kk.

1

Obrázek 1.1: Konvergenční faktor Richardsonovy metody v závislosti naω.

Čím menší bude konvergenční faktorρ= maxi∈{1,2,...,n}{|1−ωλi|}, tím rych- leji bude metoda konvergovat. Vzhledem k tomu, že vlastní čísla maticeAjsou daná, můžeme konvergenční faktor ovlivnit pouze vhodnou volbou parametru ω. Odvození ideální hodnoty ω ilustrujme na Obrázku 1.1. Jsou na něm zná- zorněny funkcef₁(ω) =|1−ωλ₁|af_n(ω) =|1−ωλ_n|. Protože směrnice funkcí f_i(ω) =|1−ωλ_i|jsou určeny vlastními čísly maticeAa ta jsou seřazena od ne- jmenšího po největší, budou grafy všech funkcíf_i, i= 2, . . . , n−1,ležet ”mezi”

grafyf₁ af_n (v šedě vyznačené oblasti). Funkci max

i∈{1,2,...,n}{|1−ωλ_i|}

tedy můžeme vykreslit jako červeně zvýrazněnou lomenou čáru tvořenou částí funkcef1 a částí funkcefn. Z grafu této funkce tedy rovnou můžeme odvodit:

1. Interval, ve kterém musíω ležet. Aby metoda konvergovala, musí platit ρ = maxi∈{1,2,...,n}{|1−ωλ_i| < 1. Červená funkce tedy musí ležet pod zakreslenou konstantní funkcíρ= 1. Levý krajní bod intervalu je 0, pravý určíme jako průsečík příslušné části funkcef_n s konstantní funkcí 1:

−(1−ωλ_n) = 1 ⇒ ω= 2 λn

. Metoda tedy konverguje proω∈(0,2/λn).

2. Optimálníωje bod, ve kterém červeně vyznačená funkce dosahuje minima.

Tento bod dostaneme jako průsečík funkcíf₁a f_n: 1−ω_optλ1=−(1−ω_optλn) ⇒ ω_opt= 2

λ1+λn

.

Zjistili jsme tedy, že nejlepší konvergence dosáhneme, zvolíme-li ω_opt =

2

λ₁+λ_n. Konvergenční faktor bude v tomto případě ρ_opt= 1−ω_optλ₁= 1− 2λ₁

λ1+λn

= λ₁+λ_n−2λ₁ λ1+λn

=

=λn−λ1

λn+λ1 1 λ₁

1 λ1

=

λ_n λ₁ −1

λ_n

λ1 + 1 = κ(A)−1 κ(A) + 1, kdeκ(A) =λn/λ1je číslo podmíněnosti maticeA.

Můžeme také odvodit, kolik iterací je třeba, abychom dosáhli požadované relativní změny normy rezidua. Hledáme tedyk, pro které platí

kr^kk

kr⁰k ≤ε, tedy kr^kk ≤εkr⁰k Využijme toho, žekr^kk ≤ρ^k_optkr⁰ka přepišme nerovnici na

ρ^k_optkr⁰k ≤εkr⁰k.

Vykrácením normy a zlogaritmováním obou stran nerovnice dostanem řešení k ≥ _log^log_ρ^ε

opt (nezapomeňte, že protože ρ_opt ∈ (0,1), je třeba otočit znaménko nerovnosti).

1.1.4 Ukončovací podmínky

Při použití iteračního řešiče většinou nemáme předem zadaný počet iterací, které mají proběhnout. Chceme výpočet ukončit ve chvíli, kdy se s odhadem řešení dostaneme dostatečně blízko přesnému řešení. Vzhledem k tomu, že přesné řešení (tedy ani přesnou chybu v dané iteraci) neznáme, musíme si pomoci jinak.

Jednou z možností je ukončit cyklus ve chvíli, kdy se s novým odhadem řešení příliš nepohneme od předchozího odhadu (tzn. kx^k+1−x^kk< ε). Tato podmínka ale nijak nebere v potaz velikost prvků v matici soustavy a vektoru pravé strany (jiná situace nastane, pokud jsou prvky matice a vektoru v řádech tisíců, jiná pokud jsou v řádech tisícin). Proto je vhodné tuto podmínku zvolit relativně např. vzhledem k normě vektoru pravé strany (tzn. kx^k+1−x^kk <

kbkε, tedykx^k+1−x^kk/kbk< ε).

Nejčastěji se ovšem k výpočtu ukončovací podmínky používá normy vektoru reziduar^k+1=b−Ax^k+1. To nám poskytuje přirozený odhad toho, jak dobře je splněna původní rovnice. Ukončovací podmínku lze tedy volit ve tvarukb− Ax^k+1k < ε. Podobně jako v předchozím případě je i zde vhodnější použít relativní změnu rezidua oproti vektoru pravé strany (kb−Ax^k+1k/kbk < ε) nebo počátečnímu reziduu (kb−Ax^k+1k/kb−Ax⁰k< ε).

1.2 Gradientní iterační metody

Věta Řešení soustavy Ax=bse symetrickou pozitivně definitní maticíA je ekvivalentní s minimalizací kvadratické formy

f(x) = 1

2x^TAx−b^Tx.

Důkaz Dokažme nejdříve implikaciAx=b⇒x= arg min_v∈Rⁿf(v). Podíve- jme se, jak se změní funkční hodnota f, posuneme-li se z bodu x o nějaký nenulový vektorc:

f(p) =f(x+c) = 1

2(x+c)^TA(x+c)−b^T(x+c) =

=1

2x^TAx+c^T Ax

|{z}

=b

+1

2c^TAc−b^Tx−b^Tc=

=1

2x^TAx−b^Tx

| {z }

=f(x)

+c^Tb−b^Tc

| {z }

=0

+1

2c^TAc=f(x) +1 2c^TAc

| {z }

>0

.

Díky pozitivní definitnostiAje výrazc^TAckladný. Posuneme-li se tedy z bodu x v libovolném směru, hodnota funkce f se zvětší. V bodě x tedy nastává minimum.

K důkazu opačné implikacex= arg min_v∈Rⁿf(v)⇒Ax=bje třeba si uvě- domit nutnou podmínku minima funkcef :Rⁿ →R, tedy nulovost gradientu:

x= arg min

v∈Rⁿf(v)⇒ ∇f(x) =o. (1.7) Lze ukázat, že pro gradient funkcef platí

∇f(x) = ∂f

∂x1

(x), . . . , ∂f

∂xn

(x) T

= 1

2A^Tx+1

2Ax−b=Ax−b.

Z podmínky (1.7) tedy vyplýváAx−b=o.

V případě symetrické pozitivně definitní matice A∈ R^n×n máme tedy dvě možnosti, jak geometricky nahlížet na řešení soustavy lineárních rovnic. První přístup je chápat každou rovnici jako předpis nadroviny v n-rozměrném pros- toru. Řešení soustavy pak odpovídá hledání průsečíku těchto rovin. Druhý přístup, který využijeme při odvozování následujících algoritmů, odpovídá min- imalizaci příslušné pozitivně definitní kvadratické formy. Grafem pozitivně definitní kvadratické formyf : Rⁿ → R je n-dimenzionální paraboloid, který má minimum (viz Obrázek 1.2 pron= 2).

1.2.1 Metoda největšího spádu

Metoda největšího spádu je iterační metoda s předpisem

x^k+1:=x^k+α^kv^k, (1.8)

Obrázek 1.2: Graf kvadratické formy s pozitivně definitní maticí A (zdroj Wikipedia)

kde v^k volíme jako směr největšího poklesu funkce f. Všimněme si, že pro gradient platí∇f(x^k) =Ax^k−ba pro reziduum k-tém krokur^k =b−Ax^k. Tedyr^k=−∇f(x^k). Protože gradient odpovídá směru největšího růstu funkce v daném bodě, reziduum je směr největšího spádu. Logicky, protože chceme dosáhnout minima dané funkce, vydáváme se v každém kroku ve směru rezidua, tedyv^k =r^k.

Otázkou je, jak daleko se v každém kroku v tomto směru vydat, tedy jak zvolit koeficient α^k. Metoda největšího spádu volí tento koeficient tak, aby v každém kroku dosáhla minima funkce f ve směru rezidua. Definujme si tedy pomocnou funkciF :R→R:

F(α) =f(x^k+αr^k) = 1

2(x^k+αr^k)^TA(x^k+αr^k)−b^T(x^k+αr^k) =

= 1

2(x^k)^TAx^k+α(x^k)^TAr^k+1

2α²(r^k)^TAr^k−b^Tx^k−αb^Tr^k Hledámeα, ve kterém tato funkce dosahuje minima, její derivace se tedy musí rovnat nule:

F⁰(α) =α(r^k)^TAr^k+ (r^k)^T Ax^k

|{z}

=(b−r^k)

−b^Tr^k =α(r^k)^TAr^k−(r^k)^Tr= 0

Odtud

α^k= (r^k)^Tr^k

(r^k)^TAr^k. (1.9)

Stejný předpis můžeme odvodit, použijeme-li místo pomocné funkce funkcif a položíme její derivaci ve směrur^k rovnu nule (vzpomeňme si, že platí ^df(x)_dh =

(∇f(x))^Th):

df(x^k+1)

dr^k = 0 (1.10)

(∇f(x^k+1))^Tr^k = 0 (1.11)

(−r^k+1)^Tr^k = 0 (1.12)

Dosazením r^k+1 = r^k −α^kAr^k do předchozí rovnice a jednoduchou úpravou dostaneme stejný předpis proα^k jako v předchozím případě (1.16). Předchozí odvození nám také prozradilo důležitou vlastnost metody největšího spádu – každý směr v^k = r^k je kolmý na předchozí směr. Jak brzy uvidíme, není to vždy žádaná vlastnost.

Algoritmus tedy počítá jednotlivé aproximace pomocí následujících předpisů:

r^k :=b−Ax^k =b−A(x^k−1+α^k−1r^k−1) =b−Ax^k−1

| {z }

r^k−1

−α^k−1Ar^k−1=

=r^k−1−α^k−1Ar^k−1 α^k := (r^k)^Tr^k

(r^k)^TAr^k x^k+1 :=x^k+α^kr^k

Díky úpravě předpisu pro výpočetr^kjsme ušetřili jedno násobení matice-vektor (Ax^k) – výsledekAr^k−1 si totiž můžeme zapamatovat z předchozí iterace.

Ukažme si nyní, že metoda konverguje. Konvergenci budeme dokazovat v tzv. energetické normě k · kA, tedy normě indukované skalárním součinem (Ax,x) =x^TAx:kxkA=√

x^TAx.

ke^k+1k²_A= (e^k+1)^TAe^k+1 = (e^k+α^kr^k)^TA(e^k+α^kr^k) =

= (e^k)^TAe^k+ 2α^k(r^k)^T Ae^k

|{z}

=−r^k

+(α^k)²(r^k)^TAr^k =

=ke^kk²_A−2 (r^k)^Tr^k (r^k)^TAr^k

| {z }

=α^k

(r^k)^Tr^k+ ( (r^k)^Tr^k

(r^k)^TAr^k)²(r^k)Ar^k =

=ke^kk²_A−((r^k)^Tr^k)² (r^k)^TAr^k =

=ke^kk²_A

1− ((r^k)^Tr^k)² ((r^k)^TAr^k) ((e^k)^TAe^k)

| {z }

=(A⁻¹r^k)^TA(A⁻¹r^k)=(r^k)^TA⁻¹r^k

.

Rovnost v poslední závorce vyplývá ze vztahu (1.6) mezi reziduem a chybou.

Zjistili jsme tedy, že chyba v kroku k+ 1 je nějakým násobkem chyby v před-

chozím kroku. Podívejme se na tento násobek podrobněji:

1− ((r^k)^Tr^k)²

((r^k)^TAr^k)((r^k)^TA⁻¹r^k) ≤1− ((r^k)^Tr^k)²

kAkkr^kk²kA⁻¹kkr^kk² =

= 1− 1

kAkkA⁻¹k =q,

kde první nerovnost vyplývá z Cauchyho-Schwarzovy nerovnosti a definice mati- cové normy indukované vektorovou normou. VýrazkAkkA⁻¹knazývýme číslem podmíněnosti (viz také Richardsonova metoda) a v případě symetrické pozitivně definitní matice jej vypočítáme jako podíl největšího a nejmenšího vlastního čísla κ(A) = ^λ_λ^max

min.

Protožeκ(A) =λmax/λmin≥1, platíq= 1−1/κ(A)<1. Chyba se tedy v každém kroku zmenšuje a metoda konverguje k řešení. Navíc

ke^k+1k²_A≤

1− 1 κ(A)

ke^kk²_A

a je možné dokázat vztah mezi počáteční chybou a chybou vk-té iteraci:

ke^k+1kA≤

κ(A)−1 κ(A) + 1

k

ke⁰kA. (1.13)

Metoda tedy konverguje k řešení pomalu. Zejména je-li κ(A) velké, je zlomek ve výrazu (1.13) blízký 1 a dostaneme velmi pomalou konvergenci.

Celý algoritmus můžeme napsat na několika řádcích:

functionsteepestd_descent(A,b,x⁰) r⁰=b−Ax⁰

k= 0

whilekr^kk/kr⁰k> εdo αk= ((r^k)^Tr^k)/((r^k)^TAr^k) x^k+1=x^k+αkr^k

r^k+1=r^k−αkAr^k k=k+ 1

end while end function

Poznámka Z (1.12) vyplývá, že směrr^k+1, ve kterém se v každém kroku vy- dáme, je kolmý na předchozí směr. To může vést k situacím jako na Obrázku 1.3, kdy procházíme „cik-cak“ údolí tvořené minimalizovanou kvadratickou formou.

Vydáváme se tedy několikrát v tom samém směru a postupně zkracujeme krok. Ideální by bylo najít metodu, která by minimalizovala celkovou chybu v daném směru vždy pouze jednou (v každé iteraci vynulovala chybu v daném směru). Taková metoda by byla schopná vyřešit soustavu onneznámých během n iterací. Zkusme si tedy pro A ∈ R^n×n zvolit n navzájem ortogonálních směrůd⁰,d¹, . . . ,dⁿ⁻¹ (např. rovnoběžných s osami souřadnic) a použít před- pis: x^k+1=x^k+αkd^k.

x

x⁰

Obrázek 1.3: Pomalá konvergence metody největšího spádu (vrstevnice mini- malizované kvadratické formy a jednotlivé směryr^k).

Jak ale získatαk? Ilustrujme to ve dvou dimenzích na Obrázku 1.4. Na za- čátku (Obrázek 1.4 vlevo) si zvolíme nějaký počáteční odhad řešeníx⁰; počáteční chybu e⁰ pak lze rozložit na dvě kolmé složky rovnoběžné s osami souřadnic (e⁰ = e⁰_x+e⁰_y). V první iteraci si zvolme vektor d⁰ rovnoběžný s osou x (Obrázek 1.4 vpravo). Z obrázku je patrné, že abychom vynulovali chybu ve směru osyx, musíme bodx¹ zvolit tak, aby vektor nové chybye¹ byl kolmý k d⁰.

Obecně tedy chceme, aby nová chyba byla kolmá k aktuálnímu směru:

(d^k)^Te^k+1= 0.

Z definice chyby dostaneme

0 = (d^k)^T(x^k+1−x) = (d^k)^T(x^k+αkd^k−x) = (d^k)^T(e^k+αkd^k) a odtud

αk=−(d^k)^Te^k (d^k)^Td^k.

Je jasné, že jsme si příliš nepomohli. Abychom vypočítaliαk, potřebovali by- chom znát chybue^k. Tu ale bez znalosti přesného řešení nespočítáme. Jak si ale ukážeme v příští lekci, pokud požadavek na ortogonalitu směrůd^knahradíme za tzv. A-ortogonalitu, budeme skutečně schopni odvodit metodu, která (alespoň teoreticky) dokonverguje k řešení běhemniterací.

1.2.2 Metoda sdružených gradientů

Definice Buď A symetrická pozitivně definitní matice. Nenulové vektory {pi}ⁿ⁻¹_i=0 nazveme sdružené (A-ortogonální), platí-li

∀i, j∈ {0,2, . . . , n−1}:i6=j ⇒p^T_iApj= 0.

Dva navzájemA-ortogonální vektory někdy označujemea⊥A b.

x

x0

e⁰

x y

e_x⁰

e_y⁰

x

x0

x y

e¹ x1

d⁰

Obrázek 1.4: Ortogonální složky počáteční chyby (vlevo). Chceme-li vynulovat složku chyby ve směrud⁰, musí býte¹ kolmý kd⁰(vpravo).

Lemma Prvky množiny sdružených vektorů jsou lineárně nezávislé.

Důkaz Chceme ukázat, že platí

α0p0+α1p1+. . .+α_n−1p_n−1= 0⇒ ∀i:αi= 0.

Přenásobme rovnici výrazem(Apk)^T, k∈ {0,1, . . . , n−1} zleva:

α0p^T_kAp0+α1p^T_kAp1+. . .+αn−1p^T_kApn−1= 0.

Protože jsou vektory sdružené, všechny členy v této sumě až na jeden jsou nulové:

α_kp^T_kAp_k = 0.

Současně víme, že matice A je symetrická pozitivně definitní a součin p^T_kApk

musí být kladný. Takžeαk = 0.

Vraťme se k poznámce na konci předchozí kapitoly a ortogonální směry nahraďme zaA-ortogonální. Předpokládejme tedy, že známenA-ortogonálních směrů{d⁰,d¹. . . ,dⁿ⁻¹}a podobně jako na konci předchozí kapitoly se pokusme odvodit iterační metodu ve tvarux^k+1=x^k+α_kd^k. Požadavek na ortogonalitu nové chyby k aktuálnímu směru nahraďme za požadavek na jejich vzájemnou A-ortogonalitu, tzn. e^k+1 ⊥_A d^k. Ukažme si, že toto nastane v minimu naší kvadratické formyf(x) = 1/2x^TAx−b^Txve směru d^k. Hledejme tedy nový bodx^k+1 tak, aby se v něm derivacef ve směrud^k rovnala nule (nutná pod- mínka minima):

df(x^k+1)

dd^k = gradf(x^k+1)^Td^k =−( r^k+1

| {z }

=−Ae^k+1

)^Td^k= (e^k+1)^TAd^k = 0.

Chybae^k+1 je tedy v bodu minima skutečněA-ortogonální ke směrud^k. Nyní snadno můžeme získat předpis pro koeficientα_k(použijeme již několikrát zmíněný vztahe^k+1=e^k+α_kd^k):

(e^k+1)^TAd^k= (e^k+α_kd^k)^TAd^k= 0⇒α_k =−(e^k)^TAd^k

(d^k)^TAd = (r^k)^Td^k (d^k)^TAd

Věta Algoritmus spočte přesné řešeníxv nejvýšenkrocích.

Důkaz Vyjádříme si počáteční chybue⁰v bázi tvořenéA-ortogonálními (tedy nezávislými) vektory{dⁱ}ⁿ⁻¹_i=0:

e⁰=

n−1

X

j=0

δjd^j.

Souřadniceδj najdeme tak, že rovnici přenásobíme zleva výrazem(d^k)^TA, k∈ {0,1, . . . , n−1} a využijeme toho, že většina členů výsledné sumy bude díky A-ortogonalitě nulová.

(d^k)^TAe¹=

n−1

X

j=0

δj(d^k)^TAd^j=δk(d^k)^TAd^k. Odtud

δk = (d^k)^TAe¹

(d^k)^TAd^k = (d^k)^TA(e¹+Pk−1 i=0 α_idⁱ)

(d^k)^TAd^k = (d^k)^TAe^k

(d^k)^TAd^k =− (d^k)^Tr^k (d^k)^TAd^k. Zde si uvědomme, že přičtením sumy Pk−1

i=0 αidⁱ jsme výraz nijak nezměnili (všechna dⁱ v sumě jsou A-ortogonální k d^k před závorkou, součin je tedy nulový). Dále rovnost e⁰ +Pk−1

i=0 αidⁱ = e^k snadno vyplyne z již známého vztahue^k+1=e^k+αkd^k.

Všimněme si, žeδk =−αk. Vyjádřeme si tedy chybu vi-té iteraci jako

eⁱ =e⁰+

i−1

X

j=0

αjd^j=

n−1

X

j=0

δjd^j+

i−1

X

j=0

αjd^j=

n−1

X

j=0

δjd^j−

i−1

X

j=0

δjd^j=

n−1

X

j=i

δjd^j.

V každé iteraci se tedy odstraní jedna složka počáteční chybye⁰. Poniteracích je každá složka chyby v bázi{di}ⁿ⁻¹_i=0 vynulovaná⇒eⁿ =o.

Zbývá ukázat, jak najít jednotlivé A-ortogonální směry {dⁱ}ⁿ⁻¹_i=0. Před- pokládejme, že známe {d⁰,d¹, . . . ,d^k} a hledejme d^k+1. Vyjdeme z vektoru rezidua r^k+1 a pomocí Gramova-Schmidtova procesu (připomeňte si) jej A- ortogonalizujeme vůči předchozím směrům. Hledejme nový vektor ve tvaru

d^k+1=r^k+1−

k

X

j=0

β_k,jd^j. (1.14)

Určeme koeficientyβk,jtak, aby výsledný vektor bylA-ortogonální k předchozím vektorům{d⁰,d¹, . . . ,d^k}. Přenásobme proto rovnost zleva výrazem(dⁱ)^TA, i∈ {0,1, . . . , k}

(dⁱ)^TAd^k+1= (dⁱ)^TAr^k+1−

k

X

j=0

βk,j(dⁱ)^TAd^j

a využijme toho, že díkyA-ortogonalitě je výraz na levé straně rovnosti stejně jako většina členů sumy napravo roven nule. Tedy

0 = (dⁱ)^TAr^k+1−β_k,i(dⁱ)^TAdⁱ ⇒β_k,i= (dⁱ)^TAr^k+1 (dⁱ)^TAdⁱ . Dosazením vypočtených koeficientů do (1.14) dostaneme

d^k+1=r^k+1−

k

X

j=0

(d^j)^TAr^k+1

(d^j)^TAd^j d^j. (1.15) Všimněme si, že abychom pomocí tohoto vzorce vypočítali nový směr, museli bychom si pamatovat všechny předchozí směry, což by z paměťového hlediska bylo značně nevýhodné. Ukážeme si ale nyní, že všechny koeficientyβk,j kromě βk,k jsou nulové. Dokažme tedy následující tři lemmata.

Lemma Reziduumr^k+1je ortogonální kd⁰,d¹ . . . ,d^k. Důkaz Přenásobme r^k+1vektorydⁱ, i∈ {0,1, . . . , k}.

(dⁱ)^T r^k+1

| {z }

=−Ae^k+1

=−(dⁱ)^TAe^k+1=−(dⁱ)^TA

n−1

X

j=k+1

δjd^j= 0,

kde poslední rovnost vyplývá zA-ortogonality všech vektorůd^jv sumě k vektoru dⁱ.

Lemma Reziduumr^k+1je ortogonální kr⁰,r¹, . . . ,r^k. Důkaz Využijme předpis (1.14) a vyjádřeme rⁱ:

rⁱ=dⁱ+

i−1

X

j=0

βi,jd^j.

Pak díky předchozímu lemmatu dostaneme pro všechnai≤k:

(r^k+1)^Trⁱ= (r^k+1)^T(dⁱ+

i−1

X

j=0

βi,jd^j) = (r^k+1)^Tdⁱ+

i−1

X

j=0

βi,j(r^k+1)^Td^j= 0.

Lemma Nechť je dána množina reziduí{rⁱ}ⁿ⁻¹_i=0 a sdružených směrů{dⁱ}ⁿ⁻¹_i=0. Potom platí:

j < k: (r^k+1)^TAd^j= 0 j=k: (r^k+1)^TAd^j6= 0

Důkaz Pro reziduum platí

r^j+1=b−Ax^j+1=b−A(x^j+αjd^j) =r^j−αjAd^j Odtud

Ad^j =− 1

α_j(r^j+1−r^j) a

(r^k+1)^TAd^j = (r^k+1)^T(1

α_j(r^j−r^j+1)) = 1

α_j(r^k+1)^Tr^j− 1

α_j(r^k+1)^Tr^j+1 Proj < kje tento výraz díky ortogonalitě reziduí nulový. Proj =kdostaneme

(r^k+1)^TAd^k= 1 αk

(r^k+1)^Tr^k− 1 αk

(r^k+1)^Tr^k+1=− 1 αk

(r^k+1)^Tr^k+1. (1.16) Vraťme se nyní zpět k předpisu pro výpočet koeficientů βk,i v Gramově- Schmidtově procesu. Díky předchozímu lemmatu jsme zjistili, že pro j = 0, . . . , k−1 jsou čitatelé v sumě v (1.15) nulové. Předpis se tedy zjednoduší na

d^k+1=r^k+1−(d^k)^TAr^k+1 (d^k)^TAd^k d^k.

K výpočtu nového směru nám stačí znát nové reziduum a předchozí směr. Pro- tože všechny ostatní koeficienty v dané sumě jsou nulové, budeme pro jednodu- chost psát místoβk,kjenβk.

Využijeme-li rovnosti r^k+1 = r^k −α_kAd^k, tedy Ad^k = _α¹

k(r^k −r^k+1), můžeme přepsat čitatele na

(r^k+1)^TAd^k = 1 αk

(r^k+1)^T(r^k−r^k+1) =− 1 αk

(r^k+1)^Tr^k+1. Podobně lze upravit jmenovatele

(d^k)^TAd^k= (r^k−β_k−1d^k−1)^TAd^k= 1

α_k(r^k)^T(r^k−r^k+1) = 1

α_k(r^k)^Tr^k. Po dosazení do předpisu proβ_k dostaneme

βk =−(r^k+1)^Tr^k+1 (r^k)^Tr^k

V podobném duchu můžeme upravit i čitatele v předpisuα_k

(r^k)^Td^k = (r^k)^Td^k =r^k−β_k−1d^k−1= (r^k)^Tr^k−β_k,k−1(r^k)^Td^k−1

| {z }

=0

,

čímž získáme

α_k = (d^k)^Tr^k

(d^k)^TAd^k = (r^k)^Tr^k (d^k)^TAd^k.

Pokud použijme tyto předpisy pro αk, βk, ušetříme násobení matice-vektor a skalární součin.

Zapišme nyní celý algoritmus:

functionconjuate_gradient(A,b, x⁰) d⁰=r⁰=b−Ax⁰

k= 0

whilekr^kk/kr⁰k> εdo α_k= _(d^(r_k^k₎⁾_T^T_Ad^rk_k x^k+1=x^k+αkd^k r^k+1=r^k−α_kAr^k βk= ^(rk+1_(rk⁾₎^T_T^r_r^k+1_k d^k+1=r^k+1+β_kd^k k=k+ 1

end while end function

Doplňme ještě informace o konvergenci. Chybu můžeme odhadnout pomocí vztahů

ke^k+1kA≤ κ(A)−1 κ(A) + 1ke^kkA

a

ke^kkA≤2

pκ(A)−1 pκ(A) + 1

!k

ke⁰kA.

Můžeme také odhadnout maximální počet iterací nutných k dosažení relativní přesnostiε(tedyke^kk ≤εke⁰k):

k≤ d1 2

pκ(A) ln2 εe.

Porovnejme to s metodou největšího spádu, pro kterou platí k≤ d1

2κ(A) ln1 εe.

Rozdíl mezi metodami je dobře patrný na Obrázku 1.5.

Na závěr poznamenejme, že metoda sdružených gradientů patří mezi tzv.

Krylovovské metody. Tyto metody generují posloupnost Krylovových podpros- torů ve tvaruK^k(A,b) = span

b,Ab,A²b, . . . ,A^kb a v každé iteraci minimal- izují chybu v energetické normě na daném podprostoru, tedy

x^k+1= arg min

x∈x⁰+K^k(A,r⁰)

kx−xkA,

kdexje přesné řešení.

10^-15 10^-10

10^-5

Relativni presnost 0

20 40 60 80 100 120 140 160 180 200

Max. pocet iteraci

Sdruzene gradienty Nejvetsi spad

Obrázek 1.5: Porovnání maximálního počtu iterací nutných k dosažení dané relativní přesnosti (κ(A) = 10).

1.3 Předpodmínění

V předchozích kapitolách jsme viděli, že rychlost konvergence iteračních metod k řešení závisí na čísle podmíněnosti κ(A) matice soustavy. Předpodmínění je nenáročná ekvivalentní úprava soustavy, jejímž cílem je snížení čísla pod- míněnosti, tedy snížení počtu potřebných iterací.

Předpokládejme, že známe matici M takovou, že platí κ(M⁻¹A) κ(A).

Pak místo původní soustavyAx=bmůžeme řešit ekvivalentní systém

M⁻¹Ax=M⁻¹b. (1.17)

Vzhledem k menšímu číslu podmíněnosti matice M⁻¹Aby měl oproti původní soustavě klesnout počet iterací nutných k jejímu vyřešení. MaticiMnazýváme levým předpodmiňovačem.

Je-li matice soustavyA symetrická pozitivně definitní, je výhodné původní soustavu řešit metodou největšího spádu nebo metodou sdružených gradientů.

Po předpodmiňovači M pak budeme chtít, aby byl také symetrický pozitivně definitní. Problémem však je, že součin dvou symetrických pozitivně definitních matic obecně nemusí být symetrická pozitivně definitní matice. Systém (1.17) tedy můžeme vyřešit pomocí iteračního řešiče, který nevyžaduje tyto vlastnosti matice, ale na rozdíl od původní soustavy Ax = b nemůžeme použít metodu největšího spádu ani metodu sdružených gradientů.

Pokusme se tedy odvodit takový tvar předpodmínění, abychom zachovali symetrii a pozitivní definitnost matice systému. Využijeme toho, že každou sy- metrickou pozitivně definitní matici můžeme rozložit (např. pomocí Choleského rozkladu) na součinM=LL^T, kdeLje dolní trojúhelníková matice. Dále využi- jme toho, že maticeM⁻¹A a L⁻¹AL^−T mají stejná vlastní čísla (a tedy stejné číslo podmíněnosti). Je-li totižv vlastní vektor maticeM⁻¹As vlastním číslem λ, pakL^Tv je vlastní vektorL⁻¹AL^−T se stejným vlastním číslemλ:

(L⁻¹AL^−T)(L^Tv) =L⁻¹A L^−TL^T

| {z }

=I

v=L⁻¹Av= (L^TL^−T

| {z }

=I

)L⁻¹Av

=L^TL^−TL⁻¹

| {z }

=M⁻¹

Av=L^TM⁻¹Av

| {z }

=λv

=λL^Tv Přenásobme nyní původní soustavuAx=bzleva maticíL⁻¹:

L⁻¹Ax=L⁻¹b

Abychom zachovali symetrii a pozitivní definitnost, musíme maticiApřenásobit L^−T zprava. “Vložme” tedy meziAaxjednotkovou matici ve tvaruI=L^−TL^T, tím soustavu nijak nezměníme:

L⁻¹AL^−TL^Tx=L⁻¹b.

Místo původní soustavy Ax = b a místo soustavy (1.17) řešme ekvivalentní systém

(L⁻¹AL^−T)(L^Tx) =L⁻¹b.

Zaveďme značení Aˆ = L⁻¹AL^−T,xˆ = L^Tx,ˆb = L⁻¹b. Soustavu pak můžeme zapsat jako

Aˆˆx= ˆb. (1.18)

V tomto případě je matice soustavy Aˆ = L⁻¹AL^−T symetrická a pozitivně definitní, k řešení (1.18) tedy můžeme použít metodu největšího spádu či metodu sdružených gradientů. Po nalezení xˆ získáme řešení původní soustavy jako x=L^−Tx.ˆ²

Poznamenejme, že nevýhodou tohoto přístupu je nutnost znalosti rozkladu M=LL^T. Vhodnou manipulací s předpisy jednotlivých iteračních metod se však této potřeby zbavíme. Demonstrujme to pro jednoduchost na Richardsonově metodě, která má při aplikaci na (1.18) předpis

ˆ

x^k+1= ˆx^k+ω(ˆb−Aˆˆx^k),

2Inverzní maticeL⁻¹aL^−T samozřejmě nemusíme explicitně sestavovat. Např. násobení w= L⁻¹uje ekvivalentní s řešením soustavy Lw = u. Matice Lje dolní trojúhelníková, řešení tedy snadno získáme pomocí dopředné substituce. Obdobně postupujeme i v případě násobení maticíL^−T. Není také třeba a není výhodné, abychom explicitně sestavovali matici L⁻¹AL^−Tjakou maticový součin. SoučinL⁻¹AL^−Tuvypočítáme jako posloupnost tří po sobě jdoucích součinů matice-vektor.

tedy

L^Tx^k+1=L^Tx^k+ω(L⁻¹b−L⁻¹AL^−TL^Tx^k).

Přenásobme tento předpis zleva maticíL^−T:

L^−TL^Tx^k+1=L^−TL^Tx^k+ω(L^−TL⁻¹b−L^−TL⁻¹AL^−TL^Tx^k).

Odtud snadno získáme předpis předpodmíněné Richardsonovy metody, který již nevyžaduje znalost rozkladuM:

x^k+1=x^k+ωM⁻¹(b−Ax^k).

Podobným způsobem bychom mohli upravit předpisy metody největšího spádu či metody sdružených gradientů. Všimněme si, že při použití tohoto předpisu potřebujeme aplikovatM⁻¹(tedy řešit soustavu s maticí M). Předpodmiňovač tedy musíme volit tak, abychom toto řešení byli schopni najít rychle a efektivně.

Jak tedy zvolit M? Požadavky na předpodmiňovač se dají shrnout takto:

inverze M by měla aproximovat inverzi A (M⁻¹ ≈ A⁻¹), řešení soustavy s M by mělo být snadné a efektivní a pokud je matice A symetrická a pozitivně definitní, měl by i předpodmiňovač být symetrický pozitivně definitní. Dva extrémní případy jsou:

• M =A – v takovém případě inverze M aproximuje inverzi A přesně, ale aplikace samotného předpodmínění je stejně náročná jako řešení původní úlohy;

• M = I – v tomto případě je sice aplikace předpodmiňovače triviální, nedosáhneme ale žádného zlepšení čísla podmíněnosti (získáme původní nepředpodmíněný algoritmus).

Mezi těmito nepoužitelnými extrémy existuje celé spektrum sofistikovaných metod předpodmínění, pro jednoduchost však zmiňme dva dobře známé přís- tupy:

• diagonální předpodmiňovač – volíme M = diagA. V tomto případě je výpočet inverzeMvelmi snadný. Tento předpodmiňovač je efektivní např.

pro soustavy s diagonálně dominantní maticí.

• neúplný LU/Choleského rozklad – je vhodný pro soustavy s řídkou maticí.

Problémem klasického LU/Choleského rozkladu je, že přestože je matice A řídká, mohou být její faktory L,U plné nebo více zaplněné. Princip neúplného LU/Choleského rozkladu je jednoduchý - použijeme stejný al- goritmus jako u klasického LU/Choleského rozkladu, ale nenulový prvek do matice L nebo U uložíme pouze tehdy, existuje-li nenulový prvek na stejné pozici v rozkládané matici A. Dostaneme neúplné faktory L,˜ U˜ a příslušné předpodmíňovače definujeme jakoM= ˜LU˜ neboM= ˜L˜L^T. Řešit soustavu sMje pak snadné, protože známe její rozklad na součin trojúhel- níkových matic a můžeme použít dopřednou a zpětnou substituci.

References

[1] Trefethen, L. N, Bau, D. Numerical Linear Algebra. SIAM. 1997.

[2] Schewchuk, J. R. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. 1994. Dostupné zhttps://www.cs.cmu.edu/~quake-papers/painless- conjugate-gradient.pdf

[3] Lukáš, D. Zápisky z přednášek. Dostupné zhttps://homel.vsb.cz/~luk76