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ń A⁢x=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, x0, jest dowolnie zadanym wektorem w RN (na przykład, x0=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 xk będziemy dobierać w taki sposób, by xk∈x0+Kk⁢r0,A oraz spełniało pewien dodatkowy warunek, na przykład — minimalizowało pewną miarę błędu na przestrzeni Kryłowa

Kk⁢r0,A=span⁢r0,A⁢r0,…,Ak-1⁢r0,

gdzie r0=b-A⁢x0 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 Kk zamiast Kk⁢r0,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:

K0⊆K1⊆…⊆Kk⊆Kk+1⊆…⊆RN.

Łatwo pokazać, że nie tylko same przybliżenia xk, 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∈x0+Kk⁢r0,A, to x-x*∈Kk+1⁢x0-x*,A oraz b-A⁢x∈Kk+1⁢r0,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∈x0+Kk⁢r0,A, to błąd x-x* jest postaci I+c1⁢A+…+ck⁢Ak⁢x0-x* i analogicznie residuum b-A⁢x jest postaci I-c1⁢A-…-ck⁢Ak⁢r0, gdzie c1,…,ck są pewnymi rzeczywistymi współczynnikami.

Ćwiczenie 6.1

Udowodnij powyższe stwierdzenie i wniosek.

Wskazówka:

Wystarczy zauważyć, że x∈x0+Kk⁢r0,A jest postaci x=x0+∑i=0k-1ci+1⁢Ai⁢r0.

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 xB będziemy oznaczali normę energetyczną wektora x, indukowaną przez symetryczną, dodatnio określoną macierz B,

xB=xT⁢B⁢x.

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

Niech B∈RN będzie pewną macierzą symetryczną i dodatnio określoną i niech k-ta iteracja metody Kryłowa, xk∈x0+Kk, będzie określona przez jeden z warunków:

minimalizacji błędu:

xk-x*B≤x-x*B⁢⁢∀x∈x0+Kk,

albo
minimalizacji residuum:

rkB≤b-A⁢xB⁢⁢∀x∈x0+Kk.

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 xk jest dobrze określony przez warunek minimalizacji, gdyż jest zdefiniowany jako rozwiązanie pewnego liniowego zadania najmniejszych kwadratów.

Jeśli przez Vk oznaczymy macierz, której kolumnami są wektory tworzące bazę Kk, to każdy x∈x0+Kk daje się zapisać w postaci x=x0+Vk⁢a, gdzie a jest wektorem o tylu współrzędnych, jaki jest wymiar Kk (oznaczmy go dk=dim⁡Kk). Podstawiając do minimalizowanego wyrażenia, dostajemy zadanie znalezienia ak∈Rdk takiego, że

x0+Vk⁢ak-x*B≤x0+Vk⁢a-x*B⁢⁢∀a∈Rdk;

(6.1)

(wtedy xk=x0+Vk⁢ak).

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

B1/2⁢x0+Vk⁢ak-x*2≤B1/2⁢x0+Vk⁢a-x*2⁢⁢∀a∈Rdk.

(6.2)

Jest to standardowe zadanie najmniejszych kwadratów,

M⁢ak-g2=min!,

gdzie M=B1/2⁢Vk jest macierzą prostokątną pełnego rzędu, natomiast g=B1/2⁢x*-x0. A zatem xk jest wyznaczone jednoznacznie.

Na mocy wniosku 6.1, dowolny x∈x0+Kk spełnia x-x*=p⁢A⁢x0-x*, gdzie p jest pewnym wielomianem stopnia co najwyżej k takim, że p⁢0=1. Jeśli oznaczymy przez 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

xk-x*B=minp∈P~k⁡p⁢A⁢x0-x*B.

(6.3)

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

Dowód drugiej części twierdzenia, dotyczącej przypadku, gdy xk 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 xk znów jest dobrze określony jako rozwiązanie pewnego liniowego zadania najmniejszych kwadratów postaci

B1/2⁢r0-A⁢Vk⁢a2=min!

oraz na mocy wniosku 6.1 zachodzi

rkB=minp∈P~k⁡p⁢A⁢r0B.

(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 B~. Wystarczy skorzystać z pierwszej części twierdzenia, biorąc B=A-T⁢B~⁢A-1, gdyż

rkB~2=rkT⁢B~⁢rk=x*-xkT⁢A-T⁢B~⁢A-1⁢x*-xk=xk-x*B.

Wniosek 6.3

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

Dowód

Jeśli Kk=Kk-1, to oczywiście xk-1=xk=xk+1=…=xN=…. Ponieważ jednak musi być xN=x* (skończony stop po N iteracjach!) to oznacza, że xk-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=AT, oraz jest dodatnio określona,

xT⁢A⁢x>0⁢⁢∀x≠0.

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

yA2=yT⁢A⁢y.

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

xk-x*A≤x-x*A⁢⁢∀x∈x0+Kk.

(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 Vk jest macierzą, której kolumny tworzą bazę Kk, to xk jest dane wzorem xk=x0+Vk⁢ak, gdzie ak spełnia układ równań

VkT⁢A⁢Vk⁢a=VkT⁢A⁢x*-x0=VkT⁢r0.

(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 xk jest istotnie obliczalny: do jego wyznaczenia nie jest nam efektywnie potrzebna znajomość rozwiązania!

6.1.1. Implementacja

Aby wyznaczyć xk, 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 xk.

Ponieważ A jest symetryczna, istnieje baza ortogonalna w RN złożona z wektorów własnych q1,…,qN:

A⁢qi=λi⁢qi,⁢i=1,…,N.

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

Q=q1|q2|…|qN,

mamy, że Q jest macierzą ortogonalną, QT⁢Q=I=Q⁢QT, a ponadto A ma rozkład:

A=Q⁢Λ⁢QT,

gdzie

Λ=λ1⋱λN.

Gdyby baza przestrzeni Kk była A–ortogonalna (z powodów historycznych, jej elementy oznaczymy p0,…,pk-1 tak, że Vk=p0|p1|…pk-1), tzn. piT⁢A⁢pj=0 dla i≠j, to wtedy macierz równań normalnych byłaby diagonalna,

VkT⁢A⁢Vk=⋮…piT⁢A⁢pj…⋮=p0T⁢A⁢p0⋱pk-1T⁢A⁢pk-1.

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

xk=x0+∑i=0k-1piT⁢r0piT⁢A⁢pi⁢pi.

(6.7)

Zatem potrzebna jest nam skuteczna metoda wyznaczania bazy ortogonalnej w przestrzeni Kk…. Oczywiście, ze względu na koszt obliczeniowy i pamięciowy, generowanie i następnie ortogonalizacja oryginalnego zestawu wektorów r0,A⁢r0,…,Ak-1⁢r0 rozpinających Kk 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, rk, jest prostopadłe do Kk. Ponadto rk∈Kk+1.

Dowód

Uzasadnienie pierwszej części łatwo wynika z układu równań normalnych (6.6), określającego pośrednio xk. Rzeczywiście, ponieważ xk-x0=Vk⁢ak, to z (6.6) wynika, że VkT⁢A⁢xk-x0=VkT⁢A⁢x*-x0. Upraszczając wyrazy z x0, dostajemy VkT⁢A⁢x*-xk=0, czyli VkT⁢rk=0.

Druga część wynika natychmiast z faktu, że xk∈x0+Kk, skąd A⁢xk∈A⁢x0+A⁢Kk i w konsekwencji, odejmując stronami od b, dochodzimy do rk∈r0+A⁢Kk. Tymczasem z definicji przestrzeni Kryłowa r0+A⁢Kk⊆Kk+1.

∎

Z powyższego wynika, że jeśli rk-1≠0, to Kk-1⊂Kk, a więc dopóki nie trafimy w rozwiązanie dokładne, x*, kolejne przestrzenie Kryłowa w metodzie C⁢G tworzą ściśle wstępujący ciąg przestrzeni, K0⊂K1⊂…⁢Kk-1⊂Kk⊆RN.

W dalszm ciągu założymy więc, że rk-1≠0 — a więc, że xk-1≠x*. Przypuśćmy, że mamy już zadaną bazę A–ortogonalną p0,…,pk-1 przestrzeni Kk-1 i znamy xk-1, rk-1. Naszym celem będzie wyznaczenie xk, rk oraz pk. Z zależności (6.7) mamy, że

xk=xk-1+αk⁢pk-1,

(6.8)

gdzie

αk=pk-1T⁢r0pk-1T⁢A⁢pk-1.

(6.9)

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

rk=rk-1-αk⁢A⁢pk-1.

(6.10)

Potrzeba nam jeszcze zależności rekurencyjnej pozwalającej wyznaczyć pk — ostatni wektor bazy ortogonalnej dla Kk+1 (wcześniejsze znamy z założenia indukcyjnego). Ponieważ z lematu 6.1 wynika, że Kk+1=span⁢rk,p0,…,pk-1, znaczy to, że

pk=rk+βk⁢pk-1+γk⁢pk-2+….

Mnożąc skalarnie tę równość przez A⁢pk-1 dostajemy z założenia A-ortogonalności

βk=-rkT⁢A⁢pk-1pk-1T⁢A⁢pk-1

(6.11)

i podobnie, że γk oraz wszystkie następne współczynniki są równe zero (ponieważ z lematu o residuach rk jest ortogonalne do Kk, a A⁢pj∈Kk dla j≤k-2). Zatem ostatecznie dostajemy kolejną elegancką zależność rekurencyjną, tym razem na wektory bazy A–ortogonalnej dla Kk+1:

pk=rk+βk⁢pk-1.

(6.12)

Ponieważ p0=r0, tym samym zależności (6.8)—(6.11) stanowią domknięty układ: startując z zadanego x0, 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-A⁢x;

ρ0=r22, β=0, k=1

while not stop

begin

p=r+β⁢p

w=A⁢p

α=ρk-1pT⁢w

x=x+α⁢p

r=r-α⁢w

ρk=r22

β=ρkρ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ń Kk), 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

αk=rk-122pk-1T⁢A⁢pk-1.

(6.13)

Ponieważ z założenia xj=x0+∑i=0j-1αi+1⁢pi dla każdego j≤k, to (mnożąc obustronnie tę równość przez A i odejmując stronami od b) zachodzi także rj=r0-∑i=0j-1αi+1⁢A⁢pi. Mnożąc skalarnie tę równość przez pj i uwzględniając A–ortogonalność kierunków pi dochodzimy do wniosku, że

pjT⁢rj=pjT⁢r0.

Z drugiej zaś strony, mnożąc (6.12) skalarnie przez rj otrzymujemy

pjT⁢rj=rj22-βj-1⁢pj-1T⁢rj=rj22,

ponieważ rj jest prostopadłe do Kj, w której zawarty jest wektor pj-1. Ostatecznie więc pjT⁢r0=rj22 dla każdego j≤k. Biorąc j=k-1, z (6.9) otrzymujemy (6.13).

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

rkT⁢rk=rkT⁢rk-1-αk⁢rkT⁢A⁢pk-1.

Ponieważ z lematu o ortogonalności residuów rk-1∈Kk oraz rk jest ortogonalne do Kk, to rkT⁢rk-1=0, więc podstawiając do powyższego wzoru uzyskane przed chwilą nowe wyrażenie na αk dostajemy

rk22=-rk-122pk-1T⁢A⁢pk-1⁢rkT⁢A⁢pk-1=rk-122⁢βk.

Stąd i z (6.13) już wynika wzór na współczynnik βk,

βk=rk22rk-122.

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 pk 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,

xk-xA≤2⁢⁢κ-1κ+1k⁢⁢x0-xA,

gdzie κ=cond2⁡A=λmax⁢A/λmin⁢A.

Dowód

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

pk⁢A⁢x*-x0A≤pk⁢A2⁢x*-x0A

oraz

pk⁢A2=maxλ∈σ⁢A⁡pk⁢λ,

zatem wystarczy oszacować wartości wybranego wielomianu pk∈P~k (nasza norma błędu jest i tak nie większa). Niech M=λmax⁢A i m=λmin⁢A. Jako pk w (6.3) weźmy przeskalowany k-ty wielomian Czebyszewa,

pk⁢z=Tk⁢M+m-2⁢zM-mTk⁢M+mM-m.

Rzeczywiście, pk∈P~k. Ponadto, ponieważ wielomiany Czebyszewa spełniają zależność

Tk⁢x≤1⁢⁢ dla ⁢x∈-1,1,

maxz∈m,M⁡pk⁢z≤1Tk⁢M+mM-m=1Tk⁢κ+1κ-1.

Należy więc oszacować Tk⁢κ+1κ-1. Ponieważ κ+1κ-1>1, to skorzystamy ze wzoru

Tk⁢x=12⁢x+x2-1k+12⁢x-x2-1k⁢⁢ dla ⁢x≥1.

W szczególności więc, biorąc x=κ+1κ-1 mamy

Tk⁢x≥12⁢x+x2-1k=12⁢κ+1κ-1k.

∎

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⁢λ=λ-λ1⁢⋯⁢λ-λm.

Quiz 6.1

najszybszego spadku
Tak
Nie
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
gradientów sprzężonych
Tak
Nie
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=BT⁢B+p⁢I,

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) ω=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+p⁢I dla p>0 (brak symetrii) tak dobranego, by Asym>0 oraz A=BT⁢B+p⁢I 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 TN. 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 x0+Kk, residuum (stąd nazwa) mierzone w normie euklidesowej:

rk2=b-A⁢xk2≤b-A⁢x2⁢⁢∀x∈x0+Kk.

(6.14)

Stwierdzenie 6.4 (GMRES jako metoda bezpośrednia)

Zadanie minimalizacji (6.14) ma jednoznaczne rozwiązanie. Jeśli Vk jest macierzą, której kolumny tworzą bazę Kk, to xk jest dane wzorem xk=x0+Vk⁢ak, gdzie ak jest rozwiązaniem zadania najmniejszych kwadratów względem a∈Rdim⁡Kk,

r0-A⁢Vk⁢a2=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⁢Λ⁢X-1 i Λ jest macierzą diagonalną⁵Oczywiście, na diagonali Λ znajdują się wartości własne macierzy A. W ogólności więc, zarówno X, jak i Λ mogą być zespolone.), to k-ta iteracja GMRES spełnia

rk2≤cond2⁡X⁢minp∈P~k⁡maxλ∈σ⁢A⁡p⁢λ.

Dowód

Z (6.4) wynika, że

rk2≤minp∈P~k⁡p⁢A2⁢⁢r02,

więc wystarczy oszacować p⁢A2:

p⁢A2=p⁢X⁢Λ⁢X-12=X⁢p⁢Λ⁢X-12≤X2⁢p⁢Λ2⁢X-12=cond2⁡X⁢maxλ∈σ⁢A⁡p⁢λ.

∎

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 AT⁢A=A⁢AT, to wtedy X-1=XT i w konsekwencji cond2⁡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-λ1⁢⋯⁢z-λm/λ1⁢⋯⁢λ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ć xk. Oznaczamy, jak zwykle, przez Vk macierz, której kolumny tworzą bazę Kk. Wtedy, jeśli nie osiągnięto jeszcze rozwiązania, to dim⁡Kk=k, zatem Vk∈RN×k. Aby wyznaczyć współczynniki ak reprezentacji xk w bazie Vk, xk=x0+Vk⁢ak, musimy rozwiązać względem a∈Rk liniowe zadanie najmniejszych kwadratów,

r0-A⁢Vk⁢a2=min!

(6.16)

Wydawać by się mogło, że sprawa jest tu prosta, bo A⁢Vk 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 xk potrzebna nam będzie znajomość nie samego ak, ale także Vk, sytuacja robi się nieco niezręczna. Dlatego po raz kolejny spróbujemy spreparować bazę Vk i jednocześnie wykorzystać fakt, że A⁢Vk rozpina podprzestrzeń w Kk+1 (rozpiętej przez Vk+1).

Zwróćmy uwagę na to, że jeśli v1,…,vj są bazą ortogonalną w Kj dla j=1,…,k, to aby rozszerzyć ją do Kk+1 wystarczy zortogonalizować — zamiast Ak⁢r0 — wektor A⁢vk. Stosując zmodyfikowany algorytm ortogonalizacji Grama-Schmidta do A⁢vk:

Metoda Arnoldiego wyznaczania bazy ortogonalnej

v1=r0/r02

for j=1 to k-1 begin

hi⁢j=viT⁢A⁢vj, dla i=1,…,j

vj+1=A⁢vj-∑i=1jhi⁢j⁢vi {A⁢wj to wektor, który będziemy ortogonalizować}

hi+1,j=vj+12

vj+1=vj+1/hi+1,j {unormowanie vj+1}

end

dostajemy wprost z relacji A⁢vk=∑i=1j+1hi⁢j⁢vi związek:

A⁢Vk=Vk+1⁢Hk,

(6.17)

gdzie Hk jest macierzą Hessenberga rozmiaru k+1×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

r0-A⁢Vk⁢a2=r0-Vk+1⁢Hk⁢a2=Vk+1T⁢r0-Hk⁢a2

— w ostatniej równości skorzystaliśmy z ortogonalności Vk+1. Ale przecież r0=β⁢v1 (gdzie β=r02), zaterm Vk+1T⁢r0=β,0,…,0T=β⁢e1! Tak więc, zadanie najmniejszych kwadratów (6.16) dla macierzy rozmiaru N×k sprowadziliśmy do zadania najmniejszych kwadratów:

β⁢e1-Hk⁢a2=min!

(6.18)

dla macierzy Hessenberga Hk, rozmiaru k+1×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

r0=b-A⁢x0, V1=r0/r02, k = 0

while not stop

{wykonaj ortogonalizację metodą Arnoldiego}

hk=VkT⁢A⁢vk

vk+1=A⁢vk-Vk⁢hk

hk+1,k=vk+12

vk+1=vk+1/hk+1,k

{rozszerz bazę Vk i macierz Hk}

Vk+1=Vkvk+1

Hk=Hk-1hkhk+1,k

{znajdź minimum zadania najmniejszych kwadratów}

rk+1=mina⁡β⁢e1-Hk⁢a2

{zwiększ licznik iteracji}

k=k+1

end

{rozwiąż ostatnie zadanie najmniejszych kwadratów}

ak-1=arg⁢mina⁡β⁢e1-Hk-1⁢a2

xk=x0+Vk-1⁢ak-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, rk), nie musimy znajdować ani ak, ani xk. Dzięki temu, koszt jednej iteracji GMRES jest tylko rzędu O⁢k⁢N flopów (gdybyśmy wyznaczali ak na każdej iteracji, dodawałoby to dodatkowe O⁢k2 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 Vk. Dlatego uzupełnia się go o warunkową reortogonalizację Vk: 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=α⁢IX0β⁢I,

gdzie α,β∈R∖0, a X jest macierzą n×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 Kk. 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.