Zagadnienia

5.1 Macierze rozrzedzone
- 5.1.1 Przykłady macierzy rozrzedzonych
5.2 Macierze rozrzedzone: implementacja
- Format współrzędnych
- Format spakowanych kolumn (lub wierszy)
5.3 Metody stacjonarne rozwiązywania układów równań liniowych
5.4 Metody blokowe
5.5 Metody projekcji
- 5.5.1 Metoda najszybszego spadku
- 5.5.2 Metoda najmniejszego residuum

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 regularna¹Zazwyczaj 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).

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)²Mniej 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).

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

w formacie współrzędnych,
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 Toeplitza³Macierzą 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.

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

function [x, info, iter, resid] = jacobi(A, b, tol=1e-8, maxit=200, x0)
%
% metoda Jacobiego na macierzy A
%
if (nargin < 5)
	x0 = zeros(size(b));
end
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy na diagonale, scisle dolny trojkat i scisle gorny: 
% A = L + diag(D) + U; LpU = L + U
D = diag(A); LpU = tril(A,-1) + triu(A,1);

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

N = 2500; % 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);
spy(B); pause(0.5);

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;
	semilogy(resid,['b-;p=',num2str(p),';']); pause(0.5);
end

% dla porownania - solver bezposredni

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

Ć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).

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

function [x, info, iter, resid] = jacobi(A, b, tol=1e-8, maxit=1000, x0)
%
% metoda Jacobiego na macierzy A
%
if (nargin < 5)
	x0 = zeros(size(b));
end
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy na diagonale, scisle dolny trojkat i scisle gorny: 
% A = L + diag(D) + U; LpU = L + U
D = diag(A); LpU = tril(A,-1) + triu(A,1);

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

N = 25; % rozmiar macierzy
b = rand(N,1);
e = ones(N,1); B = spdiags([-e, 2*e, -e], [-1,0,1], N,N);

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

% testy dla roznego stopnia diagonalnej dominacji macierzy A

prange = [0, 1, 5, 10, 50]; % zestaw wartosci testowych parametru p
spy(B); pause(0.5);

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;
	semilogy(resid,['b-;p=',num2str(p),';']); pause(0.5);
end

% dla porownania - solver bezposredni

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

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.

%
% 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

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

N = 2500; % 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);
spy(B); pause(0.5);

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); 
	iterg(k) = its; timeg(k) = toc;
	mi = max(iter(k),iterg(k));
	semilogy(1:mi,[resid(1:mi),residg(1:mi)],['b-;p=',num2str(p),' Jacobi;';'r-;p=',num2str(p),' (Gauss-Seidel);']); 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 Gaussa-Seidela (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d | %4d | %10.5f\n', [prange; iterg; timeg]);
printf('Metoda bezposrednia (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d %4d %10.5f\n', [prange(end); 0; timed]);

Ć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)⁴W 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.

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

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

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

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

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

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

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

prange = linspace(1,1.9,5); % zestaw wartosci testowych parametru omega
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

hold off;
X = A\b;

[x, info, its, resid] = jacobi(A, b); iter = its; 
semilogy(1:its, resid(1:its),['.k;Jacobi;']); hold on; pause(0.5);
[x, info, its, residg] = sor(A, b, omega); iter = its; 
semilogy(1:its, residg(1:its),['-k;omega=',num2str(omega),' (SOR opt);']); hold on; pause(0.5);
for k = 1:length(prange)
	p = prange(k);
	tic; [x, info, its, residg] = sor(A, b, p); 
	iterg(k) = its; timeg(k) = toc;
	semilogy(1:its, residg(1:its),['-',num2str(k),';omega=',num2str(p),' (SOR);']);
	pause(0.5);
end
hold off;

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.

Sprawdź, że jest to stacjonarna metoda iteracyjna oparta na rozszczepieniu $A=M-Z$ , gdzie $M=\dfrac{1}{\tau}\, I$ .
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}$ .
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:

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.$
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.$

%
% 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 na diagonale, scisle dolny trojkat i scisle gorny: 
% A = L + diag(D) + U;
D = diag(A); L = tril(A,-1); U = triu(A,1);

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

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

function [x, info, iter, resid] = gaussseidel(A, b, tol=1e-8, maxit=100, x0)
%
% metoda Gaussa-Seidela 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 na diagonale, dolny trojkat i scisle gorny: 
% A = L + U;
L = tril(A); U = triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = L\(b - U*x);
	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 = 1500; % 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;
spy(B); pause(0.5);

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] = gaussseidel(A, b); 
	iterg(k) = its; timeg(k) = toc;
	tic; [x, info, its, resids] = steepestdescent(A, b); 
	iters(k) = its; times(k) = toc;
	mi = max([iter(k),iterg(k),iters(k)]);
	semilogy(1:mi,[resid(1:mi), residg(1:mi), resids(1:mi)],['b-;p=',num2str(p),' Jacobi;';'r-;p=',num2str(p),' (Gauss-Seidel);';'g-;p=',num2str(p),' (najszybszy spadek);']); 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 Gaussa-Seidela (p, iter, czas)\n');
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 bezposrednia (p, iter, czas)\n');
printf('  p  | iter |   czas\n====================================\n');
printf('%4d %4d %10.5f\n', [prange(end); 0; timed]);

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

diagonalna
Tak

Nie
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
trójdiagonalna
Tak

Nie
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
pełna, tzn. bez zerowych elementów
Tak

Nie
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.