Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Matematyka obliczeniowa II – 5. Proste metody iteracyjne rozwiązywania układów równań liniowych – MIM UW

Zagadnienia

5. Proste metody iteracyjne rozwiązywania układów równań liniowych

Zajmiemy się przybliżonymi metodami rozwiązywania układów równań liniowych

Ax=b

z nieosobliwą rzeczywistą macierzą A. Jak już sugerowaliśmy we wstępie, wraz z coraz większymi modelami pojawiającymi się w praktyce obliczeniowej, częściej zachodzi potrzeba rozwiązywania zadań algebry liniowej o bardzo wielkiej liczbie niewiadomych. Gdy A jest dużego rozmiaru, metody bezpośrednie (wykorzystujące na przykład rozkład LU metodą eliminacji Gaussa lub inne rozkłady macierzy, np. QR) mogą być zbyt kosztowne lub mniej wygodne w użyciu. Przykładowo, jeśli niewiadomych jest kilka milionów, a macierz nie ma żadnej specjalnej struktury, to korzystając z eliminacji Gaussa musielibyśmy wykonać rzędu 10^{{19}} działań arytmetycznych (i zarezerwować rzędu 10^{{13}} bajtów — czyli 10000 GB — pamięci na czynniki rozkładu). Zakładając optymistycznie, że nasz komputer jest w stanie wykonać 10^{{10}} działań na sekundę (czyli jest w stanie realnie osiągnąć wydajność 10 Gflopów/s), dawałoby to 10^{9} sekund, czyli około 32 lat… Rozsądną alternatywą może być wtedy rozwiązanie układu w sposób przybliżony, ale za to znacznie tańszy lub wygodniejszy. Jednym ze sposobów przybliżonego rozwiązywania układów równań liniowych są metody iteracyjne, szczególnie użyteczne wówczas, gdy A jest macierzą rozrzedzoną.

5.1. Macierze rozrzedzone

Macierze, w których bardzo wiele elementów jest zerowych, nazywamy macierzami rozrzedzonymi lub, potocznie, rzadkimi. Dla odróżnienia, macierze, które nie są rzadkie, nazwiemy gęstymi — przykładem takiej macierzy jest macierz, której wszystkie elementy są niezerowe. Intuicyjnie możemy spodziewać się, że praktyczne zadania liniowe wielkiego wymiaru będą prowadziły właśnie do macierzy rozrzedzonej, gdyż związki pomiędzy niewiadomymi w pojedynczym równaniu zadanego układu nie będą raczej dotyczyć wszystkich, tylko wybranej (nielicznej) grupy innych niewiadomych.

Wykorzystanie rozrzedzenia macierzy nie tylko może doprowadzić do algorytmów istotnie szybszych od ich analogonów dla macierzy gęstych (jak pamiętamy, standardowe algorytmy oparte na rozkładzie macierzy gęstej, np. LU, potrzebują O(N^{3}) działań arytmetycznych), ale wręcz może być jedynym realnym sposobem na to, by niektóre zadania w ogóle stały się rozwiązywalne przy obecnym stanie techniki obliczeniowej!

Gdy A jest rozrzedzona, mnożenie takiej macierzy przez wektor jest bardzo tanie (koszt jest proporcjonalny do liczby niezerowych elementów macierzy). Dlatego, konstruując metodę przybliżonego rozwiązywania układu, warto oprzeć się na iteracji, których głównym składnikiem jest operacja mnożenia przez A lub jej część.

5.1.1. Przykłady macierzy rozrzedzonych

Jak zobaczymy poniżej, rzeczywiście łatwo można spotkać realne zadania matematyki stosowanej, w których macierz wymiaru N ma tylko O(N) niezerowych elementów. Omówimy tu dwie klasy takich zadań: dyskretyzacje równań różniczkowych oraz łańcuchy Markowa z dyskretną przestrzenią stanów.

Równania różniczkowe cząstkowe

Jednym ze szczególnie ważnych źródeł układów równań z macierzami rozrzedzonymi są równania różniczkowe cząstkowe (pochodzące np. z modeli pogody, naprężeń w konstrukcji samochodu, przenikania kosmetyków do głębszych warstw skóry, itp.).

Przykład 5.1 (Dyskretyzacja jednowymiarowego laplasjanu)

Rozważmy modelowe eliptyczne równanie różniczkowe z jednorodnym warunkiem brzegowym Dirichleta

\displaystyle-u^{{\prime\prime}}(x)=f(x)\qquad\forall x\in(0,1), (5.1)
\displaystyle u(0)=u(1)=0, (5.2)

w którym u jest szukanym rozwiązaniem, a f — zadaną funkcją. Aby znaleźć jego przybliżone rozwiązanie, możemy na przykład zastąpić równanie różniczkowe odpowiadającym mu równaniem różnicowym,

\displaystyle-\dfrac{u_{{i-1}}-2u_{i}+u_{{i+1}}}{h^{2}}=f(x_{i})\qquad\forall i=1,\ldots,N, (5.3)
\displaystyle u_{0}=u_{{N+1}}=0, (5.4)

gdzie u_{i} ma odpowiadać wartości rozwiązania w węźle dyskretyzacji x_{i}, u_{i}\approx u(x_{i}), oraz x_{i}=i\cdot h, przy czym h=1/(N+1). Odpowiedzią na pytanie, jak dobrze (i czy w ogóle) takie rozwiązanie dyskretne aproksymuje rozwiązanie dokładne, zajmujemy się na wykładzie z Numerycznych równań różniczkowych; tutaj zobaczymy, do jakiego zagadnienia algebraicznego prowadzi równanie różnicowe (5.3).

Jest to oczywiście równanie liniowe na współczynniki U_{N}=(u_{1},\ldots,u_{N})^{T},

T_{N}\, U_{N}=h^{2}F_{N},

gdzie

T_{N}=\begin{pmatrix}2&-1&&&\\
-1&2&-1&&\\
&\ddots&\ddots&\ddots&\\
&&-1&2&-1\\
&&&-1&2\end{pmatrix}_{{N\times N}} (5.5)

oraz F_{N}=(f(x_{1}),\ldots,f(x_{N}))^{T}. Ponieważ interesują nas dobre przybliżenia (to znaczy przypadek, gdy h jest bardzo małe), znaczy to, że N będzie bardzo duże. Nasza macierz jest rzeczywiście macierzą rozrzedzoną, bo ma nieco mniej niż 3N niezerowych elementów. Jej współczynnik wypełnienia (ang. density), czyli stosunek liczby niezerowych do liczby wszystkich elementów macierzy, jest rzędu 3/N i maleje ze wzrostem N.

Otrzymana powyżej macierz, aczkolwiek ma bardzo wiele cech charakterystycznych dla dyskretyzacji ,,poważniejszych” równań różniczkowych, wydaje się jednak zbyt trywialna, by traktować ją metodą inną niż eliminacja Gaussa (uwzględniając jej pasmową strukturę): rzeczywiście, jest to macierz trójdiagonalna (a przy tym: symetryczna i dodatnio określona), zatem właściwie zrealizowany algorytm eliminacji Gaussa (lub lepiej: rozkładu LDL^{T}) pozwoli nam wyznaczyć rozwiązanie układu T_{N}U_{N}=h^{2}F_{N} kosztem O(N), a więc — z dokładnością do stałej — optymalnym.

Jednak, gdy przejdziemy do wyższych wymiarów i być może dodatkowo regularną siatkę węzłów dyskretyzacji zastąpimy na przykład nieregularną siatką elementu skończonego, uzyskamy macierze o znacznie bardziej skomplikowanej strukturze, które już nie tak łatwo poddają się eliminacji Gaussa.

Przykład 5.2 (Dyskretyzacja wielowymiarowego laplasjanu)

Przez analogię z poprzednim przykładem, rozważmy równanie Poissona w obszarze \Omega\subset\mathbb{R}^{d}, z jednorodnym warunkiem Dirichleta:

\displaystyle-\Delta u(x)=f(x)\qquad\forall x\in\Omega, (5.6)
\displaystyle u_{{|_{{\partial\Omega}}}}=0, (5.7)

w którym u:\Omega\rightarrow\mathbb{R} jest szukanym rozwiązaniem, a f:\Omega\rightarrow\mathbb{R} — zadaną funkcją. \Delta oznacza operator Laplace'a,

\Delta\equiv\dfrac{\partial^{2}}{\partial x_{1}^{2}}+\ldots+\dfrac{\partial^{2}}{\partial x_{d}^{2}}.

W dalszym ciągu przyjmiemy, dla ominięcia dodatkowych trudności, że \Omega jest kostką jednostkową: \Omega=(0,1)^{d}.

Aby znaleźć przybliżone rozwiązanie, znów możemy na przykład zastąpić równanie różniczkowe odpowiadającym mu równaniem różnicowym. Wybierzemy tu najmniej efektywną, ale za to koncepcyjnie najprostszą metodę dyskretyzacji, opartą na równomiernej siatce węzłów w \Omega.

Przykładowo, jeśli d=2, to weźmiemy siatkę węzłów p_{{ij}}=(x_{i},y_{j})\in\Omega, gdzie x_{i}=ih_{x},y_{j}=jh_{y}, przy czym h_{x}=a/(n_{x}+1) jest krokiem siatki w kierunku x i analogicznie h_{y}=b/(n_{y}+1), zob. rysunek LABEL:regulargrid.pstex_t.

Rys. 5.1. Regularna siatka na \Omega.

Zastępując pochodne ilorazami różnicowymi [8]

\Delta u(x_{i},y_{j})\approx\Delta _{h}u_{{ij}}=\frac{1}{h^{2}}(4u_{{i,j}}-u_{{i+1,j}}-u_{{i-1,j}}-u_{{i,j+1}}-u_{{i,j-1}}), (5.8)

dostajemy równanie różnicowe

\displaystyle-\dfrac{u_{{i-1,j}}-2u_{{i,j}}+u_{{i+1,j}}}{h_{x}^{2}}-\dfrac{u_{{i,j-1}}-2u_{{i,j}}+u_{{i,j+1}}}{h_{y}^{2}}=f_{{ij}}\qquad\forall i=1,\ldots,n_{x},\; j=1,\ldots,n_{y}. (5.9)

W pozostałych węzłach siatki rozwiązanie przyjmuje wartość równą zero na mocy warunku brzegowego Dirichleta. Więcej na temat metod dyskretyzacji takich i podobnych równań różniczkowych można dowiedzieć się na wykładzie Numeryczne równania różniczkowe.

W dalszym ciągu, dla uproszczenia zapisu będziemy zakładać, że n_{x}=n_{y}=n i w konsekwencji h_{x}=h_{y}=h. W takim układzie mamy do wyznaczenia N=n^{2} niewiadomych:

U_{N}=(\underbrace{u_{{1,1}},u_{{2,1}},u_{{3,1}},\ldots,u_{{N_{x},1}}}_{{\text{niewiadome pierwszej kolumny}}},\underbrace{u_{{1,2}},u_{{2,2}},u_{{3,2}},\ldots,u_{{N_{x},2}}}_{{\text{drugiej kolumny}}},\ldots,\underbrace{u_{{1,N_{y}}},u_{{2,N_{y}}},u_{{3,N_{y}}},\ldots,u_{{N_{x},N_{y}}}}_{{\text{ostatniej kolumny}}})^{T}.

które spełniają liniowy układ równań o postaci macierzowej

P_{N}\, U_{N}=h^{2}F_{N},

gdzie P_{N} jest tym razem pięciodiagonalną macierzą o blokowej strukturze,

P_{N}=\begin{pmatrix}T_{n}+2I&-I&&&\\
-I&T_{n}+2I&-I&&\\
&\ddots&\ddots&\ddots&\\
&&-I&T_{n}+2I&-I\\
&&&-I&T_{n}+2I\end{pmatrix}_{{N\times N}}. (5.10)

W kolejnych diagonalnych blokach P_{N} znajdują się macierze trójdiagonalne T_{n}+2I, gdzie T_{n} jest znaną już nam macierzą dyskretyzacji jednowymiarowego laplasjanu. Używając więc symbolu Kroneckera dla iloczynu tensorowego macierzy mamy przy naszej numeracji niewiadomych

P_{N}=I\otimes T_{n}+T_{n}\otimes I.

W przypadku tej macierzy, prosta strategia polegająca na eliminacji Gaussa (z wykorzystaniem dodatkowo faktu, że P_{N} jest pasmowa, o szerokości pasma 2n), dałaby nam szansę rozwiązać to równanie kosztem O(n^{4}), co jest, przy dużych wielkościach n, wartością trudną do zaakceptowania.

Prowadząc analogicznie dyskretyzację na jednorodnej siatce trójwymiarowego zadania Poissona na kostce jednostkowej (d=3), dostalibyśmy N=n^{3} oraz zadanie z macierzą siedmiodiagonalną (o szerokości pasma 2n^{2}):

S_{N}=I\otimes I\otimes T_{n}+I\otimes T_{n}\otimes I+T_{n}\otimes I\otimes I. (5.11)

Nawet dla średnich wartości n, układ ten trudno rozwiązać metodą bezpośrednią, opartą na prostym rozkładzie pasmowej macierzy (tutaj na przykład: rozkładzie Cholesky'ego), zarówno ze względu na koszt obliczeniowy, jak i koszt pamięciowy. Rzeczywiście, koszt rozkładu Cholesky'ego uwzględniającego to, że w czynnikach rozkładu poza pasmem nie pojawią się elementy niezerowe wyniesie O(n^{4}) w przypadku dwuwymiarowego równania Poissona i O(n^{7}) w przypadku trójwymiarowego równania Poissona (dlaczego?). To oczywiście mniej niż pesymistyczne O(N^{3})=O(n^{9}), ale gdy n jest duże, może być wystarczająco zachęcające do tego, by poszukać tańszych rozwiązań, dających sensowne rozwiązanie w sensownym czasie.

Przykład 5.3 (ewolucyjne równanie dyfuzji)

Rozważmy równanie

\displaystyle u_{t}(x,t)-\Delta _{x}u(x,t)=f(x,t) \displaystyle\text{dla }(x,t)\in\Omega\times(0,T),
\displaystyle u(x,t)=0 \displaystyle\text{dla }(x,t)\in\partial\Omega\times[0,T).

Dla uproszczenia przyjmiemy, że \Omega jest kwadratem jednostkowym. Wprowadzając dyskretyzację tego równania po zmiennej przestrzennej jak poprzednio, a po zmiennej czasowej — niejawny schemat Eulera (patrz rozdział o schematach różnicowych Numerycznych równań różniczkowych), otrzymujemy, po uwzględnieniu warunku brzegowego, liniowy układ równań algebraicznych na przybliżenie U_{N}^{k} rozwiązania w chwili t_{k},

\frac{U_{N}^{k}-U_{N}^{{k-1}}}{\tau}+P_{N}U_{N}^{k}=f_{N}^{k},

gdzie P_{N} jest określone jak w poprzednim przykładzie oraz

f_{N}^{k}=\begin{pmatrix}f(x_{{ij}},t_{k})\end{pmatrix}

i węzły w \Omega numerujemy zgodnie z regułą podaną powyżej. Grupując niewiadome po lewej stronie, dostajemy równanie

(I+\tau P_{N})U_{N}^{k}=U_{N}^{{k-1}}+\tau f_{N}^{k},

a więc takie, w którym macierz układu jest dalej rozrzedzona, ale — dla dostatecznie małych \tau — jest ,,małym zaburzeniem” macierzy jednostkowej.

Gdy siatka nie jest regularna1Zazwyczaj w zaawansowanych metodach dyskretyzacji równań różniczkowych stosuje się właśnie takie siatki, często — adaptacyjnie dopasowane do przebiegu rozwiązania., macierze dyskretyzacji także tracą regularną strukturę (w powyższych przykładach była ona odwzwierciedlona w tym, że niezerowe elementy macierzy układały się wzdłuż nielicznych diagonali). Macierze bez regularnej struktury są jeszcze trudniejsze dla metod bezpośrednich ze względu na rosnące kłopoty w uniknięciu wypełnienia czynników rozkładu macierzy.

Przykład 5.4 (Macierz z kolekcji Boeinga)

Spójrzmy na macierz sztywności dla modelu silnika lotniczego, wygenerowaną swego czasu w zakładach Boeinga i pochodzącą z dyskretyzacji pewnego równania różniczkowego cząstkowego metodą elementu skończonego. Przykład pochodzi z kolekcji Tima Davisa. Jest to mała macierz, wymiaru 8032 (w kolekcji spotkasz równania z milionem i więcej niewiadomych).

Opis zastępczy
Rys. 5.2. Struktura niezerowych elementów macierzy bcsstk38.

Jej współczynnik wypełnienia (to znaczy, stosunek liczby niezerowych do wszystkich elementów macierzy) wynosi jedynie 0.006, a więc macierz jest bardzo rozrzedzona: w każdym wierszu są średnio tylko 44 niezerowe elementy.

We wspomnianej kolekcji znajdzuje się jeszcze więcej macierzy rozrzedzonych, pochodzących z najrozmaitszych realnych modeli, nie tylko opartych na równaniach różniczkowych. Choć sporo z nich — tak jak tutaj przedstawiona — będzie symetrycznych, nie jest to regułą; inne układy równań ze zgromadzonego zbioru nie muszą mieć żadnej z pożądanych przez numeryka cech: regularnej/pasmowej struktury, symetrii, ani dodatniej określoności. Jedną z nich przedstawiamy poniżej.

Symulacje działania układów elektronicznych

Współczesne układy elektroniczne mogą zawierać miliony części (tranzystorów, rezystorów, kondensatorów), a ich projektowania jest procesem bardzo złożonym i kosztownym. Co więcej, w miarę postępów opracowywania, coraz trudniej o dokonanie w nich poważniejszych zmian projektu. Z tego powodu konieczna jest komputerowa symulacja ich działania. W ogólności prowadzi to do ogromnego układu równań różniczkowo–algebraicznych, ale także bada się modele (prawa Kirchhoffa!), sprowadzające się do rozwiązania układu równań liniowych wielkiego rozmiaru (N na poziomie kilku milionów).

Wizualizacja grafu macierzy Circuit5M o rozmiarze powyżej pięciu milionów

(Źródło: The University of Florida Sparse Matrix Collection)

Duże łańcuchy Markowa

Modele wielostanowych systemów kolejkowych (np. routera obsługującego wiele komputerów) także prowadzą do gigantycznych układów równań z macierzami rozrzedzonymi o specyficznej strukturze. Jeśli taki układ może z danego stanu przejść tylko do niewielu innych, macierz reprezentująca łańcuch Markowa będzie rozrzedzona.

Jednym z zadań, w którym rozważa się bodaj największy ,,naturalnie” pojawiający się łańcuch Markowa jest zadanie wyznaczania tzw. PageRank — wektora wartościującego strony internetowe w wyszukiwarce Google (zob. artykuł ,,Miara ważności” K. Diksa w Delcie 8/2008): macierz ta jest oryginalnie bardzo rozrzedzona, całkowicie nieregularna i niesymetryczna, a jej rozmiar jest równy liczbie indeksowanych stron WWW (a więc mniej więcej rzędu 10^{9}).

Macierz jest bardzo rzadka, bez wyraźnie regularnej struktury
Rys. 5.3. Struktura niezerowych elementów macierzy PageRank dla portalu http://www.mimuw.edu.pl.
Ćwiczenie 5.1

Wykaż, że zadanie T_{N}x=b można rozwiązać, korzystając z eleminacji Gaussa bez wyboru elementu głównego kosztem O(N).

Ćwiczenie 5.2

Niech macierz A będzie macierzą pasmową rozmiaru N\times N, o szerokości pasma 2k:

a_{{ij}}=0,\qquad\text{ gdy }|i-j|>k.

Wykaż, że rozkład LU (bez wyboru elementu głównego) macierzy A można wyznaczyć kosztem O(Nk^{2}) działań arytmetycznych. Jak zmieni się odpowiedź, gdy trzeba będzie dokonać wyboru elementu głównego w kolumnie?

Ćwiczenie 5.3

Zbadaj koszt rozkładu LU macierzy P_{N}, w którym wykorzystuje się jedynie informację o tym, że P_{N} jest macierzą pasmową.

Rozwiązanie: 

Mamy N=n^{2} i szerokość pasma jest równa n. Zatem na mocy poprzedniego zadania, koszt jest rzędu O(n^{4}).

5.2. Macierze rozrzedzone: implementacja

Zanim przystąpimy do omawiania metod rozwiązywania układów równań liniowych z macierzami rozrzedzonymi, warto zapoznać się ze sposobami reprezentacji (formatami) macierzy rozrzedzonych. Ponieważ macierze rozrzedzone mają dużo zerowych elementów, ogólną zasadą jest zapamiętywanie tylko tych różnych od zera, w połączeniu z informacją o ich lokalizacji w macierzy. Spośród wielu struktur danych wygodnych dla przechowywania macierzy rozrzedzonych, opisanych m.in. w monografii [13], największą popularnością wśród numeryków cieszą się dwa: format współrzędnych oraz format spakowanych kolumn (lub wierszy)2Mniej popularny format diagonalny spotkamy m.in. w reprezentacji macierzy pasmowych w LAPACKu.. Liczbę niezerowych elementów macierzy A rozmiaru N\times M będziemy oznaczali, przez analogię do funkcji MATLABa, \text{nnz}(A).

Format współrzędnych

Jest to najprostszy sposób reprezentacji macierzy rozrzedzonych. Do zapamiętania macierzy A rozmiaru N\times M, liczącej \text{nnz}(A) niezerowych elementów, wykorzystujemy trzy wektory: I, J — oba typu int — oraz V, typu double, wszystkie o długości \text{nnz}(A), przy czym zachodzi

A_{{I[k],J[k]}}=V[k],\qquad\forall k=0,\ldots,\text{nnz}(A)-1,

zob. rysunek 5.4.

Schemat rozmieszczenia elementów macierzy w wektorach
Rys. 5.4. Format współrzędnych reprezentacji macierzy rozrzedzonych.
Przykład 5.5 (Macierz jednowymiarowego laplasjanu w formacie współrzędnych)

Zdefiniujemy niezbędne struktury danych i następnie wypełnimy je tak, by reprezentować macierz T_{N} jednowymiarowego laplasjanu (5.5) w formacie współrzędnych. Macierz T_{N} jest kwadratowa rozmiaru N\times N, o co najwyżej 3N niezerowych elementach.

Poniższa struktura:

typedef struct {
int I;
int J;
double V;
} SpElem;
będzie odpowiadała jednemu elementowi reprezentowanej macierzy. Cała macierz T_{N} zmieści się zatem w tablicy 3N elementów typu SpElem, na które przydzielamy miejsce:

SpElem *T;

/* alokujemy miejsce na 3N elementów */
T = (SpElem *) calloc(3*N, sizeof(SpElem));

Zwróćmy uwagę na to, że do przydzielenia pamięci użyliśmy tym razem funkcji calloc(), a nie malloc(), jak zwykle. To dla większego bezpieczeństwa, bo calloc() wypełnia zerami przydzieloną pamięć.

Teraz możemy wypełnić tablicę. Jak widzimy z dalszego ciągu kodu, zmienną T_NNZ inkrementujemy od początkowej wartości 0 za każdym razem, gdy do macierzy T zapisujemy kolejny element — w ten sposób będziemy dokładnie wiedzieli, ile w rzeczywistości (potencjalnie) niezerowych elementów wpisaliśmy do zaalokowanej tablicy.

int nnz = 0; /* bieżąca liczba wpisanych do tablicy elementów */

/* wypełniamy macierz */
for(i = 1; i <= N-1; i++)
{
	/* element diagonalny */
	T[nnz].I = i;
	T[nnz].J = i;
	T[nnz].V = 2;
	nnz++;
	/* element naddiagonalny w i-tym wierszu */
	T[nnz].I = i;
	T[nnz].J = i+1;
	T[nnz].V = -1;
	nnz++;
	/* element poddiagonalny w i-tej kolumnie */
	T[nnz].I = i+1;
	T[nnz].J = i;
	T[nnz].V = -1;
	nnz++;
}
/* ostatni element, (N,N) */
T[nnz].I = N;
T[nnz].J = N;
T[nnz].V = 2;
nnz++;

Nie zawsze taka forma implementacji formatu współrzędnych będzie możliwa: często wybór implementacji będzie narzucony przez biblioteki, jakich zechcemy używać.

Format spakowanych kolumn (lub wierszy)

Format współrzędnych nie narzucał żadnego uporządkowania elementów macierzy — można było je umieszczać w dowolnej kolejności. Z drugiej strony, narzucenie sensownego porządku mogłoby wspomóc realizację wybranych istotnych operacji na macierzy, na przykład, aby wygodnie było realizować działanie (prawostronnego) mnożenia macierzy przez wektor, wygodnie byłoby przechowywać elementy macierzy wierszami. Tak właśnie jest zorganizowany format spakowanych wierszy (CSR, ang. Compressed Sparse Row). Analogicznie jest zdefiniowany format spakowanych kolumn (CSC, Compressed Sparse Column), którym zajmiemy się bliżej.

Podobnie jak w przypadku formatu współrzędnych, macierz w formacie CSC jest przechowywana w postaci trzech wektorów. V jest wektorem typu double o długości \text{nnz}(A), zawierającym kolejne niezerowe elementy macierzy A wpisywane kolumnami, I jest wektorem typu int, także o długości \text{nnz}(A), zawierającym numery wierszy macierzy A, odpowiadających elementom z V. Natomiast zamiast tablicy J, jak to było w formacie współrzędnych, mamy krótszy wektor typu int, P, o długości M+1, zawierający na (j-1) miejscu indeks pozycji w V, od której w V rozpoczynają się elementy j-tej kolumny macierzy A (zob. rysunek 5.5).

Schemat rozmieszczenia elementów macierzy w wektorach
Rys. 5.5. Format CSC reprezentacji macierzy rozrzedzonych.

W ten sposób, wartości kolejnych elementów j-tej kolumny macierzy A znajdują się w tablicy V na pozycjach o indeksach P[j-1]P[j]-1. Konsekwentnie, zawsze P[0]=0, a ostatni element, P[M], jest równy \text{nnz}(A). Dzięki temu, takie operacje, jak na przykład wydobycie jednej kolumny z macierzy A, lub mnożenie przez wektor, są łatwiejsze do implementacji w schludnej, bezwariantowej pętli [13].

Ostatecznie mamy więc zależność

A_{{I[k],j}}=V[k],\qquad j=1,\ldots,M,\, k=P[j-1],\ldots,P[j]-1.

Z tego rodzaju formatu reprezentacji macierzy rzadkich (z drobnymi modyfikacjami) korzystają np. pakiety Octave, MATLAB i UMFPACK.

Przykład 5.6 (Procedura mnożenia macierzy w formacie CSC przez wektor)

Mając macierz w formacie CSC, zadaną tablicami I[NNZ], P[M+1], V[NNZ], możemy napisać zgrabną procedurę mnożenia jej przez wektor x[M]. Wynik będziemy zapisywać w wektorze y[N]:

for(i = 0; i < N; i++)
	y[i] = 0.0; /* zerujemy wektor wynikowy */

for(j = 0; j < M; j++)
	for(k = P[j]; k < P[j+1]; k++)
		y[I[k]] += V[k]*x[j];
Ćwiczenie 5.4

Zapisz w pseudokodzie procedurę mnożenia macierzy rozrzedzonej przez wektor, gdy jest ona reprezentowana

  1. w formacie współrzędnych,

  2. w formacie CSR.

5.3. Metody stacjonarne rozwiązywania układów równań liniowych

Teraz zajmiemy się prostymi i historycznie najstarszymi metodami iteracyjnego rozwiązywania układu równań

Ax=b,

gdzie nieosobliwa macierz A jest (zapewne wielkiego) rozmiaru N. Są one bodaj najprostsze w analizie i implementacji, ale — jak można się domyślić — w praktyce najmniej efektywne. Z drugiej strony, stanowią one ważny składnik jednej z najszybszych metod rozwiązywania niektórych trudnych układów równań (por. rozdział 8.3).

Rozważane przez nas w tym rozdziale metody opierają się na rozkładzie macierzy A na część ,,łatwo odwracalną”, M, i ,,resztę”, Z=M-A. Dokładniej, równanie Ax=b można zapisać jako zadanie punktu stałego

Mx=Zx+b,

a jeśli M jest nieosobliwa, to równoważnie:

x=M^{{-1}}(Zx+b).

Ponieważ jest to zadanie znajdowania punktu stałego postaci x=\Phi(x), można spróbować zastosować doń metodę iteracji prostej Banacha:

x_{{k+1}}=M^{{-1}}(Zx_{k}+b). (5.12)

Takie metody nazywamy stacjonarnymi metodami iteracyjnymi.

Ćwiczenie 5.5

Wykaż, że każdą stacjonarną metodę iteracyjną można zapisać w wygodniejszej (dlaczego?) do praktycznej implementacji postaci

x_{{k+1}}=x_{k}+M^{{-1}}r_{k},

gdzie r_{k}=b-Ax_{k}.

Rozwiązanie: 
x_{{k+1}}=M^{{-1}}(Zx_{k}+b)=M^{{-1}}((M-A)x_{k}+b)=x_{k}+M^{{-1}}(b-Ax_{k}).

Jest to wygodniejsza postać, bo wymaga operowania wektorem residuum r_{k}, na którym zapewne będziemy opierać kryterium stopu iteracji (więc i tak wtedy byśmy go wyznaczali)). Nie wymaga też jawnego Z.

Metody stacjonarne można zapisać w ogólnej postaci

x_{{k+1}}\,=\, Bx_{{k}}\,+\, c, (5.13)

gdzie macierz B oraz wektor c są tak dobrane, by punkt stały x^{*} równania (5.13) był rozwiązaniem równania wyjściowego Ax=b.

Dla stacjonarnej metody iteracyjnej, B=M^{{-1}}Z=I-M^{{-1}}A oraz c=M^{{-1}}b. Macierz B nazywamy macierzą iteracji, gdyż zachodzi

x_{{k+1}}-x^{*}\,=\, B(x_{{k}}-x^{*})=\ldots=B^{{k+1}}(x_{{0}}-x^{*}).

Jasne jest, że z twierdzenia Banacha o punkcie stałym wynika od razu następujący warunek zbieżności:

Stwierdzenie 5.1 (warunek wystarczający zbieżności metody stacjonarnej)

Jeśli \| B\|<1, to metoda (5.13) jest zbieżna co najmniej liniowo do x^{*} dla dowolnego x_{0}\in R^{N}. Tutaj \|\cdot\| oznacza dowolną normę macierzową indukowaną przez normę wektorową.

Dowód

Przy naszych założeniach, mamy do czynienia z kontrakcją \Phi(x)=Bx\,+\, c, która odwzorowuje kulę K=\{ x\in R^{N}:\| x-x^{*}\|\leq\| x_{0}-x^{*}\|\} w siebie:

\|\Phi(x)-x^{*}\|=\|\Phi(x)-\Phi(x^{*})\|=\| Bx-Bx^{*}\|\leq\| B\|\| x-x^{*}\|<\| x-x^{*}\|.

Ponieważ, z założenia, stała Lipschitza dla \Phi wynosi \| B\|<1, to z twierdzenia Banacha o kontrakcji (por. twierdzenie 9.1) wynika teza stwierdzenia.

Wniosek 5.1

Jeśli \| I-M^{{-1}}A\|<1, to ciąg określony przez stacjonarną metodę iteracyjną jest zbieżny do x^{*}=A^{{-1}}b oraz

\| x_{{k+1}}-x^{*}\|\leq\| I-M^{{-1}}A\|\| x_{{k}}-x^{*}\|.
Dowód

Aby stwierdzić zbieżność, wystarczy w poprzednim stwierdzeniu podstawić B=M^{{-1}}(M-A)=I-M^{{-1}}A oraz x=x_{k}. Oszacowanie błędu jest konsekwencją własności kontrakcji.

Warunek konieczny i dostateczny zbieżności tej iteracji dla dowolnego wektora startowego x_{0} podaje poniższe twierdzenie. Podkreślmy, że chodzi nam tu o zbieżność metody do rozwiązania układu, gdy startujemy z dowolnego przybliżenia początkowego x_{0}.

Twierdzenie 5.1 (warunek konieczny i dostateczny zbieżności metody stacjonarnej)

Niech B=I-M^{{-1}}A będzie macierzą iteracji metody (5.13). Ciąg x_{k} określony tą metodą jest zbieżny do rozwiązania x^{*} dla dowolnego x_{0}\in R^{N} wtedy i tylko wtedy, gdy

\rho(B)=\max\{|\lambda|:\lambda\text{ jest wartością własną }B\}<1.
Dowód

Na mocy (5.13) otrzymujemy równanie błędu e_{k}=x_{k}-x^{*},

e_{{k}}=Be_{{k-1}}=\ldots=B^{k}e_{0}.

Konieczność warunku jest oczywista: jeśli metoda jest zbieżna dla każdego x_{0}, to w szczególności dla takiego, że e_{0} jest równe wektorowi własnemu macierzy B odpowiadającego zadanej (dowolnej) wartości własnej \lambda. W konsekwencji,

e_{k}=B^{k}e_{0}=\lambda^{k}e_{0}\rightarrow 0,

a to oznacza, że |\lambda|<1.

To, że warunek wymieniany w tezie twierdzenia jest także wystarczający dla zbieżności metody, wynika z faktu (zob. [11]), że dla dowolnego \epsilon>0 istnieje norma macierzowa indukowana przez normę wektorową \|\cdot\| _{\epsilon} taka, że

\| B\| _{\epsilon}\leq\rho{B}+\epsilon.

Skoro więc z założenia \rho(B)<1, to można dobrać takie \epsilon, by \| B\| _{\epsilon}<1 i w konsekwencji

\| e_{k}\| _{\epsilon}\leq\| B\| _{\epsilon}^{k}\| e_{0}\| _{\epsilon},

skąd x_{k}\rightarrow x^{*}.

Ćwiczenie 5.6

W powyższym dowodzie, jeśli \lambda jest zespolona, to odpowiadający jej wektor własny też jest zespolony, a więc — formalnie — nie może być uznany za kontrprzykład zbieżności, bo ta ma być ,,dla każdego x_{0}\in\mathbb{R}^{N}”, a nie — ,,dla każdego x_{0}\in\mathbb{C}^{N}”. Uzupełnij tę lukę pokazując, że jeśli nawet jeśli \lambda jest zespolona, to gdy |\lambda|\geq 1 można wskazać taki rzeczywisty wektor x_{0}, że zbieżność nie będzie zachodzić.

Zaletą stacjonarnych metod iteracyjnych jest ich prostota powodująca, że są one wyjątkow wdzięczne do szybkiego zaprogramowania. Zobaczymy to na przykładzie kilku klasycznych metod stacjonarnych: Jacobiego, Gaussa–Seidela i SOR. Wszystkie będą bazować na podziale macierzy A na trzy części: diagonalną, ściśle dolną trójkątną i ściśle górną trójkątną:

\displaystyle A=L+D+U, (5.14)

gdzie

L=\begin{pmatrix}0&&&&\\
a_{{21}}&0&&&\\
a_{{31}}&a_{{32}}&0&&\\
\vdots&\vdots&\ddots&0&\\
a_{{N1}}&a_{{N2}}&\cdots&\cdots&0\end{pmatrix},\quad D=\begin{pmatrix}a_{{11}}&&&&\\
&a_{{22}}&&&\\
&&\ddots&&\\
&&&\ddots&\\
&&&&a_{{NN}}\end{pmatrix},\quad U=\begin{pmatrix}0&a_{{12}}&a_{{13}}&\cdots&a_{{1N}}\\
&0&a_{{23}}&\cdots&a_{{2N}}\\
&&0&\ddots&\vdots\\
&&&0&\vdots\\
&&&&0\end{pmatrix}

5.3.1. Metoda Jacobiego

Biorąc w (5.12) M=D, gdzie D jest macierzą diagonalną składającą się z wyrazów stojących na głównej przekątnej macierzy A (zob. (5.14)), otrzymujemy (o ile na przekątnej macierzy A nie mamy zera) metodę iteracyjną

x_{{k+1}}\,=D^{{-1}}(b-(L+U)x_{{k}}),

zwaną metodą Jacobiego.

Rozpisując ją po współrzędnych, dostajemy układ rozszczepionych równań (numer iteracji wyjątkowo zaznaczamy w postaci górnego indeksu):

x^{{(k+1)}}_{i}=\frac{1}{a_{{ii}}}\left(b_{i}-\sum _{{j\neq i}}a_{{ij}}x^{{(k)}}_{j}\right),

co znaczy dokładnie tyle, że w i-tym równaniu wyjściowego układu przyjmujemy za współrzędne x wartości z poprzedniej iteracji i na tej podstawie wyznaczamy wartość x_{i}.

Widzimy więc, że metoda rzeczywiście jest banalna w implementacji, a dodatkowo jest w pełni równoległa: każdą współrzędną nowego przybliżenia możemy wyznaczyć niezależnie od pozostałych.

Twierdzenie 5.2 (O zbieżności metody Jacobiego)

W metodzie Jacobiego warunek dostateczny zbieżności, \| D^{{-1}}(L+U)\| _{\infty}<1, jest spełniony np. wtedy, gdy macierz A ma dominującą przekątną, tzn. gdy

|a_{{ii}}|\,>\,\sum _{{j\neq i}}|a_{{ij}}|,\qquad\forall i=1,\ldots,N. (5.15)
Dowód

Rzeczywiście, ponieważ wyraz (i,j) macierzy M^{{-1}}Z wynosi 0 dla i=j oraz a_{{ij}}/a_{{ii}} dla i\neq j, to

\| M^{{-1}}Z\| _{\infty}=\max _{{1\le i\le N}}\sum _{{j=1,j\ne i}}^{N}{|a_{{ij}}|}/{|a_{{ii}}|}=\max _{{1\le i\le N}}\sum _{{j=1}}^{N}{|a_{{ij}}|}/{|a_{{ii}}|}\,-\, 1\;<\; 1,

przy czym ostatnia nierówność wynika z warunku diagonalnej dominacji.

Niestety, w wielu ważnych wypadkach metoda Jacobiego, choć zbieżna, będzie zbieżna zbyt wolno, by nas zadowolić.

Przykład 5.7 (Macierz laplasjanu)

Macierz N\times N, zwana macierzą jednowymiarowego laplasjanu

T_{N}=\begin{pmatrix}2&-1&&\\
-1&2&\ddots&\\
&\ddots&\ddots&-1\\
&&-1&2\end{pmatrix}

pojawia się w bardzo wielu zastosowaniach (patrz przykład 5.1), również jako sztucznie wprowadzane podzadanie w niektórych algorytmach numerycznych. Ta macierz jest macierzą taśmową, symetryczną i dodatnio określoną, więc układ równań z tą macierzą można bez trudu rozwiązać metodami bezpośrednimi, kosztem O(N). Jak jednak za chwilę się przekonamy, układ równań z macierzą T_{N} będzie wyjątkowo trudny dla klasycznej metody stacjonarnej.

Stosując do T_{N} metodę Jacobiego mamy M=2I oraz Z=T_{N}-M. Obliczając normę macierzy iteracji Jacobiego dostajemy \| M^{{-1}}Z\| _{\infty}=1, co — jak wiemy z twierdzenia 5.1 — nie rozstrzyga jeszcze o jej zbieżności lub niezbieżności, ale już może stanowić dla nas poważne ostrzeżenie.

Potrzebna będzie bardziej subtelna analiza. Okazuje się, że są znane wzory na wartości własne \lambda _{j} macierzy T_{N} (por. ćwiczenie 5.8):

\lambda _{j}=2\left(1-\cos\left(\frac{j\pi}{N+1}\right)\right)=4\sin^{2}\left(\frac{j\pi}{2(N+1)}\right),

(a więc, 0<\lambda _{j}<4) dla 1\leq j\leq N. W konsekwencji, wartościami własnymi M^{{-1}}Z=\frac{1}{2}Z=\frac{1}{2}T_{N}-I są liczby \mu _{i}=\frac{1}{2}\lambda _{i}-1. Ponieważ -1<\mu _{i}<1, znaczy to, że metoda Jacobiego jest zbieżna dla macierzy T_{N}.

Z drugiej strony, nie dajmy się zwieść optymizmowi matematyka (,,nie martw się, jest zbieżny…”): nietrudno sprawdzić, że \rho(M^{{-1}}Z)=\cos(\pi/(N+1))=1-\dfrac{\pi^{2}}{2(N+1)^{2}}+O(N^{{-4}}), co oznacza, że metoda Jacobiego — choć zbieżna! — dla dużych N staje się zbieżna tak wolno, że w praktyce bezużyteczna: na przykład, gdy N=100, \rho(M^{{-1}}Z)\approx 0.9995, zatem potrzebowalibyśmy około 4600 iteracji, by mieć pewność, że początkowy błąd zredukowaliśmy zaledwie… 10 razy.

Dopiero w rozdziale 8.3 przekonamy się, jak — przez niebanalne wykorzystanie głębszych informacji o samym zadaniu — można wykorzystać metodę stacjonarną do konstrukcji metody iteracyjnej o optymalnym koszcie, którą będziemy mogli stosować także do macierzy dwu- i trójwymiarowego laplasjanu.

Ćwiczenie 5.7

Wykaż, że wartości własne macierzy dwuwymiarowego laplasjanu P_{N} — zob. (5.10) — są postaci

\lambda _{{jk}}=4\left(\sin^{2}\left(\frac{j\pi}{2(n_{x}+1)}\right)+\sin^{2}\left(\frac{k\pi}{2(n_{y}+1)}\right)\right),

dla j=1,\ldots,n_{x},k=1,\ldots,n_{y}. Podobnie pokaż, że, wartości własne macierzy trójwymiarowego laplasjanu S_{N} — zob. (5.11) — są postaci

\lambda _{{jkl}}=4\left(\sin^{2}\left(\frac{j\pi}{2(n_{x}+1)}\right)+\sin^{2}\left(\frac{k\pi}{2(n_{y}+1)}\right)+\sin^{2}\left(\frac{l\pi}{2(n_{z}+1)}\right)\right),

dla j=1,\ldots,n_{x},k=1,\ldots,n_{y},l=1,\ldots,n_{z}.

Rozwiązanie: 

Niech T_{N}=V\Lambda V^{T} będzie rozkładem macierzy jednowymiarowego laplasjanu takim, że \Lambda jest macierzą diagonalną zawierającą wartości własne, a V taka, że V^{T}V=I zawiera wektory własne T_{N}. Ponieważ macierz P_{N}=I\otimes T_{N}+T_{N}\otimes I to konsekwentnie

P_{N}=\text{\TODO{dopisac z demmela}}....

Stąd wynika teza. Analogicznie postępujemy dla macierzy trójwymiarowego laplasjanu.

Ćwiczenie 5.8

Wykaż, że jeśli b\neq 0 oraz c\neq 0, to rzeczywista trójdiagonalna macierz Toeplitza3Macierzą Toeplitza nazywamy macierz o stałych współczynnikach na diagonalach.

T=\begin{pmatrix}d&b&&&\\
c&d&b&&\\
&\ddots&\ddots&\ddots&\\
&&c&d&b\\
&&&c&d\end{pmatrix}_{{N\times N}} (5.16)

ma jednokrotne wartości własne, równe

\lambda _{k}=d+2c\alpha\cos(k\theta),\qquad k=1,\ldots,N, (5.17)

gdzie

\alpha=\sqrt{\dfrac{b}{c}},\qquad\theta=\dfrac{\pi}{N+1}.

Zwróć uwagę na to, że gdy bc<0, wartości własne są zespolone. Wektor własny v_{k} macierzy T odpowiadający \lambda _{k} ma współrzędne dane wzorem

(v_{k})_{j}=\alpha^{{-j}}\sin(jk\theta),\qquad j=1,\ldots,N. (5.18)
Wskazówka: 

Dla X=\operatorname{diag}(\alpha,\alpha^{2},\ldots,\alpha^{N}) zachodzi

T=dI+c\alpha X^{{-1}}T^{0}_{N}X=X^{{-1}}\left(dI+c\alpha T^{0}\right)X,

gdzie

T^{0}=\begin{pmatrix}0&1&&&\\
1&0&1&&\\
&\ddots&\ddots&\ddots&\\
&&1&0&1\\
&&&1&0\end{pmatrix}_{{N\times N}}.

Znaczy to, że T jest podobna do macierzy dI+c\alpha T^{0} — a więc ma te same wartości własne. Nietrudno sprawdzić wprost (ze wzorów trygonometrycznych), że T^{0} ma wektory i wartości własne takie, jak przewiduje wzór.

Przykład 5.8

Niech B będzie losową kwadratową macierzą N\times N, która ma średnio 10 elementów w wierszu (zatem jej współczynnik wypełnienia wynosi 10/N) i niech

A=B+pI.

Dobierając do B wartość p możemy sterować stopniem diagonalnej dominacji macierzy A i, ostatecznie, zagwarantować sobie, że A jest diagonalnie dominująca.

W poniższym eksperymencie komputerowym możemy naocznie przekonać się, jak stopień dominacji diagonali wpływa na szybkość zbieżności metody Jacobiego. Dla porównania zobaczymy, jak z naszą macierzą radzi sobie metoda bezpośrednia: oparta na rozkładzie LU macierzy A i starająca się wykorzystać w inteligentny sposób jej rozrzedzenie.



Ćwiczenie 5.9
  • Przekonaj się, czy stopień diagonalnej dominacji (wartość parametru p w skrypcie) ma istotny wpływ na szybkość zbieżności.

  • Domyślnie testy prowadzimy dla macierzy rozmiaru 2500. Sprawdź, jak zmieni się czas wykonania skryptu, gdy N będzie jeszcze większe: która metoda: iteracyjna, czy bezpośrednia zagwarantuje nam sensowny wynik w krótszym czasie?

  • Przekonaj się, że gdy rozmiar macierzy jest niewielki, na przykład N=100, nie ma większego sensu stosować metody iteracyjnej. Przy jakim poziomie tolerancji tol metoda iteracyjna staje się konkurencyjna dla bezpośredniej (pod względem czasu działania)?

Ćwiczenie 5.10
  • Uzupełnij powyższy skrypt o wyznaczenie teoretycznej wartości współczynnika redukcji błędu w normie \|\cdot\| _{\infty} i porównaj ją z faktyczną szybkością redukcji błędu.

  • Sprawdź, że symetria A nie ma większego wpływu na szybkość zbieżności.

Przykład 5.9

Zdarzają się macierze — niestety, nie są to sztucznie generowane, akademickie przypadki — które są patologicznie trudne dla stacjonarnej metody iteracyjnej, a bardzo łatwe dla metody bezpośredniej. Jedną taką macierz już znamy: jest to macierz jednowymiarowego laplasjanu T_{N}.

W naszym teście, dla zadania z macierzą A=T_{N}+pI i różnych wartości parametru p\geq 0 zbadamy szybkość zbieżności metody Jacobiego i porównamy czas jej działania z metodą bezpośrednią (w przypadku macierzy trójdiagonalnej jest ona praktycznie optymalna).



Tym razem zmiana p wyraźnie wpływa na szybkość zbieżności, a dla p=0 zbieżność jest — jak już zdążyliśmy się przekonać także analitycznie — jest bardzo wolna.

Ćwiczenie 5.11
  • Wyjaśnij przyczynę szybkiej zbieżności metody dla dużych p.

  • Domyślnie eksperymenty prowadzimy dla macierzy rozmiaru 25. Sprawdź, jak zmieni się czas wykonania skryptu, gdy N będzie większe: która metoda zagwarantuje nam wynik w krótszym czasie? Jak odpowiedź zależy od p?

  • Jeśli zależałoby nam na obliczeniu dokładnego rozwiązania układu, powinniśmy znacząco zmniejszyć warunek stopu metody, tol, na przykład do poziomu 10^{{-16}}. Jak teraz wypada porównanie metody Jacobiego i bezpośredniej, gdy p=0?

Ćwiczenie 5.12

Rozważmy macierz

A=\begin{pmatrix}1&\alpha&\alpha\\
\alpha&1&\alpha\\
\alpha&\alpha&1\end{pmatrix}.

Wykaż, że

  • A>0\iff-1<2\alpha<2,

  • metoda Jacobiego dla tej macierzy jest zbieżna, wtedy i tylko wtedy, gdy -1<2\alpha<1.

Czy rezultat o zbieżności metody Jacobiego można uogólnić na przypadek, gdy A rozmiaru N\times N ma na diagonali same jedynki, a poza nią wszystkie jej elementy są równe \alpha?

5.3.2. Metoda Gaussa–Seidela

Heurystyka tej metody opiera się na zmodyfikowaniu metody Jacobiego tak, by w każdym momencie iteracji korzystać z najbardziej ,,aktualnych” współrzędnych przybliżenia rozwiązania x. Rzeczywiście, przecież wykonując jeden krok metody Jacobiego, czyli rozwiązując kolejno równania skalarne względem x^{{(k+1)}}_{i} dla i=1,\ldots,N:

x^{{(k+1)}}_{i}=\frac{1}{a_{{ii}}}\left(b_{i}-\sum _{{j\neq i}}a_{{ij}}x^{{(k)}}_{j}\right),

nietrudno zauważyć, że w części sumy, dla j<i, moglibyśmy odwoływać się — zamiast do ,,starych” x^{{(k)}}_{j} — do ,,dokładniejszych”, świeżo wyznaczonych, wartości x^{{(k+1)}}_{j}, tzn. ostatecznie wyznaczać

x^{{(k+1)}}_{i}=\frac{1}{a_{{ii}}}\left(b_{i}-\sum _{{j<i}}a_{{ij}}x^{{(k+1)}}_{j}-\sum _{{j>i}}a_{{ij}}x^{{(k)}}_{j}\right).

W języku rozkładu macierzy A=M-Z i iteracji x_{{k+1}}=M^{{-1}}(Zx_{k}+b) mielibyśmy więc M=L+D (dolny trójkąt macierzy A z diagonalą) oraz Z=-U (ściśle górny trójkąt A) i konsekwentnie zapis macierzowy iteracji

x_{{k+1}}=(L+D)^{{-1}}(b-Ux_{{k}}).
Twierdzenie 5.3 (O zbieżności metody Gaussa–Seidela)

Jeśli macierz A jest diagonalnie dominująca, to metoda Gaussa–Seidela jest zbieżna do x^{*} dla dowolnego wektora startowego x_{0}.

Inny wariant tej metody dostalibyśmy, biorąc za M górny trójkąt macierzy A.

Uwaga 5.1

Obie metody, Jacobiego i (zwłaszcza) Gaussa–Seidela stosuje się także czasem w prostych algorytmach rozwiązywania układów równań nieliniowych: ich zaletą jest to, że głównym składnikiem iteracji jest rozwiązywanie skalarnego równania nieliniowego na każdym kroku metody.

Metoda Gaussa–Seidela jest w wielu przypadkach rzeczywiście szybciej zbieżna od metody Jacobiego, np. tak jest w przypadku macierzy jednowymiarowego laplasjanu, patrz wniosek 5.2. Wciąż jednak, dodajmy, dla tego zadania jej zbieżność jest zbyt wolna, by ją stosować jako samodzielną metodę.

Przykład 5.10

Kontynuując przykład 5.8, porównamy szybkość zbieżności metody Gaussa–Seidela i Jacobiego na tym samym zadaniu.



Ćwiczenie 5.13

Porównaj szybkość zbieżności metod Jacobiego i Gaussa–Seidela na macierzy T_{N} jednowymiarowego laplasjanu.

Wskazówka: 

Możesz przeprowadzić porównanie, prowadząc dobrze zaplanowane eksperymenty numeryczne. W MATLABie łatwo utworzysz macierz T_{N} poleceniem

e = ones(N,1);
TN = spdiags([-e, 2*e, -e], [-1,0,1], N, N);
Jeśli interesuje Cię wynik teoretyczny, czytaj dalej.

5.3.3. Metoda SOR

Zbieżność metody Gaussa–Seidela można przyspieszyć, wprowadzając parametr relaksacji \omega i kolejne współrzędne nowego przybliżenia x_{{k+1}} wyznaczać, kombinując ze sobą poprzednie przybliżenie x^{{(k)}}_{i} oraz współrzędną nowego przybliżenia \tilde{x}^{{k+1}}_{i}, uzyskanego metodą Gaussa–Seidela:

x^{{(k+1)}}_{i}=(1-\omega)x^{{(k)}}_{i}+\omega\tilde{x}^{{k+1}}_{i}.

Gdyby \omega=1, dostalibyśmy z powrotem metodę Gaussa–Seidela. Ponieważ zwykle wybiera się \omega>1, powstałą metodę nazywa się metodą nadrelaksacji (ang. successive overrelaxation, SOR), która jest na swój sposób szczytowym osiągnięciem wśród metod stacjonarnych. Rozkład macierzy A=M-Z odpowiadający metodzie SOR z parametrem \omega jest zadany przez

M=\frac{1}{\omega}D+L,\qquad Z=\left(\frac{1}{\omega}-1\right)D-U.
Twierdzenie 5.4 (lemat Kahana o dopuszczalnych wartościach parametru relaksacji SOR)

Jeśli A ma niezerową diagonalę, a metoda SOR jest zbieżna, to musi być 0<\omega<2.

Dowód

Wystarczy zbadać promień spektralny macierzy iteracji metody SOR, która jest równa M^{{-1}}Z=(D+\omega L)^{{-1}}((1-\omega)D-\omega U). Ponieważ macierz (D+\omega L) jest macierzą trójkątną o diagonali D, to \det(D+\omega L)^{{-1}}=\det(D)^{{-1}} i w konsekwencji

\det(M^{{-1}}Z)=\det(D+\omega L)^{{-1}}\det((1-\omega)D-\omega U)=\det{D}^{{-1}}\det((1-\omega)D)=\det((1-\omega)I)=(1-\omega)^{N}.

Ze względu na to, że wyznacznik macierzy jest równy produktowi jej wartości własnych, dochodzimy do wniosku, że \rho(M^{{-1}}Z)\geq|1-\omega|, co kończy dowód.

Jako ciekawostkę odnotujmy poniższe

Twierdzenie 5.5 (Ostrowskiego i Reicha)

Jeśli A=A^{T}>0, to SOR jest zbieżna dla dowolnego \omega\in(0,2).

Dowód tego twierdzenia można znaleźć na przykład w [15, rozdział 8.3].

Dla klasy macierzy obejmującej niektóre spotykane w praktycznych zastosowaniach, można wskazać brzemienny w konsekwencje związek promienia spektralnego macierzy iteracji metody SOR i metody Jacobiego.

Definicja 5.1

Rzeczywistą macierz kwadratową A rozmiaru N, o niezerowych elementach na diagonali, nazywamy macierzą zgodnie uporządkowaną, jeśli wartości własne macierzy

J_{{\alpha}}=\alpha\tilde{L}+\dfrac{1}{\alpha}\tilde{U},

gdzie \tilde{L}=D^{{-1}}L, \tilde{U}=D^{{-1}}U, są niezależne od \alpha\in\mathbb{C}\setminus\{ 0\}.

Zauważmy, że J_{{-1}}=-\tilde{L}-\tilde{U} jest macierzą iteracji metody Jacobiego dla macierzy A.

Szczególnym przypadkiem macierzy zgodnie uporządkowanych są macierze postaci

\displaystyle A=\begin{pmatrix}D_{1}&U_{2}\\
L_{1}&D_{2}\end{pmatrix}, (5.19)

w których D_{1},D_{2} są nieosobliwymi macierzami diagonalnymi.

Rzeczywiście, dla macierzy postaci (5.19) mamy

A=L+D+U,

gdzie

L=\begin{pmatrix}0&0\\
L_{1}&0\end{pmatrix},\quad D=\begin{pmatrix}D_{1}&\\
&D_{2}\end{pmatrix},\quad U=\begin{pmatrix}0&U_{2}\\
0&0\end{pmatrix}

i w konsekwencji

J_{{\alpha}}=\begin{pmatrix}&\dfrac{1}{\alpha}D_{1}^{{-1}}U_{2}\\
\alpha D_{2}^{{-1}}L_{1}&\end{pmatrix}=\begin{pmatrix}I&\\
&\alpha I\end{pmatrix}\underbrace{\begin{pmatrix}&D_{1}^{{-1}}U_{2}\\
D_{2}^{{-1}}L_{1}&\end{pmatrix}}_{{J_{1}}}\begin{pmatrix}I&\\
&\alpha I\end{pmatrix}^{{-1}}.

Zatem macierz J_{\alpha} jest podobna do J_{1} i w konsekwencji ma te same wartości własne, co J_{1} (niezależnie od \alpha\neq 0), a więc — jest zgodnie uporządkowana.

Ćwiczenie 5.14

Wykaż, że jeśli A jest macierzą trójdiagonalną (taką, jak na przykład T_{N}) lub blokowo trójdiagonalną (taką, jak na przykład P_{N}) to A można poddać takiej permutacji P wierszy i kolumn, że uzyskana macierz PAP^{T} jest postaci (5.19)4W rzeczywistości, można wykazać więcej, że nawet bez tej permutacji macierze T_{N} i P_{N} są zgodnie uporządkowane..

Rozwiązanie: 

Rzeczywiście, w przypadku macierzy trójdiagonalnej T_{N} wystarczy wziąć najpierw zmienne o indeksach nieparzystych, a potem — o indeksach parzystych: P(1,2,\ldots,N)^{T}=(1,3,\ldots,2,4,\ldots)^{T}. Dla macierzy P_{N} wybieramy podobnie (tzw. red-black ordering).

Twierdzenie 5.6 (Vargi o związku pomiędzy zbieżnością SOR i metody Jacobiego)

Niech A będzie macierzą zgodnie uporządkowaną i niech \omega>0. Wtedy jeśli \mu jest wartością własną macierzy iteracji metody Jacobiego, to -\mu też nią jest. Ponadto, jeśli \lambda jest różną od zera wartością własną macierzy iteracji metody SOR z parametrem relaksacji \omega oraz zachodzi

(\lambda+\omega-1)^{2}=\lambda\,(\omega\,\mu)^{2}

to \mu jest wartością własną macierzy iteracji Jacobiego. Na odwrót, jeśli \mu jest wartością własną macierzy iteracji Jacobiego, to \lambda zadane powyższym wzorem jest wartością własną macierzy iteracji SOR.

Dowód

Macierz iteracji metody Jacobiego to po prostu

-(\tilde{L}+\tilde{U})=J(-1)=J(1),

(ostatnia równość na mocy założenia). Stąd oczywiście wynika część pierwsza tezy.

Macierz iteracji metody SOR ma postać

G_{\omega}=(D+\omega L)^{{-1}}((1-\omega)D-\omega U)=(I-\omega\tilde{L})^{{-1}}((1-\omega)I+\omega\tilde{U}).

Niech teraz \lambda będzie wartością własną macierzy iteracji metody SOR,

\det(\lambda I-G_{\omega})=0.

Ponieważ \det(I-\omega\tilde{L})=1, mamy równoważnie

\displaystyle 0 \displaystyle=\det\left((I-\omega\tilde{L})(\lambda I-G_{\omega})\right)=\det\left(\lambda(I-\omega\tilde{L})-(1-\omega)I-\omega\tilde{U}\right)
\displaystyle=\det\left((\lambda+\omega-1)I-\lambda\omega\tilde{L}-\omega\tilde{U}\right).

Zatem \lambda jest wartością własną G_{\omega} wtedy i tylko wtedy, gdy \dfrac{\lambda+\omega-1}{\omega} jest wartością własną macierzy \lambda\tilde{L}+\tilde{U}.

W konsekwencji, jeśli tylko \lambda\neq 0, to \dfrac{(\lambda+\omega-1)}{\omega\sqrt{\lambda}} jest wartością własną macierzy J(\sqrt{\lambda}). Na mocy założenia, J(\sqrt{\lambda}) ma te same wartości własne, co J(-1), co kończy dowód drugiej części twierdzenia. Dowód części trzeciej pozostawiamy jako ćwiczenie.

Wniosek 5.2 (O zbieżności metody Gaussa–Seidela)

Jeśli A jest zgodnie uporządkowana i promień spektralny macierzy iteracji metody Jacobiego dla A jest równy \mu<1, to promień spektralny macierzy iteracji metody Gaussa–Seidela jest równy \mu^{2}. Znaczy to, że metoda Gaussa–Seidela jest dwukrotnie szybsza od metody Jacobiego.

Dowód

Rzeczywiście, na mocy twierdzenia Vargi, każdej niezerowej wartości własnej macierzy SOR odpowiada wartość własna macierzy iteracji Jacobiego \mu, spełniająca równanie

(\lambda+\omega-1)^{2}=\lambda(\omega\mu)^{2}.

Dla metody Gaussa–Seidela, \omega=1, i w konsekwencji \lambda=\mu^{2}.

Korzystając z twierdzenia Vargi można udowodnić

Twierdzenie 5.7 (Twierdzenie Younga o optymalnym parametrze SOR)

Niech macierz A będzie symetryczną macierzą zgodnie uporządkowaną i niech \mu<1 będzie równe promieniowi spektralnemu macierzy iteracji metody Jacobiego. Wtedy optymalna wartość parametru relaksacji \omega, gwarantująca najmniejszy promień spektralny macierzy iteracji metody SOR, jest równa

\omega _{{opt}}=\dfrac{2}{1+\sqrt{1-\mu^{2}}}.

Promień spektralny macierzy iteracji SOR z tym parametrem wynosi wtedy

\left(\dfrac{\mu}{1+\sqrt{1-\mu^{2}}}\right)^{2}.

Dowód powyższego twierdzenia można znaleźć na przykład w [15, rozdział 8.3].

Jakkolwiek piękne teoretycznie, nawet twierdzenia Younga i Vargi są dość bezlitosne dla metody SOR w wielu praktycznych zastosowaniach.

Przykład 5.11

Dla macierzy jednowymiarowego laplasjanu rozmiaru N, oznaczając \xi=\pi/(N+1), mamy \mu=\cos(\xi) i w konsekwencji \omega _{{opt}}=2/(1+\sin(\xi)). Zatem, dla N=100, gdy \mu\approx 0.9995, promień spektralny macierzy iteracji SOR z optymalnym parametrem relaksacji jest równy aż (1-\sin(\xi))/(1+\sin(\xi))\approx 0.9396, a więc jest praktycznie tak samo blisko jedności, jak w przypadku metody Jacobiego!

Przykład 5.12

Kontynuując przykład 5.9, porównamy szybkość zbieżności metody SOR, Gaussa–Seidela (czyli SOR z parametrem \omega=1) i Jacobiego na tym samym, ekstremalnie trudnym dla metod iteracyjnych zadaniu. Na początek weźmiemy N=20 i zbadamy zależność szybkości zbieżności SOR od wartości \omega. Korzystając z wyniku poprzedniego przykładu, będziemy mogli także uruchomić SOR z optymalnym parametrem.



Zwróć uwagę na dramatyczną przewagę SOR z optymalnym parametrem nad metodą Gaussa–Seidela. Zbadaj, jak zmienią się wyniki, gdy wyraźnie zwiększysz N, np. do tysiąca.

Ćwiczenie 5.15

Gdy macierz iteracji B metody stacjonarnej

x_{{k+1}}=Bx_{k}+d

jest daleka od normalnej, numeryczna realizacja iteracji w arytmetyce zmiennoprzecinkowej ograniczonej precyzji może napotkać na problemy: w początkowych iteracjach residuum może znacząco rosnąć, co w efekcie może doprowadzić do praktycznej utraty zbieżności.

Zbadaj to na przykładzie–zabawce (por. [2, 10.2.3]),

B=\begin{pmatrix}\lambda&1\\
0&\lambda\end{pmatrix},

gdy \lambda\in(0.51.0). Wyznacz promień spektralny tej macierzy, jej normę, a także wykres normy x_{k} w zależności od k.

Rozwiązanie: 

Dla \lambda\approx 1 (wyraźnie większych od 1), macierz iteracji B, choć o promieniu spektralnym mniejszym niż 1, ma tę właściwość, że norma jej kolejnych potęg początkowo rośnie.

lam = 0.9995; B = [lam 1; 0 lam];
resid=[]; x = [1;1];
for i = 1:100000,
	x = B*x; resid(i) = norm(x);
end;
semilogy(resid);

Ćwiczenie 5.16

Metoda Richardsona z parametrem \tau\in\mathbb{R} jest określona wzorem

x_{{k+1}}=x_{k}+\tau(b-Ax_{k}).

Niech macierz A będzie symetryczna.

  1. Sprawdź, że jest to stacjonarna metoda iteracyjna oparta na rozszczepieniu A=M-Z, gdzie M=\dfrac{1}{\tau}\, I.

  2. Wykaż, że jeśli A jest nieokreślona, dla żadnego wyboru \tau nie jest prawdą, że metoda Richardsona jest zbieżna do x^{*} dla dowolnego x_{0}.

  3. Wykaż, że jeśli A jest dodatnio określona, to metoda Richardsona jest zbieżna dla 0<\tau<2/\lambda _{{\max}}(A), a dla

    \tau _{{\text{opt}}}=\dfrac{2}{\lambda _{{\max}}(A)+\lambda _{{\min}}(A)},

    promień spektralny macierzy iteracji M^{{-1}}Z jest najmniejszy możliwy i równy

    \rho(M^{{-1}}Z)=\| M^{{-1}}Z\| _{2}=\dfrac{\kappa-1}{\kappa+1},\qquad\text{ gdzie }\kappa=\dfrac{\lambda _{{\max}}(A)}{\lambda _{{\min}}(A)}.
Ćwiczenie 5.17

Niech A_{{\operatorname{sym}}}=\dfrac{A+A^{T}}{2} i niech istnieją stałe 0<\lambda\leq\Lambda takie, że dla każdego x\in\mathbb{R}^{N}

\lambda x^{T}x\leq x^{T}A_{{\operatorname{sym}}}x\qquad\text{oraz}\qquad\| Ax\| _{2}^{2}\leq\Lambda x^{T}A_{{\operatorname{sym}}}x.

Wykaż, że wtedy metoda Richardsona jest zbieżna w normie euklidesowej dla 0<\tau<2/\Lambda. Ponadto, jeśli \tau=1/\Lambda, to współczynnik zbieżności \rho w tej normie można oszacować przez \sqrt{1-\dfrac{\lambda}{\Lambda}}.

Ćwiczenie 5.18 (Hackbusch)

Niech A,M będą macierzami symetrycznymi i dodatnio określonymi. Dla dwóch macierzy symetrycznych X,Y będziemy pisali X>Y, jeśli X-Y>0, czyli, innymi słowy, gdy x^{T}Xx\,>\, x^{T}Yx dla każdego x\neq 0. Wykaż, że:

  1. Jeśli

    2M>A,

    to iteracja

    x_{{k+1}}=x_{k}+M^{{-1}}(b-Ax_{k})

    jest zbieżna, a dokładniej, macierz iteracji, B=I-M^{{-1}}A, spełnia

    \rho(B)=\| B\| _{A}=\| B\| _{M}<1.
  2. Jeśli dla pewnych 0<\lambda\leq\Lambda zachodzi

    \lambda W\leq A\leq\Lambda W,

    to wszystkie wartości własne B=I-M^{{-1}}A są rzeczywiste oraz 1-\Lambda\leq\lambda(B)\leq 1-\lambda.

5.4. Metody blokowe

Pewnym remedium na powolną zbieżność metod stacjonarnych jest stosowanie metod blokowych, w których iterację konstruuje się tak jak dotychczas, ale zamiast na elementach, działa się na całych blokach macierzy (naturalnie, jeśli trzeba, musimy założyć ich odwracalność). Dla poprawienia zbieżności często stosuje się wariant z zakładką, zob. rysunek 5.6. W praktyce takie metody mogą być wyraźnie skuteczniejsze od dotychczas rozważanych: rzeczywiście, w skrajnym przypadku, pojedynczego ,,bloku” — rozmiaru N — wszystkie wyżej omówione metody stacjonarne redukują się do metody bezpośredniej…

Macierz A jest podzielona na kilka, częściowo nakładających się na siebie, bloków. Te bloki następnie są --- każdy z osobna --- odwracane.
Rys. 5.6. Schemat ideowy blokowej metody Jacobiego z zakładką: zamiast standardowego D^{{-1}}r=\sum _{{i=1}}^{N}e_{i}d_{{ii}}^{{-1}}e_{i}^{T}r wyznaczamy \sum _{{i=1}}^{K}R_{i}A_{{ii}}^{{-1}}R_{i}^{T}r, gdzie R_{i} jest macierzą operatora obcięcia r do współrzędnych odpowiadających A_{{ii}}.

Ciekawym i nietrywialnym uogólnieniem metod z zakładką jest addytywna metoda Schwarza: jeden z przedstawicieli klasy tzw. metod dekompozycji obszaru, będących bardzo efektywnymi i naturalnie równoległymi operatorami ściskającymi (mówimy o nich w rozdziale 7) dla dyskretyzacji równań różniczkowych cząstkowych, [14].

5.5. Metody projekcji

Dotychczas omówione metody opierają się na czysto formalnym, mechanicznym, podziale macierzy, bezpośrednio związanym z jej strukturą. Tymczasem można zdefiniować metody iteracyjne o podobnym koszcie iteracji, ale w których kolejne przybliżenia rozwiązania będziemy wybierać tak, by w jakimś sensie zminimalizować błąd.

Metodę iteracyjną będziemy definiować w formie aktualizacji poprzedniego przybliżenia,

x_{{k+1}}=x_{{k}}+\delta _{k}.

,,Idealna poprawka” \delta _{k}, oznaczmy ją \delta^{*}, powinna spełniać więc równanie

A\delta^{*}=r_{{k}},

bo wtedy dostalibyśmy x_{{k+1}}=x^{*}. Ma ona jednak tę wadę, że do jej wyznaczenia musielibyśmy dokładnie rozwiązać nasz układ równań (co prawda z nieco inną prawą stroną, ale to bez znaczenia!) — a więc wpadlibyśmy z deszczu pod rynnę! Jednak, gdyby można było tanio rozwiązać ten układ równań w przybliżeniu, to może w ten sposób udałoby się nam określić zbieżną metodę?

Aby określić przybliżone rozwiązanie równania poprawki, najpierw przeformułujemy je nieco. Niech dwie macierze: U i V mają tę własność, że kolumny każdej z nich tworzą bazę \mathbb{R}^{N}. Wtedy rozwiązanie dokładne równania idealnej poprawki, \delta^{*}, możemy reprezentować w bazie V, \delta^{*}=Va^{*} dla pewnego a^{*}\in\mathbb{R}^{N}. Co więcej, z nieosobliwości macierzy U i V wynika, że

A\delta^{*}=r_{{k}}\iff U^{T}AVa^{*}=U^{T}r_{{k}}.

To daje nam pomysł na definicję przybliżonej poprawki idealnej: niech U_{k} i V_{k} będą macierzami pełnego rzędu, tego samego zadanego z góry rozmiaru N\times d_{k}. Wtedy poprawkę \delta _{k} określimy jako wektor

\delta _{k}=V_{k}a_{k}

taki, że a_{k}\in\mathbb{R}^{{d_{k}}} spełnia równanie

\displaystyle U_{k}^{T}AV_{k}a_{k}=U_{k}^{T}r_{{k}}. (5.20)

Jest to właśnie metoda projekcji.

Ponieważ macierz zredukowana A_{k}=U_{k}^{T}AV_{k} jest kwadratowa rozmiaru d_{k}, to \delta _{k} będzie tanie do wyznaczenia, gdy d_{k} jest niewielkie. Nazwa metody bierze się stąd, że z powyższej definicji wynika, że wektor residualny równania poprawki, r_{{k}}-A\delta _{k} jest prostopadły do podprzestrzeni rozpiętej przez kolumny U_{k}:

U_{k}^{T}(r_{{k}}-A\delta _{k})=0.

Z drugiej strony, \delta _{k} będzie dobrze określone tylko wtedy, gdy macierz zredukowana A_{k} będzie nieosobliwa — co nie zawsze musi być prawdą, nawet jeśli A jest nieosobliwa (pomyśl o macierzy \begin{pmatrix}0&1\\
1&1\end{pmatrix} i U_{k}=V_{k}=\begin{pmatrix}1\\
0\end{pmatrix}). Aby zagwarantować sobie odwracalność macierzy zredukowanej, zwykle wybiera się macierz pełnego rzędu V_{k} i w zależności od własności macierzy A dobiera się macierz U_{k}:

  • Jeśli A=A^{T}>0, to kładziemy U_{k}=V_{k}. Wtedy macierz zredukowana A_{k}=V_{k}^{T}AV_{k} jest symetryczna i dodatnio określona.

  • Jeśli A jest tylko nieosobliwa, to kładziemy U_{k}=AV_{k}. Macierz zredukowana jest postaci A_{k}=V_{k}^{T}A^{T}AV_{k}, a więc jest macierzą zredukowaną poprzedniego rodzaju, ale dla macierzy równań normalnych, A^{T}A.

Jak możemy się domyślić, metody projekcji są metodami minimalizacji, co potwierdza poniższy dwuczęściowy lemat:

Lemat 5.1

Niech V_{k} będzie zadaną macierzą N\times d_{k}, pełnego rzędu. Oznaczmy przez \mathcal{V}_{k} podprzestrzeń rozpiętą przez kolumny V_{k}.

  • Jeśli A=A^{T}>0 oraz U_{k}=V_{k}, to x_{{k+1}} określone metodą projekcji (5.20) spełnia

    \| x^{*}-x_{{k+1}}\| _{A}\leq\| x^{*}-x\| _{A}\qquad\forall x\in x_{{k}}+\mathcal{V}_{k}.
  • Jeśli A jest nieosobliwa i U_{k}=AV_{k}, to x_{{k+1}} określone metodą projekcji (5.20) spełnia

    \| b-Ax_{{k+1}}\| _{2}\leq\| b-Ax\| _{2}\qquad\forall x\in x_{{k}}+\mathcal{V}_{k}.

Dowód zostawiamy jako ćwiczenie.

Ćwiczenie 5.19

Udowodnij lemat 5.1.

Wskazówka: 

Patrz [13] lub dowód twierdzenia 6.1.

5.5.1. Metoda najszybszego spadku

Jednym z bardziej prominentnych przykładów metody projekcji jest metoda najszybszego spadku, działająca w przypadku, gdy macierz A jest symetryczna i dodatnio określona. W tej metodzie wybieramy U_{k}=V_{k}=r_{{k}}=b-Ax_{{k}}. Ponieważ wymiar przestrzeni rozpiętej przez kolumny U_{k} jest równy d_{k}=1, równanie poprawki upraszcza się do jednego równania skalarnego na a_{k}\in\mathbb{R},

a_{k}=\dfrac{r_{{k}}^{T}r_{{k}}}{r_{{k}}^{T}Ar_{{k}}}

i w konsekwencji x_{{k+1}}=x_{{k}}+a_{k}r_{{k}}.

Nazwa metody wywodzi się stąd, że wektor poprawki w tej metodzie jest proporcjonalny do residuum, które z kolei jest kierunkiem gradientu funkcjonału \phi(x)=\| x^{*}-x\| _{A}^{2} w punkcie x_{k}:

\nabla\phi(x_{k})=b-Ax_{k}=r_{k}.
Twierdzenie 5.8 (o zbieżności metody najszybszego spadku)

W metodzie najszybszego spadku,

\| x^{*}-x_{{k+1}}\| _{A}\leq\dfrac{\kappa-1}{\kappa+1}\| x^{*}-x_{{k}}\| _{A},

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

Dowód

Łatwo wykazać (por. [13, twierdzenie 5.2]), że jeśli r_{{k}}\neq 0,

\| x^{*}-x_{{k+1}}\| _{A}^{2}\leq\left(1-\dfrac{r_{{k}}^{T}r_{{k}}}{r_{{k}}^{T}Ar_{{k}}}\dfrac{r_{{k}}^{T}r_{{k}}}{r_{{k}}^{T}A^{{-1}}r_{{k}}}\right)\| x^{*}-x_{{k}}\| _{A}^{2}.

Teza wynika z lematu, którego elegancki dowód, pochodzący od Braessa, można znaleźć w [2]:

Lemat 5.2 (Kantorowicza)

Niech A=A^{T}>0. Wtedy dla dowolnego x\neq 0,

\dfrac{x^{T}Ax}{x^{T}x}\cdot\dfrac{x^{T}A^{{-1}}x}{x^{T}x}\leq\dfrac{1}{4}(\sqrt{\kappa}+\dfrac{1}{\sqrt{\kappa}})^{2}.
Ćwiczenie 5.20

Wykaż, że w metodzie najszybszego spadku zachodzi r_{{k+1}}^{T}r_{k}=0.

Rozwiązanie: 

Ponieważ x_{{k+1}}=x_{{k}}+a_{k}r_{{k}}, to

r_{{k+1}}=r_{{k}}-a_{k}Ar_{{k}},

zatem

r_{k}^{T}r_{{k+1}}=r_{k}^{T}r_{{k}}-a_{k}r_{k}^{T}Ar_{{k}}=\| r_{k}\| _{2}^{2}-a_{k}r_{k}^{T}Ar_{{k}}=0

z definicji a_{k}.

Przykład 5.13

Kontynuujemy przykład 5.8. Chcąc porównywać trzy metody: Jacobiego, Gaussa–Seidela oraz metodę najszybszego spadku, musimy zadbać o to, żeby były spełnione warunki zbieżności tej ostatniej — a więc, aby macierz A była symetryczna i dodatnio określona. Dlatego, tym razem położymy, dla dodatniego p,

A=B^{T}B+pI.


Choć A wydaje się dosyć gęsta, w rzeczywistości wciąż jest macierzą rzadką. Zwróćmy uwagę nie tylko na samą szybkość zbieżności — mierzoną liczbą iteracji — ale także na efektywność metody: ile czasu zajmuje wyznaczenie przybliżenia z zadaną dokładnością. Metoda najszybszego spadku niewątpliwie ma najtańszą iterację (a najdroższą — metoda Gaussa–Seidela), co ostatecznie przekłada się na większą efektywność metody najszybszego spadku (pomimo mniejszej szybkości zbieżności).

Ćwiczenie 5.21

Wyjaśnij, dlaczego w powyższym przykładzie, dla małych wartości p, metoda Jacobiego czasem nie jest zbieżna.

Rozwiązanie: 

Bo macierz A nie musi być diagonalnie dominująca.

Przykład 5.14

Dla macierzy jednowymiarowego laplasjanu, mamy

\kappa=\dfrac{\lambda _{{\max}}(A)}{\lambda _{{\min}}(A)}=\left(\dfrac{\sin(\pi N/2(N+1))}{\sin(\pi/2(N+1))}\right)^{2}=O(N^{2}),

zatem dla N=100 mamy współczynnik redukcji błędu na poziomie (\kappa-1)/(\kappa+1)\approx 0.9995. Znaczy to dokładnie tyle, że metoda nie nadaje się do rozwiązywania takich zadań, gdy N jest duże. Szybkość zbieżności metody najszybszego spadku jest w naszym przykładzie porównywalna z szybkością metody Jacobiego i gorsza od metody SOR z optymalnym parametrem relaksacji.

5.5.2. Metoda najmniejszego residuum

Gdy o macierzy A wiemy jedynie, że jest nieosobliwa, możemy zastosować metodę najmniejszego residuum. W tej metodzie wybieramy V_{k}=r_{{k}}=b-Ax_{{k}} oraz U_{k}=AV_{k}=Ar_{k}. Równanie poprawki znów upraszcza się do jednego równania skalarnego na a_{k}\in\mathbb{R},

a_{k}=\dfrac{r_{{k}}^{T}A^{T}r_{{k}}}{r_{{k}}^{T}A^{T}Ar_{{k}}}

i w konsekwencji x_{{k+1}}=x_{{k}}+a_{k}r_{{k}}.

Twierdzenie 5.9

Załóżmy, że macierz A_{{\operatorname{sym}}}=(A+A^{T})/2 jest dodatnio określona i oznaczmy \mu=\lambda _{{\min}}(A_{{\operatorname{sym}}})>0. Wtedy

\| r_{{k+1}}\| _{2}\leq\left(1-\dfrac{\mu^{2}}{\| A\| _{2}^{2}}\right)^{{1/2}}\| r_{{k}}\| _{2}.

A_{{\operatorname{sym}}} nazywana jest częścią symetryczną macierzy A.

Ćwiczenie 5.22

Przeprowadź dowód twierdzenia o zbieżności metody najmniejszego residuum.

Wskazówka: 

Zob. dowód twierdzenia 5.3 w [13].

Ćwiczenie 5.23

Przypuśćmy, że umiemy tanio rozwiązywać układy równań z (nieosobliwą) macierzą A rozmiaru N\times N. Niech będzie dana macierz

\tilde{A}=\begin{pmatrix}A&v\\
u^{T}&\delta\end{pmatrix},

gdzie u,v\in\mathbb{R}^{N} oraz \delta\in\mathbb{R}.

  • Wskaż warunek konieczny i dostateczny na to, by macierz \tilde{A} była nieosobliwa.

  • Podaj możliwie tani algorytm rozwiązywania układu równań z macierzą \tilde{A}.

Wskazówka: 

\tilde{A} jest nieosobliwa wtedy i tylko wtedy, gdy \delta-u^{T}A^{{-1}}v\neq 0.

Ćwiczenie 5.24 (wzór Shermana–Morrisona)

Przypuśćmy, że umiemy tanio rozwiązywać układy równań z (nieosobliwą) macierzą A rozmiaru N\times N. Niech będzie dana macierz

\tilde{A}=A-vu^{T}.

gdzie u,v\in\mathbb{R}^{N}.

  • Wskaż warunek konieczny i dostateczny na to, by macierz \tilde{A} była nieosobliwa.

  • Podaj możliwie tani algorytm rozwiązywania układu równań z macierzą \tilde{A}.

Wskazówka: 

Zauważ, że \tilde{A}x=b wtedy i tylko wtedy, gdy \begin{pmatrix}x\\
t\end{pmatrix} jest rozwiązaniem układu

\begin{pmatrix}A&v\\
u^{T}&1\end{pmatrix}\begin{pmatrix}x\\
t\end{pmatrix}=\begin{pmatrix}b\\
0\end{pmatrix}.

Patrz także stwierdzenie 11.1.

Quiz 5.1
  1. diagonalna



    Dobrze Oczywiście, macierz diagonalną możemy rozwiązać bezpośrednio kosztem N flopów — a więc optymalnym. Wiele metod jednak sprowadziłoby się do metody bezpośredniej — na przykład metoda Jacobiego…

    Źle Zastanów się, ile kosztowałoby rozwiązanie takiego układu równań metodą bezpośrednią, wykorzystującą specjalną postać macierzy

  2. trójdiagonalna



    Dobrze Rozkład LU macierzy trójdiagonalnej możemy wyznaczyć bezpośrednio kosztem O(N) flopów, podobnie z rozwiązaniem samego układu (macierz L jest dwudiagonalna, a macierz U — co najwyżej trójdiagonalna). Stałe przy tym nie są zbyt duże

    Źle Zastanów się, ile kosztowałoby rozwiązanie takiego układu równań metodą bezpośrednią, wykorzystującą specjalną postać macierzy

  3. pełna, tzn. bez zerowych elementów



    Dobrze Ponieważ macierz jest pełna, to jedno mnożenie wektora przez A kosztuje O(N^{2}) flopów, więc metoda iteracyjna będzie miała sens, gdy satysfakcjonujące przybliżenie dostaniemy po k\leq N iteracjach

    Źle Rozkład LU macierzy pełnej możemy wyznaczyć bezpośrednio kosztem O(N^{3}) flopów, co zwykle jest ponad siły współczesnego komputera, gdy N jest bardzo duże

Ćwiczenie 5.25

Czy metodę iteracyjną warto stosować do rozwiązywania układu Ax=b z macierzą N\times N, w przypadku, gdy N jest bardzo duże oraz A jest rzadka?

Rozwiązanie: 

Tutaj, jak wynika z naszych dotychczasowych rozważań, odpowiedź jest zniuansowana. Jedno mnożenie wektora przez A kosztuje zapewne O(N) flopów, więc metoda iteracyjna będzie miała sens, gdy satysfakcjonujące przybliżenie dostaniemy po k\ll N iteracjach (a to, jak wiemy, nie zawsze jest gwarantowane). Z drugiej strony, może istnieć skuteczny reordering macierzy (czyli zmiana uporządkowania niewiadomych lub równań), pozwalający tanio wyznaczyć rozwiązanie metodą bezpośrednią [13].

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.