Zagadnienia

6.1 Metoda gradientów sprzężonych (CG)
- 6.1.1 Implementacja
6.2 Metoda GMRES
- 6.2.1 Implementacja
- 6.2.2 Wersja z restartem: GMRES( $m$ )

6. Metody iteracyjne Kryłowa

W dalszym ciągu przedmiotem naszego zainteresowania będzie układ równań $Ax=b$ , którego rozwiązanie oznaczymy (jak zwykle) przez $x^{*}$ . O ile nie będzie jasno zaznaczone, że jest inaczej, będziemy przyjmować że początkowe przybliżenie rozwiązania, $x_{0}$ , jest dowolnie zadanym wektorem w $R^{N}$ (na przykład, $x_{0}=0$ ).

Tym razem rozważymy pewien wariant metody projekcji, w którym rozmiar podprzestrzeni będzie powiększał się wraz z postępem iteracji. Kolejne przybliżenie $x_{k}$ będziemy dobierać w taki sposób, by $x_{k}\in x_{0}+K_{k}(r_{0},A)$ oraz spełniało pewien dodatkowy warunek, na przykład — minimalizowało pewną miarę błędu na przestrzeni Kryłowa

$K_{k}(r_{0},A)=\mbox{span}\{ r_{0},Ar_{0},\ldots,A^{{k-1}}r_{0}\},$

gdzie $r_{0}=b-Ax_{0}$ jest residuum początkowym. (Zwróć uwagę, że przestrzeń Kryłowa jest rozpięta przez kolejne wektory metody potęgowej, o której więcej w rozdziale LABEL:sec:xxx — to nie przypadek!). Tam, gdzie nie będzie to prowadziło do nieporozumień, będziemy pomijali parametry przestrzeni Kryłowa i pisali po prostu $K_{k}$ zamiast $K_{k}(r_{0},A)$ .

Metody zbudowane zgodnie z powyższym schematem będziemy ogólnie nazywać metodami Kryłowa. Jak za chwilę zobaczymy, taki sposób konstrukcji iteracji pozwoli metodzie samodzielnie ,,dostosować się” do własności macierzy, przez co metody Kryłowa, w przeciwieństwie do metod typu relaksacji, będą w stanie wykorzystać korzystne własności macierzy do przyspieszenia zbieżności: coś, czego ani metody stacjonarne, ani proste metody projekcji nie były w stanie osiągnąć!

Jasne jest, że przestrzenie Kryłowa tworzą wstępujący ciąg podprzestrzeni:

$K_{0}\subseteq K_{1}\subseteq\ldots\subseteq K_{k}\subseteq K_{{k+1}}\subseteq\ldots\subseteq\mathbb{R}^{N}.$

Łatwo pokazać, że nie tylko same przybliżenia $x_{k}$ , ale także błędy i residua na $k$ -tym kroku metody Kryłowa można powiązać z pewnymi przestrzeniami Kryłowa:

Stwierdzenie 6.1

Jeśli $x\in x_{0}+K_{k}(r_{0},A)$ , to $x-x^{*}\in K_{{k+1}}(x_{0}-x^{*},A)$ oraz $b-Ax\in K_{{k+1}}(r_{0},A)$ .

Z tego faktu wynika, że metody przestrzeni Kryłowa są metodami wielomianowymi: zarówno błąd, jak i residuum dadzą wyrazić się jako pewien wielomian macierzy $A$ działający na, odpowiednio, błędzie albo residuum początkowym:

Wniosek 6.1

Jeśli $x\in x_{0}+K_{k}(r_{0},A)$ , to błąd $x-x^{*}$ jest postaci $(I+c_{1}A+\ldots+c_{k}A^{k})(x_{0}-x^{*})$ i analogicznie residuum $b-Ax$ jest postaci $(I-c_{1}A-\ldots-c_{k}A^{k})r_{0}$ , gdzie $c_{1},\ldots,c_{k}$ są pewnymi rzeczywistymi współczynnikami.

Ćwiczenie 6.1

Udowodnij powyższe stwierdzenie i wniosek.

Wskazówka:

Wystarczy zauważyć, że $x\in x_{0}+K_{k}(r_{0},A)$ jest postaci $x=x_{0}+\sum _{{i=0}}^{{k-1}}c_{{i+1}}A^{i}r_{0}$ .

Konstrukcja dobrej metody Kryłowa powinna więc (jeśli mielibyśmy zachować ducha metod projekcji z poprzedniego rozdziału) zmierzać w stronę wskazania takiego wielomianu, by norma błędu była jak najmniejsza — na każdym kroku iteracji będziemy wtedy musieli rozwiązać pewne zadanie optymalizacyjne. W zależności od wyboru sposobu miary błędu i warunku optymalności, dostaniemy inną metodę iteracyjną: CG, GMRES, PCR, BiCG i inne. W niniejszym wykładzie omówimy pokrótce tylko najpopularniejsze: CG dla macierzy symetrycznych i dodatnio określonych oraz GMRES dla pozostałych. Należy pamiętać, że w zależności od konkretnego zadania, inne metody mogą okazać się bardziej skuteczne od tutaj omawianych.

W praktycznej implementacji, metody przestrzeni Kryłowa są metodami, w których do ich realizacji musimy jedynie umieć wykonać mnożenie macierzy przez wektor — nie jest potrzebne odwoływanie się do poszczególnych elementów lub części macierzy (co było konieczne w metodach opartych na podziale macierzy, takich jak np. metoda Gaussa–Seidela). Metody, które wymagają wykonywania jedynie mnożenia macierzy przez wektor, a nie odwołują się do poszczególnych elementów macierzy, nazywamy metodami operatorowymi (ang. matrix-free methods). Są one szczególnie pożądane w przypadku, gdy macierz ma nieregularną strukturę i jest bardzo wielkiego rozmiaru, a także w przypadku, gdy sama macierz po prostu nie jest jawnie wyznaczona.

Innym ważnym zastosowaniem metod iteracyjnych Kryłowa jest szybka realizacja jednego kroku tzw. niedokładnej metody Newtona rozwiązywania układu równań nieliniowych, por. rozdział 10.2.4.

Ćwiczenie 6.2

Czy metoda

Jacobiego
Richardsona
najszybszego spadku

jest metodą operatorową?

Na zakończenie wskażemy ciekawą teoretyczną właściwość metod Kryłowa opartych na minimalizacji. Przez $\| x\| _{B}$ będziemy oznaczali normę energetyczną wektora $x$ , indukowaną przez symetryczną, dodatnio określoną macierz $B$ ,

$\| x\| _{B}=\sqrt{x^{T}Bx}.$

Twierdzenie 6.1 (metody Kryłowa oparte na minimalizacji błędu albo residuum w normie energetycznej)

Niech $B\in\mathbb{R}^{N}$ będzie pewną macierzą symetryczną i dodatnio określoną i niech $k$ -ta iteracja metody Kryłowa, $x_{k}\in x_{0}+K_{k}$ , będzie określona przez jeden z warunków:

minimalizacji błędu:

$\| x_{k}-x^{*}\| _{B}\leq\| x-x^{*}\| _{B}\qquad\forall x\in x_{0}+K_{k},$

albo
minimalizacji residuum:

$\| r_{k}\| _{B}\leq\| b-Ax\| _{B}\qquad\forall x\in x_{0}+K_{k}.$

Wtedy metoda jest dobrze określona i znajduje dokładne rozwiązanie $x^{*}$ najpóźniej po $N$ iteracjach (jest to tzw. własność skończonego stopu).

Powyższe twierdzenie stosuje się do wielu konkretnych metod Kryłowa znajdujących się w powszechnym użyciu. Przykładowo, w metodzie CG (którą dokładniej omówimy w rozdziale 6.1), zdefiniowanej dla macierzy $A$ symetrycznej i dodatnio określonej, iterację będziemy określać warunkiem minimalizacji błędu przy $B=A$ , natomiast w metodzie GMRES (opisanej w rozdziale 6.2) wybierzemy warunek minimalizacji residuum dla $B=I$ .

Część druga powyższego twierdzenia — o skończonym stopie metody — ma znaczenie głównie teoretyczne. Wynika to z dwóch przesłanek:

gdy $N$ jest bardzo duże, wykonanie $N$ iteracji, nawet jeśli każda kosztuje tylko $O(N)$ flopów, jest zwykle ponad nasze możliwości obliczeniowe;
w implementacji metody stosuje się rozmaite chwyty służące maksymalnemu obniżeniu czasochłonności i pamięciożerności algorytmu, które zwykle opierają się na subtelnych zależnościach pomiędzy wyznaczanymi wartościami pośrednimi; gdy algorytm realizujemy w arytmetyce skończonej precyzji, z powodu błędów zaokrągleń tracimy te relacje i z postępem iteracji będą one zwykle coraz gorzej spełnione — co w efekcie doprowadzi także do utraty własności skończonego stopu algorytmu.

Dowód

(twierdzenia 6.1) Rozważmy najpierw przypadek 6.1, tzn. minimalizacji normy błędu. Pokażemy, że w ogólności $x_{k}$ jest dobrze określony przez warunek minimalizacji, gdyż jest zdefiniowany jako rozwiązanie pewnego liniowego zadania najmniejszych kwadratów.

Jeśli przez $V_{k}$ oznaczymy macierz, której kolumnami są wektory tworzące bazę $K_{k}$ , to każdy $x\in x_{0}+K_{k}$ daje się zapisać w postaci $x=x_{0}+V_{k}a$ , gdzie $a$ jest wektorem o tylu współrzędnych, jaki jest wymiar $K_{k}$ (oznaczmy go $d_{k}=\dim{K_{k}}$ ). Podstawiając do minimalizowanego wyrażenia, dostajemy zadanie znalezienia $a_{k}\in\mathbb{R}^{{d_{k}}}$ takiego, że

$\displaystyle\| x_{0}+V_{k}a_{k}-x^{*}\| _{B}\leq\| x_{0}+V_{k}a-x^{*}\| _{B}\qquad\forall a\in\mathbb{R}^{{d_{k}}};$

(6.1)

(wtedy $x_{k}=x_{0}+V_{k}a_{k}$ ).

Na mocy założenia istnieje macierz $B^{{1/2}}$ symetryczna i dodatnio określona taka, że $(B^{{1/2}})^{2}=B$ (dlatego nazywamy ją pierwiastkiem z $B$ ), a stąd wynika, że $\| x-x^{*}\| _{B}=\| B^{{1/2}}(x-x^{*})\| _{2}$ . Znaczy to, że zadanie minimalizacji (6.1) możemy zapisać w równoważnej postaci

$\displaystyle\| B^{{1/2}}(x_{0}+V_{k}a_{k}-x^{*})\| _{2}\leq\| B^{{1/2}}(x_{0}+V_{k}a-x^{*})\| _{2}\qquad\forall a\in\mathbb{R}^{{d_{k}}}.$

(6.2)

Jest to standardowe zadanie najmniejszych kwadratów,

$\| Ma_{k}-g\| _{2}=\min!,$

gdzie $M=B^{{1/2}}V_{k}$ jest macierzą prostokątną pełnego rzędu, natomiast $g=B^{{1/2}}(x^{*}-x_{0})$ . A zatem $x_{k}$ jest wyznaczone jednoznacznie.

Na mocy wniosku 6.1, dowolny $x\in x_{0}+K_{k}$ spełnia $x-x^{*}=p(A)(x_{0}-x^{*})$ , gdzie $p$ jest pewnym wielomianem stopnia co najwyżej $k$ takim, że $p(0)=1$ . Jeśli oznaczymy przez $\tilde{P}_{k}$ zbiór wielomianów stopnia co najwyżej $k$ o wyrazie wolnym równym 1, to warunek minimalizacji błędu możemy sformułować w równoważnej postaci

$\displaystyle\| x_{k}-x^{*}\| _{B}=\min _{{p\in\tilde{P}_{k}}}\| p(A)(x_{0}-x^{*})\| _{B}.$

(6.3)

W szczególności, dla wielomianu $p_{N}(\lambda)=\det(A-\lambda I)/\det(A)$ — będącego przeskalowanym wielomianem charakterystycznym macierzy $A$ — mamy, że $p_{N}(0)=1$ oraz oczywiście $p_{N}$ jest stopnia nie większego niż $N$ , a więc $p_{N}\in\tilde{P}_{N}$ . Z drugiej strony, na mocy twierdzenia Cayley'a–Hamiltona, $p_{N}(A)=0$ , skąd $\| x_{N}-x^{*}\| _{B}=0$ . Znaczy to, że przynajmniej dla $k=N$ zadanie minimalizacji ma (jednoznaczne) rozwiązanie, którym jest $x_{N}=x^{*}$ — a więc, że metoda ma własność skończonego stopu.

Dowód drugiej części twierdzenia, dotyczącej przypadku, gdy $x_{k}$ jest zadany warunkiem minimalizacji residuum, zostawiamy jako ćwiczenie.

∎

Ćwiczenie 6.3

Uzupełnij powyższy dowód.

Rozwiązanie:

Sposób I

Dla dowodu drugiego przypadku zauważmy, że $x_{k}$ znów jest dobrze określony jako rozwiązanie pewnego liniowego zadania najmniejszych kwadratów postaci

$\| B^{{1/2}}(r_{0}-AV_{k}a)\| _{2}=\min!$

oraz na mocy wniosku 6.1 zachodzi

$\| r_{k}\| _{B}=\min _{{p\in\tilde{P}_{k}}}\| p(A)(r_{0})\| _{B}.$

(6.4)

Dalszy ciąg dowodu jest identyczny jak poprzednio.

Sposób II

Oznaczmy chwilowo macierz indukującą normę energetyczną, w której minimalizujemy residuum, symbolem $\tilde{B}$ . Wystarczy skorzystać z pierwszej części twierdzenia, biorąc $B=A^{{-T}}\tilde{B}A^{{-1}}$ , gdyż

$\| r_{k}\| _{{\tilde{B}}}^{2}=r_{k}^{T}\tilde{B}r_{k}=(x^{*}-x_{k})^{T}A^{{-T}}\tilde{B}A^{{-1}}(x^{*}-x_{k})=\| x_{k}-x^{*}\| _{B}.$

Wniosek 6.3

Jeśli po $k$ -tej iteracji metody Kryłowa o własności skończonego stopu następuje stagnacja: $K_{k}=K_{{k+1}}$ , to znaczy, że właśnie znaleziono rozwiązanie dokładne, $x_{{k}}=x^{*}$ .

Dowód

Jeśli $K_{k}=K_{{k-1}}$ , to oczywiście $x_{{k-1}}=x_{{k}}=x_{{k+1}}=\ldots=x_{N}=\ldots$ . Ponieważ jednak musi być $x_{N}=x^{*}$ (skończony stop po $N$ iteracjach!) to oznacza, że $x_{{k-1}}=x^{*}$ .

∎

6.1. Metoda gradientów sprzężonych (CG)

Algorytm CG został uznany za jeden z 20 najważniejszych algorytmów numerycznych opracowanych w XX wieku.

W niniejszym rodziale będziemy zakładać, że kwadratowa macierz rzeczywista $A$ rozmiaru $N$ jest symetryczna, $A=A^{T}$ , oraz jest dodatnio określona,

$x^{T}Ax>0\quad\forall x\neq 0.$

Przy tym założeniu, można określić normę energetyczną indukowaną przez $A$ , zadaną tożsamością

$||y||_{A}^{2}=y^{T}Ay.$

Metodę gradientów sprzężonych, w skrócie CG (ang. conjugate gradients), zdefiniujemy początkowo w sposób niejawny. Kolejne przybliżenie $x_{k}$ określimy jako wektor z podprzestrzeni afinicznej $x_{0}+K_{k}$ , minimalizujący w tej podprzestrzeni błąd w normie energetycznej indukowanej przez $A$ :

$\displaystyle\| x_{k}-x^{*}\| _{A}\leq\| x-x^{*}\| _{A}\quad\forall x\in x_{0}+K_{k}.$

(6.5)

Naturalnie, taka definicja może budzić w nas nieco wątpliwości, co do jego obliczalności (w sformułowaniu warunku minimalizacji występuje szukane przez nas rozwiązanie dokładne, $x^{*}$ ).

Stwierdzenie 6.2 (CG jako metoda bezpośrednia)

Zadanie minimalizacji (6.5) ma jednoznaczne rozwiązanie. Jeśli $V_{k}$ jest macierzą, której kolumny tworzą bazę $K_{k}$ , to $x_{k}$ jest dane wzorem $x_{k}=x_{0}+V_{k}a_{k}$ , gdzie $a_{k}$ spełnia układ równań

$\displaystyle V_{k}^{T}AV_{k}a=V_{k}^{T}A(x^{*}-x_{0})=V_{k}^{T}r_{0}.$

(6.6)

Ponadto, w arytmetyce dokładnej, metoda CG znajduje dokładne rozwiązanie w co najwyżej $N$ iteracjach.

Dowód

Jest to natychmiastowy wniosek z twierdzenia 6.1, dla przypadku minimalizacji błędu gdy $B=A$ . Zależność (6.6) to nic innego jak układ równań normalnych dla zadania najmniejszych kwadratów (6.2).

∎

Z powyższego lematu wynika (por. (6.6)), że $x_{k}$ jest istotnie obliczalny: do jego wyznaczenia nie jest nam efektywnie potrzebna znajomość rozwiązania!

6.1.1. Implementacja

Aby wyznaczyć $x_{k}$ , nie będziemy bezpośrednio rozwiązywać układu (6.6) — byłoby to, wraz z postępem iteracji, coraz bardziej kosztowne, ze względu na zwiększający się rozmiar zadania najmniejszych kwadratów. Spróbujemy znaleźć tańszy sposób wyznaczania $x_{k}$ .

Ponieważ $A$ jest symetryczna, istnieje baza ortogonalna w $\mathbb{R}^{N}$ złożona z wektorów własnych $q_{1},\ldots,q_{N}$ :

$Aq_{i}=\lambda _{i}q_{i},\quad i=1,\ldots,N.$

Oznaczając przez $Q$ macierz, której kolejne kolumny są wektorami własnymi $A$ ,

$Q=\begin{pmatrix}q_{1}&|&q_{2}&|&\ldots&|&q_{N}\end{pmatrix},$

mamy, że $Q$ jest macierzą ortogonalną, $Q^{T}Q=I=QQ^{T}$ , a ponadto $A$ ma rozkład:

$A=Q\Lambda Q^{T},$

gdzie

$\Lambda=\begin{pmatrix}\lambda _{1}&&\\ &\ddots&\\ &&\lambda _{N}\end{pmatrix}.$

Gdyby baza przestrzeni $K_{k}$ była $A$ –ortogonalna (z powodów historycznych, jej elementy oznaczymy $p_{0},\ldots,p_{{k-1}}$ tak, że $V_{k}=\begin{pmatrix}p_{0}&|&p_{1}&|&\ldots&p_{{k-1}}\end{pmatrix}$ ), tzn. $p_{i}^{T}Ap_{j}=0$ dla $i\neq j$ , to wtedy macierz równań normalnych byłaby diagonalna,

$V_{k}^{T}AV_{k}=\begin{pmatrix}&\vdots&\\ \ldots&p_{i}^{T}Ap_{j}&\ldots\\ &\vdots&\\ \end{pmatrix}=\begin{pmatrix}p_{0}^{T}Ap_{0}&&\\ &\ddots&\\ &&p_{{k-1}}^{T}Ap_{{k-1}}\\ \end{pmatrix}.$

Wtedy kolejną iterację można wyznaczyć z jawnego wzoru:

$\displaystyle x_{{k}}=x_{0}+\sum _{{i=0}}^{{k-1}}\dfrac{p_{i}^{T}r_{0}}{p_{i}^{T}Ap_{i}}p_{i}.$

(6.7)

Zatem potrzebna jest nam skuteczna metoda wyznaczania bazy ortogonalnej w przestrzeni $K_{k}$ …. Oczywiście, ze względu na koszt obliczeniowy i pamięciowy, generowanie i następnie ortogonalizacja oryginalnego zestawu wektorów $\{ r_{0},Ar_{0},\ldots,A^{{k-1}}r_{0}\}$ rozpinających $K_{k}$ nie ma większego sensu. W zamian, wykorzystamy specjalne własności wektorów otrzymywanych w trakcie działania metody.

Lemat 6.1 (o ortogonalności residuów)

Residuum na $k$ -tym kroku, $r_{k}$ , jest prostopadłe do $K_{k}$ . Ponadto $r_{k}\in K_{{k+1}}$ .

Dowód

Uzasadnienie pierwszej części łatwo wynika z układu równań normalnych (6.6), określającego pośrednio $x_{k}$ . Rzeczywiście, ponieważ $x_{k}-x_{0}=V_{k}a_{k}$ , to z (6.6) wynika, że $V_{k}^{T}A(x_{k}-x_{0})=V_{k}^{T}A(x^{*}-x_{0})$ . Upraszczając wyrazy z $x_{0}$ , dostajemy $V_{k}^{T}A(x^{*}-x_{k})=0$ , czyli $V_{k}^{T}r_{k}=0$ .

Druga część wynika natychmiast z faktu, że $x_{k}\in x_{0}+K_{k}$ , skąd $Ax_{k}\in Ax_{0}+AK_{k}$ i w konsekwencji, odejmując stronami od $b$ , dochodzimy do $r_{k}\in r_{0}+AK_{k}$ . Tymczasem z definicji przestrzeni Kryłowa $r_{0}+AK_{k}\subseteq K_{{k+1}}$ .

∎

Z powyższego wynika, że jeśli $r_{{k-1}}\neq 0$ , to $K_{{k-1}}\subset K_{k}$ , a więc dopóki nie trafimy w rozwiązanie dokładne, $x^{*}$ , kolejne przestrzenie Kryłowa w metodzie $CG$ tworzą ściśle wstępujący ciąg przestrzeni, $K_{0}\subset K_{1}\subset\ldots K_{{k-1}}\subset K_{k}\subseteq\mathbb{R}^{N}$ .

W dalszm ciągu założymy więc, że $r_{{k-1}}\neq 0$ — a więc, że $x_{{k-1}}\neq x^{*}$ . Przypuśćmy, że mamy już zadaną bazę $A$ –ortogonalną $\{ p_{0},\ldots,p_{{k-1}}\}$ przestrzeni $K_{{k-1}}$ i znamy $x_{{k-1}}$ , $r_{{k-1}}$ . Naszym celem będzie wyznaczenie $x_{k}$ , $r_{k}$ oraz $p_{k}$ . Z zależności (6.7) mamy, że

$\displaystyle x_{k}=x_{{k-1}}+\alpha _{k}p_{{k-1}},$

(6.8)

gdzie

$\displaystyle\alpha _{k}=\dfrac{p_{{k-1}}^{T}r_{0}}{p_{{k-1}}^{T}Ap_{{k-1}}}.$

(6.9)

Obkładając (6.8) macierzą $A$ i odejmując obustronnie od $b$ dostajemy dodatkowo zależność rekurencyjną na residua,

$\displaystyle r_{k}=r_{{k-1}}-\alpha _{k}Ap_{{k-1}}.$

(6.10)

Potrzeba nam jeszcze zależności rekurencyjnej pozwalającej wyznaczyć $p_{{k}}$ — ostatni wektor bazy ortogonalnej dla $K_{{k+1}}$ (wcześniejsze znamy z założenia indukcyjnego). Ponieważ z lematu 6.1 wynika, że $K_{{k+1}}=\text{span}\{ r_{{k}},p_{0},\ldots,p_{{k-1}}\}$ , znaczy to, że

$p_{{k}}=r_{{k}}+\beta _{k}p_{{k-1}}+\gamma _{k}p_{{k-2}}+\ldots.$

Mnożąc skalarnie tę równość przez $Ap_{{k-1}}$ dostajemy z założenia $A$ -ortogonalności

$\displaystyle\beta _{k}=-\dfrac{r_{{k}}^{T}Ap_{{k-1}}}{p_{{k-1}}^{T}Ap_{{k-1}}}$

(6.11)

i podobnie, że $\gamma _{k}$ oraz wszystkie następne współczynniki są równe zero (ponieważ z lematu o residuach $r_{{k}}$ jest ortogonalne do $K_{{k}}$ , a $Ap_{j}\in K_{{k}}$ dla $j\leq k-2$ ). Zatem ostatecznie dostajemy kolejną elegancką zależność rekurencyjną, tym razem na wektory bazy $A$ –ortogonalnej dla $K_{{k+1}}$ :

$\displaystyle p_{{k}}=r_{{k}}+\beta _{k}p_{{k-1}}.$

(6.12)

Ponieważ $p_{0}=r_{0}$ , tym samym zależności (6.8)—(6.11) stanowią domknięty układ: startując z zadanego $x_{0}$ , jesteśmy w stanie wyznaczać kolejne przybliżenia.

Okazuje się, że powyższe wzory można jeszcze bardziej wymasować, otrzymując w końcu bardzo zwarty i tani algorytm:

Metoda CG

$r=b-Ax$ ;

$\rho _{0}=\| r\| _{2}^{2}$ , $\beta=0$ , $k=1$

while not stop

begin

$p=r+\beta p$

$w=Ap$

$\alpha=\dfrac{\rho _{{k-1}}}{p^{T}w}$

$x=x+\alpha p$

$r=r-\alpha w$

$\rho _{k}=\| r\| _{2}^{2}$

$\beta=\dfrac{\rho _{{k}}}{\rho _{{k-1}}}$

$k=k+1$

end

Jak widać, całą iterację da się wykonać, przechowując w pamięci tylko kilka wektorów (a nie, jak możnaby się obawiać, całą przestrzeń $K_{k}$ ), a najdroższym jej elementem jest mnożenie macierzy przez wektor.

Ćwiczenie 6.4

Opierając się na wzorach (6.8)—(6.11), wyprowadź powyższą postać algorytmu CG.

Rozwiązanie:

Najpierw wykażemy, że

$\displaystyle\alpha _{k}=\dfrac{\| r_{{k-1}}\| _{2}^{2}}{p_{{k-1}}^{T}Ap_{{k-1}}}.$

(6.13)

Ponieważ z założenia $x_{j}=x_{0}+\sum _{{i=0}}^{{j-1}}\alpha _{{i+1}}p_{i}$ dla każdego $j\leq k$ , to (mnożąc obustronnie tę równość przez $A$ i odejmując stronami od $b$ ) zachodzi także $r_{{j}}=r_{0}-\sum _{{i=0}}^{{j-1}}\alpha _{{i+1}}Ap_{i}$ . Mnożąc skalarnie tę równość przez $p_{{j}}$ i uwzględniając $A$ –ortogonalność kierunków $p_{i}$ dochodzimy do wniosku, że

$p_{{j}}^{T}r_{{j}}=p_{{j}}^{T}r_{{0}}.$

Z drugiej zaś strony, mnożąc (6.12) skalarnie przez $r_{{j}}$ otrzymujemy

$p_{{j}}^{T}r_{{j}}=\| r_{{j}}\| _{2}^{2}-\beta _{{j-1}}p_{{j-1}}^{T}r_{{j}}=\| r_{{j}}\| _{2}^{2},$

ponieważ $r_{{j}}$ jest prostopadłe do $K_{{j}}$ , w której zawarty jest wektor $p_{{j-1}}$ . Ostatecznie więc $p_{j}^{T}r_{{0}}=\| r_{{j}}\| _{2}^{2}$ dla każdego $j\leq k$ . Biorąc $j=k-1$ , z (6.9) otrzymujemy (6.13).

Teraz wyprowadzimy prostszą reprezentację współczynnika $\beta _{k}$ . Z rekurencyjnej zależności pomiędzy residuami (6.10) wynika, że

$r_{{k}}^{T}r_{k}=r_{k}^{T}r_{{k-1}}-\alpha _{k}r_{k}^{T}Ap_{{k-1}}.$

Ponieważ z lematu o ortogonalności residuów $r_{{k-1}}\in K_{{k}}$ oraz $r_{k}$ jest ortogonalne do $K_{k}$ , to $r_{k}^{T}r_{{k-1}}=0$ , więc podstawiając do powyższego wzoru uzyskane przed chwilą nowe wyrażenie na $\alpha _{k}$ dostajemy

$\| r_{{k}}\| _{2}^{2}=-\dfrac{\| r_{{k-1}}\| _{2}^{2}}{p_{{k-1}}^{T}Ap_{{k-1}}}r_{k}^{T}Ap_{{k-1}}=\| r_{{k-1}}\| _{2}^{2}\beta _{k}.$

Stąd i z (6.13) już wynika wzór na współczynnik $\beta _{k}$ ,

$\beta _{k}=\dfrac{\| r_{{k}}\| _{2}^{2}}{\| r_{{k-1}}\| _{2}^{2}}.$

Dla dużych $N$ , traktowanie CG jako metody bezpośredniej nie miałoby większego sensu — nie dość, że wykonanie aż $N$ iteracji mogłoby być zadaniem ponad możliwości naszego komputera, to jeszcze dodatkowo algorytm wykorzystuje bardzo specyficzne relacje pomiędzy wektorami, a całość jest przecież w praktyce realizowana w arytmetyce zmiennoprzecinkowej o ograniczonej precyzji, w której te relacje nie zachodzą (w sposób dokładny). Prowadzi to do tego, że w miarę postępu iteracji na przykład wektory $p_{k}$ są coraz mniej ortogonalne i tym samym metoda nie musi dotrzeć do dokładnego rozwiązania.

Dlatego w praktyce znacznie bardziej właściwe wydaje się potraktowanie metody CG (i innych metod Kryłowa) jako ,,czystej” metody iteracyjnej i oszacowanie szybkości redukcji błędu podobnie, jak czyniliśmy to w przypadku metod stacjonarnych.

Twierdzenie 6.2 (o zbieżności CG jako metody iteracyjnej)

Po $k$ iteracjach metody CG,

$||x_{k}-x||_{A}\leq 2\,\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{k}\,||x_{0}-x||_{A},$

gdzie $\kappa=\operatorname{cond}_{2}(A)=\lambda _{{\max}}(A)/\lambda _{{\min}}(A)$ .

Dowód

Skorzystamy z własności (6.3). Zauważmy, że

$\| p_{k}(A)(x^{*}-x_{0})\| _{A}\leq\| p_{k}(A)\| _{2}\|(x^{*}-x_{0})\| _{A}$

oraz

$\| p_{k}(A)\| _{2}=\max _{{\lambda\in\sigma(A)}}|p_{k}(\lambda)|,$

zatem wystarczy oszacować wartości wybranego wielomianu $p_{k}\in\tilde{P}_{k}$ (nasza norma błędu jest i tak nie większa). Niech $M=\lambda _{{\max}}(A)$ i $m=\lambda _{{\min}}(A)$ . Jako $p_{k}$ w (6.3) weźmy przeskalowany $k$ -ty wielomian Czebyszewa,

$p_{k}(z)=\dfrac{T_{k}\left(\dfrac{M+m-2z}{M-m}\right)}{T_{k}\left(\dfrac{M+m}{M-m}\right)}.$

Rzeczywiście, $p_{k}\in\tilde{P}_{k}$ . Ponadto, ponieważ wielomiany Czebyszewa spełniają zależność

$|T_{k}(x)|\leq 1\qquad\text{ dla }x\in[-1,1],$

$\max _{{z\in[m,M]}}|p_{k}(z)|\leq\dfrac{1}{T_{k}\left(\dfrac{M+m}{M-m}\right)}=\dfrac{1}{T_{k}\left(\dfrac{\kappa+1}{\kappa-1}\right)}.$

Należy więc oszacować $T_{k}\left(\dfrac{\kappa+1}{\kappa-1}\right)$ . Ponieważ $\dfrac{\kappa+1}{\kappa-1}>1$ , to skorzystamy ze wzoru

$T_{k}(x)=\dfrac{1}{2}(x+\sqrt{x^{2}-1})^{k}+\dfrac{1}{2}(x-\sqrt{x^{2}-1})^{k}\qquad\text{ dla }|x|\geq 1.$

W szczególności więc, biorąc $x=\dfrac{\kappa+1}{\kappa-1}$ mamy

$T_{k}(x)\geq\dfrac{1}{2}(x+\sqrt{x^{2}-1})^{k}=\dfrac{1}{2}\left(\dfrac{\sqrt{\kappa}+1}{\sqrt{\kappa}-1}\right)^{k}.$

∎

Stwierdzenie 6.3

Jeśli macierz $A$ ma $m$ różnych wartości własnych, to metoda CG w arytmetyce dokładnej znajdzie rozwiązanie dokładne $x^{*}$ w co najwyżej $m$ iteracjach.

Ćwiczenie 6.5

Udowodnij powyższe stwierdzenie.

Wskazówka:

Rozważ odpowiednio przeskalowany wielomian $p(\lambda)=(\lambda-\lambda _{1})\cdots(\lambda-\lambda _{m})$ .

Przykład 6.1

Kontynuujemy przykład 5.13. Chcąc porównywać cztery metody: Jacobiego, SOR, metodę najszybszego spadku oraz sprzężonych gradientów, będziemy korzystać z macierzy

$A=B^{T}B+pI,$

gdzie $p>0$ oraz $B$ jest losową macierzą rozrzedzoną. Zwiększanie parametru $p$ nie tylko poprawia diagonalną dominację, ale także poprawia uwarunkowanie $A$ . Jako parametr relaksacji dla SOR wybraliśmy (strzelając w ciemno) $\omega=1.3$ .

%
% rozwiazania zadania z macierza rozrzedzona
% 
disp('Matematyka obliczeniowa II');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = jacobi(A, b, tol=1e-8, maxit=100, x0)
%
% metoda Jacobiego na macierzy A
%
if (nargin < 5)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy: 
% A = M - Z
M = diag(A); Z = tril(A,-1) + triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = (b - Z*x)./M;
	iter = iter+1;
	resid(iter) = norm(b-A*x);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = sor(A, b, omega=1, tol=1e-8, maxit=100, x0)
%
% metoda SOR (gdy omega=1, to Gaussa-Seidela) na macierzy A
%
if (nargin < 6)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy: 
% A = M - Z;
M = tril(A,-1) + spdiag(diag(A)/omega);
Z = spdiag((1/omega - 1)*diag(A)) - triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = M \ (Z*x + b);
	iter = iter+1;
	resid(iter) = norm(b-A*x);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = steepestdescent(A, b, tol=1e-8, maxit=100, x0)
%
% metoda najszybszego spadku na macierzy A
% wymaga A = A' > 0
%
if (nargin < 5)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0; r = b - A*x;
resid(iter) = norm(r);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	w = A*r;
	a = (r'*r)/(r'*w);
	x = x + a*r;
	r = r - a*w;
	iter = iter+1;
	resid(iter) = norm(r);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 250; % rozmiar macierzy
b = rand(N,1);

% testy dla roznego stopnia diagonalnej dominacji macierzy A

prange = [1, 5, 10, 50]; % zestaw wartosci testowych parametru p
B = sprand(N, N, 10/N); B = B'*B;
% e = ones(N,1); B = spdiags([-e, 2*e, -e], [-1,0,1], N,N);
spy(B); pause(0.5);

omega = 1.3; % parametr dla SOR
for k = 1:length(prange)
	p = prange(k);
	A = B + p*speye(N);
	tic; [x, info, its, resid] = jacobi(A, b); 
	iter(k) = its; time(k) = toc;
	tic; [x, info, its, residg] = sor(A, b, omega); 
	iterg(k) = its; timeg(k) = toc;
	tic; [x, info, its, resids] = steepestdescent(A, b); 
	iters(k) = its; times(k) = toc;
	tic; [x, info, relres, its, residc] = pcg(A, b, 1e-8, 100); 
	iterc(k) = its; timec(k) = toc; 
	tmp = zeros(size(resid)); tmp(1:length(residc)) = residc; residc=tmp; 
	mi = max([iter(k),iterg(k),iters(k),iterc(k)]);
	semilogy(1:mi,[resid(1:mi), residg(1:mi), resids(1:mi), residc(1:mi)],[['b-;p=',num2str(p),' Jacobi;']; ['r-;p=',num2str(p),' (SOR(',num2str(omega),');']; ['g-;p=',num2str(p),' (najszybszy spadek);']; ['k-;p=',num2str(p),' (CG);']]); pause(0.5);
end

% dla porownania - solver bezposredni

tic; x = A\b; timed = toc;

printf('Metoda Jacobiego (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d | %4d | %10.5f\n', [prange; iter; time]);
printf('Metoda SOR(%f) (p, iter, czas)\n', omega);
printf('  p  | iter |   czas\n====================================\n');
printf('%4d | %4d | %10.5f\n', [prange; iterg; timeg]);
printf('Metoda najszybszego spadku (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d | %4d | %10.5f\n', [prange; iters; times]);
printf('Metoda CG (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d | %4d | %10.5f\n', [prange; iterc; timec]);
printf('Metoda bezposrednia (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d %4d %10.5f\n', [prange(end); 0; timed]);
printf('\nRozmiar zadania N = %d\n', N);

Zwróćmy uwagę na wyraźną przewagę metody CG nad pozostałymi. Sprawdź, czy podobnie jest dla większych wartości $N$ .

Ćwiczenie 6.6

Sprawdź w przykładzie 6.1, czy faktycznie uwarunkowanie macierzy $A$ wpływa na szybkość zbieżności metody CG i najszybszego spadku. Aby zbadać uwarunkowanie macierzy, możesz skorzystać z polecenia cond(A), albo wykorzystać estymator uwarunkowania dostępny w pcg.

Ćwiczenie 6.7

Sprawdź, modyfikując kod przykładu 6.1, czy jeśli $A$ nie będzie symetryczna (lub nie będzie dodatnio określona), wpłynie to istotnie na szybkość zbieżności metody CG i najszybszego spadku. Wypróbuj m.in. $A=B+pI$ dla $p>0$ (brak symetrii) tak dobranego, by $A_{{\operatorname{sym}}}>0$ oraz $A=B^{T}B+pI$ dla $p<0$ takiego, żeby $A$ miało i dodatnie, i ujemne wartości własne.

Przykład 6.2

Chcąc porównywać cztery metody: Jacobiego, SOR, metodę najszybszego spadku oraz sprzężonych gradientów dla macierzy jednowymiarowego laplasjanu $T_{N}$ . Jako parametr relaksacji dla SOR wybraliśmy wartość optymalną, zgodnie z przykładem 5.11.

%
% rozwiazania zadania z macierza rozrzedzona
% 
disp('Matematyka obliczeniowa II');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = jacobi(A, b, tol=1e-10, maxit=1000, x0)
%
% metoda Jacobiego na macierzy A
%
if (nargin < 5)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy: 
% A = M - Z
M = diag(A); Z = tril(A,-1) + triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = (b - Z*x)./M;
	iter = iter+1;
	resid(iter) = norm(b-A*x);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = sor(A, b, omega=1, tol=1e-10, maxit=1000, x0)
%
% metoda SOR (gdy omega=1, to Gaussa-Seidela) na macierzy A
%
if (nargin < 6)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy: 
% A = M - Z;
M = tril(A,-1) + spdiag(diag(A)/omega);
Z = spdiag((1/omega - 1)*diag(A)) - triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = M \ (Z*x + b);
	iter = iter+1;
	resid(iter) = norm(b-A*x);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = steepestdescent(A, b, tol=1e-8, maxit=1000, x0)
%
% metoda najszybszego spadku na macierzy A
% wymaga A = A' > 0
%
if (nargin < 5)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0; r = b - A*x;
resid(iter) = norm(r);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	w = A*r;
	a = (r'*r)/(r'*w);
	x = x + a*r;
	r = r - a*w;
	iter = iter+1;
	resid(iter) = norm(r);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nrange = [10, 100, 200, 400]; % zestaw wartosci testowych parametru p
for k = 1:length(Nrange)
	N = Nrange(k); % rozmiar macierzy
	b = rand(N,1);
	e = ones(N,1); A = spdiags([-e, 2*e, -e], [-1,0,1], N,N);
	spy(A); pause(0.5);
	
	omega = 2/(1+sin(pi/(N+1))); % optymalny parametr dla SOR, gd macierz jest 1 wym laplasjanu
	% omega = 1.3;
	tic; [x, info, its, resid] = jacobi(A, b); 
	iter(k) = its; time(k) = toc;
	tic; [x, info, its, residg] = sor(A, b, omega); 
	iterg(k) = its; timeg(k) = toc; omegag(k) = omega;
	tic; [x, info, its, resids] = steepestdescent(A, b); 
	iters(k) = its; times(k) = toc;
	tic; [x, info, relres, its, residc] = pcg(A, b, 1e-8, 1000); 
	iterc(k) = its; timec(k) = toc; 
	tmp = zeros(size(resid)); tmp(1:length(residc)) = residc; residc=tmp; 
	mi = max([iter(k),iterg(k),iters(k),iterc(k)]);
	semilogy(1:mi,[resid(1:mi), residg(1:mi), resids(1:mi), residc(1:mi)],[['b-;N=',num2str(N),' Jacobi;']; ['r-;N=',num2str(N),' (SOR(',num2str(omega),');']; ['g-;N=',num2str(N),' (najszybszy spadek);']; ['k-;N=',num2str(N),' (CG);']]); pause(1.5);
end

% dla porownania - solver bezposredni

tic; x = A\b; timed = toc;

Zwróćmy uwagę na wyraźną przewagę metody CG nad pozostałymi metodami iteracyjnymi. Jednak i tak nie wytrzymuje ona konkurencji z metodą bezpośrednią. Ta sytuacja dramatycznie zmieni się, gdy będziemy rozważali dyskretyzacje dwu- lub trójwymiarowego operatora Laplace'a. O ile wtedy metoda CG wciąż ma kłopoty z szybką zbieżnością, to metoda bezpośrednia (typu rozkładu LU) staje się całkowicie bezużyteczna.

6.2. Metoda GMRES

Metoda GMRES (ang. Generalized Minimum RESidual) nie wymaga ani symetrii, ani dodatniej określoności macierzy, jest więc bardziej uniwersalna, choć też w realizacji bardziej kosztowna od CG. Jej szczegółowe omówienie, w tym — wyprowadzenie subtelniejszych oszacowań szybkości zbieżności — wykracza niestety poza ramy niniejszego wykładu (zainteresowani mogą skonsultować m.in. podręcznik [13]).

W tej metodzie, przybliżenie na $k$ -tej iteracji ma minimalizować, w przestrzeni afinicznej $x_{0}+K_{k}$ , residuum (stąd nazwa) mierzone w normie euklidesowej:

$\displaystyle\| r_{{k}}\| _{2}=\| b-Ax_{{k}}\| _{2}\leq\| b-Ax\| _{2}\quad\forall x\in x_{0}+K_{k}.$

(6.14)

Stwierdzenie 6.4 (GMRES jako metoda bezpośrednia)

Zadanie minimalizacji (6.14) ma jednoznaczne rozwiązanie. Jeśli $V_{k}$ jest macierzą, której kolumny tworzą bazę $K_{k}$ , to $x_{k}$ jest dane wzorem $x_{k}=x_{0}+V_{k}a_{k}$ , gdzie $a_{k}$ jest rozwiązaniem zadania najmniejszych kwadratów względem $a\in\mathbb{R}^{{\dim K_{k}}}$ ,

$\displaystyle\| r_{0}-AV_{k}a\| _{2}=\min!$

(6.15)

Ponadto, w arytmetyce dokładnej, metoda GMRES znajduje rozwiązanie $x^{*}$ w co najwyżej $N$ iteracjach.

Dowód

Jest to bezpośredni wniosek z twierdzenia 6.1, dla przypadku minimalizacji residuum gdy $B=I$ .

∎

Twierdzenie 6.3 (o zbieżności metody GMRES jako metody iteracyjnej)

Jeśli macierz $A$ jest diagonalizowalna (to znaczy $A=X\Lambda X^{{-1}}$ i $\Lambda$ jest macierzą diagonalną⁵Oczywiście, na diagonali $\Lambda$ znajdują się wartości własne macierzy $A$ . W ogólności więc, zarówno $X$ , jak i $\Lambda$ mogą być zespolone.), to $k$ -ta iteracja GMRES spełnia

$\| r_{k}\| _{2}\leq\operatorname{cond}_{2}(X)\min _{{p\in\tilde{P}_{k}}}\max _{{\lambda\in\sigma(A)}}|p(\lambda)|.$

Dowód

Z (6.4) wynika, że

$\| r_{k}\| _{2}\leq\min _{{p\in\tilde{P}_{k}}}\| p(A)\| _{2}\,\| r_{0}\| _{2},$

więc wystarczy oszacować $\| p(A)\| _{2}$ :

$\| p(A)\| _{2}=\| p(X\Lambda X^{{-1}})\| _{2}=\| Xp(\Lambda)X^{{-1}}\| _{2}\leq\| X\| _{2}\| p(\Lambda)\| _{2}\| X^{{-1}}\| _{2}=\operatorname{cond}_{2}(X)\max _{{\lambda\in\sigma(A)}}|p(\lambda)|.$

∎

Powyższe twierdzenie jest mało praktyczne, ze względu na to, że zazwyczaj nie znamy macierzy $X$ . Ale jeśli na przykład macierz $A$ jest normalna, to znaczy $A^{T}A=AA^{T}$ , to wtedy $X^{{-1}}=X^{T}$ i w konsekwencji $\operatorname{cond}_{2}(X)=1$ .

Ćwiczenie 6.8

Niech $A$ będzie macierzą diagonalizowalną. Wykaż, że jeśli $A$ ma $m$ różnych wartości własnych, to GMRES w arytmetyce dokładnej osiągnie rozwiązanie w co najwyżej $m$ krokach.

Rozwiązanie:

Wynika z poprzedniego twierdzenia, wystarczy wziąć $p(z)=(z-\lambda _{1})\cdots(z-\lambda _{m})/\lambda _{1}\cdots\lambda _{m}$ .

6.2.1. Implementacja

Dla ostatecznej skuteczności metody iteracyjnej ważna jest także jej efektywna implementacja. Załóżmy więc, że z powodzeniem wykonaliśmy $k-1$ kroków metody GMRES i teraz chcemy wyznaczyć $x_{k}$ . Oznaczamy, jak zwykle, przez $V_{k}$ macierz, której kolumny tworzą bazę $K_{k}$ . Wtedy, jeśli nie osiągnięto jeszcze rozwiązania, to $\dim K_{k}=k$ , zatem $V_{k}\in\mathbb{R}^{{N\times k}}$ . Aby wyznaczyć współczynniki $a_{k}$ reprezentacji $x_{k}$ w bazie $V_{k}$ , $x_{k}=x_{0}+V_{k}a_{k}$ , musimy rozwiązać względem $a\in\mathbb{R}^{k}$ liniowe zadanie najmniejszych kwadratów,

$\displaystyle\| r_{0}-AV_{k}a\| _{2}=\min!$

(6.16)

Wydawać by się mogło, że sprawa jest tu prosta, bo $AV_{k}$ możemy w miarę łatwo wyznaczyć, a potem należałoby skorzystać z jednej z metod rozwiązywania standardowego zadania najmniejszych kwadratów. Jednak jeśli uwzględnić drobny fakt, że do wyznaczenia $x_{k}$ potrzebna nam będzie znajomość nie samego $a_{k}$ , ale także $V_{k}$ , sytuacja robi się nieco niezręczna. Dlatego po raz kolejny spróbujemy spreparować bazę $V_{k}$ i jednocześnie wykorzystać fakt, że $AV_{k}$ rozpina podprzestrzeń w $K_{{k+1}}$ (rozpiętej przez $V_{{k+1}})$ .

Zwróćmy uwagę na to, że jeśli $\{ v_{1},\ldots,v_{j}\}$ są bazą ortogonalną w $K_{j}$ dla $j=1,\ldots,k$ , to aby rozszerzyć ją do $K_{{k+1}}$ wystarczy zortogonalizować — zamiast $A^{k}r_{0}$ — wektor $Av_{k}$ . Stosując zmodyfikowany algorytm ortogonalizacji Grama-Schmidta do $Av_{k}$ :

Metoda Arnoldiego wyznaczania bazy ortogonalnej

$v_{1}=r_{0}/\| r_{0}\| _{2}$

for $j=1$ to $k-1$ begin

$h_{{ij}}=v_{i}^{T}Av_{j}$ , dla $i=1,\ldots,j$

$v_{{j+1}}=Av_{j}-\sum _{{i=1}}^{j}h_{{ij}}v_{i}$ { $Aw_{j}$ to wektor, który będziemy ortogonalizować}

$h_{{i+1,j}}=\| v_{{j+1}}\| _{2}$

$v_{{j+1}}=v_{{j+1}}/h_{{i+1,j}}$ {unormowanie $v_{{j+1}}$ }

end

dostajemy wprost z relacji $Av_{k}=\sum _{{i=1}}^{{j+1}}h_{{ij}}v_{i}$ związek:

$\displaystyle AV_{k}=V_{{k+1}}H_{k},$

(6.17)

gdzie $H_{k}$ jest macierzą Hessenberga rozmiaru $(k+1)\times k$ :

Opisana tu zmodyfikowana metoda Grama–Schmidta ortogonalizacji bazy przestrzeni Kryłowa nosi nazwę metody Arnoldiego.

Jest to istotny postęp na drodze ku efektywnemu rozwiązaniu zadania minimalizacji (6.16), gdyż podstawiając doń (6.17) mamy, że

$\| r_{0}-AV_{k}a\| _{2}=\| r_{0}-V_{{k+1}}H_{k}a\| _{2}=\| V_{{k+1}}^{T}r_{0}-H_{k}a\| _{2}$

— w ostatniej równości skorzystaliśmy z ortogonalności $V_{{k+1}}$ . Ale przecież $r_{0}=\beta v_{1}$ (gdzie $\beta=\| r_{0}\| _{2}$ ), zaterm $V_{{k+1}}^{T}r_{0}=(\beta,0,\ldots,0)^{T}=\beta e_{1}$ ! Tak więc, zadanie najmniejszych kwadratów (6.16) dla macierzy rozmiaru $N\times k$ sprowadziliśmy do zadania najmniejszych kwadratów:

$\displaystyle\|\beta e_{1}-H_{k}a\| _{2}=\min!$

(6.18)

dla macierzy Hessenberga $H_{k}$ , rozmiaru $(k+1)\times k$ (a więc: prawie–trójkątnej!)⁶Macierz Hessenberga można — przy pomocy (ortogonalnych) przekształceń Givensa — doprowadzić liniowym kosztem do postaci górnej trójkątnej, zob. [9, rozdział 3.5]..

Stąd dostajemy bazową wersję implementacji algorytmu GMRES:

Metoda GMRES, wersja bazowa

$r_{0}=b-Ax_{0}$ , $V_{1}=r_{0}/\| r_{0}\| _{2}$ , k = 0

while not stop

{wykonaj ortogonalizację metodą Arnoldiego}

$h_{k}=V_{k}^{T}Av_{k}$

$v_{{k+1}}=Av_{k}-V_{k}h_{k}$

$h_{{k+1,k}}=\| v_{{k+1}}\| _{2}$

$v_{{k+1}}=v_{{k+1}}/h_{{k+1,k}}$

{rozszerz bazę $V_{k}$ i macierz $H_{k}$ }

$V_{{k+1}}=\begin{pmatrix}V_{k}&v_{{k+1}}\end{pmatrix}$

$H_{k}=\begin{pmatrix}H_{{k-1}}&h_{k}\\ &h_{{k+1,k}}\end{pmatrix}$

{znajdź minimum zadania najmniejszych kwadratów}

$r_{{k+1}}=\min _{a}\|\beta e_{1}-H_{k}a\| _{2}$

{zwiększ licznik iteracji}

$k=k+1$

end

{rozwiąż ostatnie zadanie najmniejszych kwadratów}

$a_{{k-1}}=\arg\min _{a}\|\beta e_{1}-H_{{k-1}}a\| _{2}$

$x_{{k}}=x_{0}+V_{{k-1}}a_{{k-1}}$ {dopiero teraz musimy wyznaczyć przybliżone $x$ }

Zwróćmy uwagę na pewien drobiazg obniżający koszt iteracji: dopóki nie wyznaczymy rozwiązania z żądaną dokładnością (co możemy określić jedynie poprzez analizę residuum, $r_{k}$ ), nie musimy znajdować ani $a_{k}$ , ani $x_{k}$ . Dzięki temu, koszt jednej iteracji GMRES jest tylko rzędu $O(kN)$ flopów (gdybyśmy wyznaczali $a_{k}$ na każdej iteracji, dodawałoby to dodatkowe $O(k^{2})$ flopów).

Powyższy algorytm w praktyce jest poddawany pewnym drobnym modyfikacjom, związanym z tym, że choć zmodyfikowany algorytm Grama–Schmidta ma numeryczne własności wyraźnie lepsze od klasycznego, to wciąż w realizacji w arytmetyce skończonej precyzji będzie tracił ortogonalność kolumn $V_{k}$ . Dlatego uzupełnia się go o warunkową reortogonalizację $V_{k}$ : szczegóły opisane są w [9] (tam także są przykłady ukazujące wagę tego dodatkowego kroku) oraz w [13].

Ćwiczenie 6.9

Przeanalizuj zbieżność metody GMRES dla macierzy kwadratowej

$A=\begin{pmatrix}\alpha I&X\\ 0&\beta I\end{pmatrix},$

gdzie $\alpha,\beta\in\mathbb{R}\setminus\{ 0\}$ , a $X$ jest macierzą $n\times m$ .

6.2.2. Wersja z restartem: GMRES( $m$ )

Pewną wadą metody GMRES jest konieczność zapamiętania na $k$ -tym kroku wszystkich wektorów bazy ortogonalnej przestrzeni Kryłowa $K_{k}$ . Znaczy to, że nakład obliczeń i, co gorsza, ilość pamięci potrzebnej na wykonanie jednego kroku metody rosną wraz z liczbą iteracji. Aby zniwelować tę przypadłość, zwykle stosuje się restarty metody GMRES co $m$ (rzędu kilkunastu–kilkudziesięcu) iteracji: po $m$ iteracjach przerywane są obliczenia, a wyznaczone w ten sposób rozwiązanie przybliżone jest używane jako początkowe przybliżenie do kolejnego wywołania metody. W ten sposób co prawda spowalnia się nieco szybkość zbieżności samej metody, ale za to utrzymujemy w ryzach jej zasobożerność: to często jest rozsądny kompromis.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Matematyka obliczeniowa II

Zagadnienia

6. Metody iteracyjne Kryłowa

Stwierdzenie 6.1

Wniosek 6.1

Ćwiczenie 6.1

Ćwiczenie 6.2

Twierdzenie 6.1 (metody Kryłowa oparte na minimalizacji błędu albo residuum w normie energetycznej)

Dowód

Ćwiczenie 6.3

Wniosek 6.3

Dowód

6.1. Metoda gradientów sprzężonych (CG)

Stwierdzenie 6.2 (CG jako metoda bezpośrednia)

Dowód

6.1.1. Implementacja

Lemat 6.1 (o ortogonalności residuów)

Dowód

Metoda CG

Ćwiczenie 6.4

Twierdzenie 6.2 (o zbieżności CG jako metody iteracyjnej)

Dowód

Stwierdzenie 6.3

Ćwiczenie 6.5

Quiz 6.1

Przykład 6.1

Ćwiczenie 6.6

Ćwiczenie 6.7

Przykład 6.2

6.2. Metoda GMRES

Stwierdzenie 6.4 (GMRES jako metoda bezpośrednia)

Dowód

Twierdzenie 6.3 (o zbieżności metody GMRES jako metody iteracyjnej)

Dowód

Ćwiczenie 6.8

6.2.1. Implementacja

Metoda Arnoldiego wyznaczania bazy ortogonalnej

Metoda GMRES, wersja bazowa

Ćwiczenie 6.9

6.2.2. Wersja z restartem: GMRES()

6.2.2. Wersja z restartem: GMRES( $m$ )