Zagadnienia

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:

  1. minimalizacji błędu:

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

    albo

  2. 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=0k=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],

to

\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}).

Quiz 6.1
  1. najszybszego spadku



    Dobrze Szybkość zbieżności metody najszybszego spadku jest proporcjonalna do współczynnika uwarunkowania macierzy

    Źle Zastanów się, jak szybkość zbieżności metody zależy od współczynnika uwarunkowania

  2. gradientów sprzężonych



    Dobrze Szybkość zbieżności metody CG jest proporcjonalna do pierwiastka współczynnika uwarunkowania macierzy

    Źle Zastanów się, jak szybkość zbieżności metody zależy od współczynnika uwarunkowania

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.



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.



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ą5Oczywiś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!)6Macierz 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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.