MN08

Wielkie układy równań liniowych

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Wraz z coraz większymi modelami pojawiającymi się w praktyce obliczeniowej, coraz częściej zachodzi potrzeba rozwiązywania zadań algebry liniowej, w której macierze są co prawda wielkiego wymiaru, ale najczęściej rozrzedzone, to znaczy jest w nich bardzo dużo zer. Bardzo często zdarza się, że macierz wymiaru $N$ ma tylko $O(N)$ niezerowych elementów. Wykorzytanie tej specyficznej własności macierzy nie tylko prowadzi do algorytmów istotnie szybszych od ich analogów dla macierzy gęstych (to znaczy takich, które (w założeniu) mają $N^2$ elementów), ale wręcz są jedynym sposobem na to, by niektóre zadania w ogóle stały się rozwiązywalne przy obecnym stanie techniki obliczeniowej!

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

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.

Z reguły zadania liniowe wielkiego wymiaru będą miały strukturę macierzy rozrzedzonej, gdyż najczęściej związki pomiędzy niewiadomymi w równaniu nie dotyczą wszystkich, tylko wybranej grupy.

Przykład: 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. Pochodzi z kolekcji Tima Davisa. Jest to mała macierz, wymiaru 8032 (w kolekcji spotkasz równania z milionem i więcej niewiadomych).

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.

Reprezentacja macierzy rzadkich

Zacznijmy od sposobu przechowywania macierzy rozrzedzonych. Naturalnie, nie ma sensu przechowywać wszystkich zerowych jej elementów: wystarczy ograniczyć się do zachowania tych niezerowych! W ten sposób zmniejszamy zarówno wymagania pamięciowe, jak i liczbę operacji zmiennoprzecinkowych potrzebnych do prowadzenia działań na macierzy (np. w przypadku mnożenia macierzy przez wektor, nie będziemy mnożyć przez zera!).

Format współrzędnych

Do zapamiętania macierzy $A$ wymiaru $N\times M$ o $NNZ$ niezerowych elementów wykorzystujemy trzy wektory: I, J --- oba typu int --- oraz V, typu double, wszystkie o długości $NNZ$, przy czym

$A_{I[k],J[k]}=V[k],\mathrm{\forall }k=0,\dots ,NNZ-1.$

Format AIJ

W tym formacie wprowadzane są macierze rzadkie do Octave'a i MATLABa:

A = sparse(I,J,V,N,N);

Jednak wewnętrznie przechowywane są w innym formacie, o którym poniżej.

Przykład

Pokażemy jak w Octave wprowadzić macierz rozrzedzoną.

octave:10> I = [1, 1, 1, 2, 3, 1, 4]; octave:11> J = [1, 3, 2, 2, 3, 4, 4]; octave:12> V = [1, 2, 3, 4, 5, 6, 7]; octave:13> N = 4; octave:14> A = sparse(I,J,V,N,N) A = Compressed Column Sparse (rows = 4, cols = 4, nnz = 7) (1, 1) -> 1 (1, 2) -> 3 (2, 2) -> 4 (1, 3) -> 2 (3, 3) -> 5 (1, 4) -> 6 (4, 4) -> 7 octave:15> spy(A);

Strukturę jej niezerowych elementów ukaże nam polecenie spy(A). Tak właśnie zostały wygenerowane obrazki macierzy w niniejszym wykładzie.

Format spakowanych kolumn (wierszy)

Format współrzędnych nie narzucał żadnego uporządkowania elementów macierzy --- można było je umieszczać w dowolnej kolejności. 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, dobrze byłoby przechowywać elementy macierzy wierszami. Tak właśnie jest zorganizowany format spakowanych wierszy (Compressed Sparse Row, CSR). Analogicznie jest zdefiniowany format spakowanych kolumn (Compressed Sparse Column, CSC), 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: AV jest wektorem typu double o długości $NNZ$, zawierającym kolejne niezerowe elementy macierzy $A$ wpisywane kolumnami, AI jest wektorem typu int o długości $NNZ$, zawierającym numery wierszy macierzy $A$ odpowiadających elementom z AV. Natomiast zamiast tablicy J, jak to było w formacie współrzędnych, mamy krótszy wektor typu int, AP, o długości $M$, zawierający na $j$-tym miejscu indeks pozycji w AV, od którego rozpoczynają się w AV elementy $j$-tej kolumny macierzy $A$.

Format CSC

Mamy więc zależność, przy założeniu, że $AP[0]=0$,

$A_{AI[AP[j]+k],j+1}=AV[AP[j]+k],j=0,\dots ,M-1,k=0,\dots ,AP[j+1]-1.$

Taki (z drobnymi modyfikacjami) format macierzy wykorzystują np. pakiety ARPACK i UMFPACK (a także, wewnętrznie, Octave i MATLAB).

Format diagonalny

Znacznie mniej uniwersalny niż poprzednie i dlatego rzadziej spotykany. Kolejne diagonale macierzy przechowujemy w kolejnych wierszach macierzy $P\times N$, gdzie $P$ jest liczbą niezerowych diagonal w $A\in R^{N\times N}$.

Szczególnie wygodny do reprezentacji macierzy taśmowych. Wykorzystywany m.in. przez funkcję LAPACKa DGBSV służącą rozwiązywaniu równań z macierzami taśmowymi.

Uwagi praktyczne

Mnożenie macierzy w formacie AIJ przez wektor jest kilka razy wolniejsze w porównaniu do macierzy w formacie CSC/CSR (z tych dwóch oczywiście (dlaczego?) CSR jest najszybszy). Co więcej, okazuje się, że w typowych implementacjach, mnożenie macierzy rozrzedzonej (reprezentowanej np. w formacie CSC) przez wektor jest mało efektywne, mniej więcej na poziomie dziesięciu procent możliwości obliczeniowych procesora.

Jeśli już poważnie myśleć o przyspieszeniu mnożenia macierzy przez wektor (np. gdy chcemy stosować iteracyjną metodę rozwiązywania układu równań z tą macierzą), warto rozważyć inne formaty --- łączące w sobie podział macierzy na bloki (tak, jak w algorytmach BLAS) i przechowywanie (w zasadzie) tylko niezerowych elementów macierzy.

Macierze specjalne

Zajmiemy się teraz zadaniem rozwiązywania układu równań liniowych

$Ax=b,$

ale w sytuacji, gdy macierz $A$ jest rozrzedzona i dużego wymiaru. Dokonamy przeglądu kilku rodzajów algorytmów mających na celu wykorzystanie rozrzedzenia macierzy dla obniżenia kosztu wyznaczenia rozwiązania układu.

Należy pamiętać, że z reguły najlepsze wyniki uzyskuje się, gdy konkretny algorytm dobierze się do konkretnej macierzy. W zastosowaniach pojawiają się m.in. macierze rzadkie o bardzo szczególnej strukturze, dla nich warto stosować wyspecjalizowane algorytmy.

Macierze taśmowe

Macierz $A$ taka, że dla $|i-j|\ge d$, $a_{ij}=0$, nazywamy macierzą taśmową z rozstępem $d$, o szerokości pasma $2d+1$.

Łatwo sprawdzić, że algorytm eliminacji Gaussa (bez wyboru elementu głównego) nie spowoduje dodatkowego wypełnienia w takiej macierzy (a więc da się wykonać w miejscu). W przypadku konieczności wyboru elementu głównego, pesymistyczne oszacowanie rozstępu macierzy rozkładu LU jest równe $\max \{2d,N\}$ --- tak więc, dla niezbyt dużych $d$ wciąż wynikowa macierz jest taśmowa.

W szczególności, gdy macierz jest taśmowa z pasmem o rozstępie $d$ i jednocześnie diagonalnie dominująca, wtedy rozkład LU takiej macierzy da się wykonać w miejscu kosztem $O(dN)$, czyli liniowym.

W LAPACKu zaimplementowano szybki solver równań z macierzami taśmowymi, DGBSV, ale wymagający specjalnego sposobu przechowywania macierzy, wykorzystującego format diagonalny.

Macierze trójdiagonalne

Szczególnym przypadkiem macierzy taśmowych są macierze trójdiagonalne, tzn. taśmowe o rozstępie $d=1$:

$A=\left(\begin{array}{ccccc}{a}_1& c_1& & & \\ b_2& a_2& c_2& & \\ & b_3& a_3& \ddots & \\ & & \ddots & \ddots & c_{N-1}\\ & & & b_N& a_N\end{array}\right)$

Zadanie rozwiązywania równań z taką macierzą

$Ax=e$

jest bardzo często spotykane, więc warto przytoczyć algorytm w całej okazałości --- popularnie zwany algorytmem przeganiania (niektóre źródła nazywają go algorytmem Thomasa).

Zacznijmy jednak od stwierdzenia, kiedy macierz trójdiagonalna nie wymaga wyboru elementu głównego:

Stwierdzenie

Jeśli macierz $A$ ma słabo dominującą diagonalę, tzn.

$|a_i|\ge |b_i|+|c_i|,1\le i\le N,$

($b_1=0=c_N$) i przynajmniej dla jednego indeksu "$i$" mamy powyżej ostrą nierówność "$>$", to algorytm przeganiania jest wykonalny bez przestawień wierszy. Ponadto wymaga on $7N$ operacji arytmetycznych, a więc jest prawie optymalny.

Algorytm Metoda przeganiania

---

<math>\displaystyle d_1</math> = <math>\displaystyle a_1</math>; 
<math>\displaystyle f_1</math> = <math>\displaystyle e_1</math>;
for (i = 2; i <= N; i++)
{
	<math>\displaystyle l</math> = <math>\displaystyle b_i</math>/<math>\displaystyle a_{i-1}</math>;
	<math>\displaystyle d_i</math> = <math>\displaystyle a_i</math> - <math>\displaystyle l</math> * <math>\displaystyle c_{i-1}</math>;
	<math>\displaystyle f_i</math> = <math>\displaystyle e_i</math> - <math>\displaystyle l</math> * <math>\displaystyle f_{i-1}</math>;
} 
<math>\displaystyle x_1</math> = <math>\displaystyle f_N</math>;
for (i = N-1; i >= 1; i--)
	<math>\displaystyle x_i</math> = <math>\displaystyle f_i</math> - <math>\displaystyle c_i</math> * <math>\displaystyle x_{i+1}</math>;

Metody bezpośrednie

Przykład: Strzałka Wilkinsona

Rozważmy układ równań z macierzą diagonalnie dominującą postaci

$A=\left(\begin{array}{ccccc}*& *& *& \cdots & *\\ *& *& & & \\ *& & *& & \\ ⋮& & & \ddots & \\ *& & & & *\end{array}\right)$

gdzie $*$ oznacza jakiś niezerowy element. Łatwo sprawdzić, że chociaż wyjściowa macierz jest rozrzedzona, to zastosowanie do niej eliminacji Gaussa powoduje, że w wyniku dostajemy gęste czynniki rozkładu.

Tymczasem wystarczy odwrócić kolejność równań i numerację niewiadomych (co dla macierzy jest równoznaczne z odwróceniem porządku wierszy i kolumn, korzystając z pewnej (jakiej?) macierzy permutacji $P$):

$\tilde{A}=PA{P}^T=\left(\begin{array}{ccccc}*& & & & *\\ & \ddots & & & ⋮\\ & & *& & *\\ & & & *& *\\ *& \cdots & *& *& *\end{array}\right)$

Wtedy okazuje się, że rozkład naszej macierzy nie powoduje już wypełnienia czynników rozkładu!

Właśnie na tym polega główny problem w rozwiązywaniu układów z macierzami rzadkimi metodami bezpośrednimi: jak maksymalnie wykorzystać rozrzedzenie macierzy tak, by czynniki rozkładu były możliwie mało wypełnione. Albowiem wiedząc to będziemy mogli ograniczyć się jedynie do fizycznego wyznaczenia wartości niezerowych elementów macierzy rozkładu. Ponadto wymagania pamięciowe algorytmu nie będą istotnie wykraczać ponad ilość pamięci potrzebnej na przechowanie danych (tzn. macierzy).

W ogólnym przypadku rozwiązanie takiego zadania jest trudne i większość algorytmów opiera się na pewnych heurystykach, które jednak w praktyce warto wspomóc wcześniejszą analizą konkretnego układu równań jaki mamy rozwiązać. Najczęściej dąży się do takiego przenumerowania równań i niewiadomych, by w efekcie z góry przewidzieć, gdzie wystąpią zera w macierzach rozkładu --- i, by takich zer było jak najwięcej (by wypełnienie było jak najmniejsze)! Na architekturach z pamięcią hierarchiczną dąży się także do tego, by w trakcie rozkładu można było korzystać z BLAS Level 3, a więc permutacje wierszy i kolumn macierzy muszą to także brać pod uwagę.

Stosuje się kilka strategii wyznaczania korzystnych permutacji (reorderingu), z których warto wymienić

• przybliżone algorytmy minimalnego stopnia (approximate minimum degree, AMD), np. AMD;
• techniki podziału grafów na (prawie) rozłączne składowe nested dissection, ND, np. METIS.

Używają ich biblioteki takie jak UMFPACK, czy HSL.

W Octave mamy do dyspozycji także kilka procedur generujących takie permutacje, w tym: colamd (AMD dla macierzy niesymetrycznych) oraz symamd (AMD dla macierzy symetrycznych). Większy wybór oferuje MATLAB, jednak należy bezwzględnie pamiętać o jednym: nie ma uniwersalnej metody reorderingu i dla konkretnej macierzy może istnieć specjalna metoda, która da oszałamiające rezultaty, podczas gdy standardowe podejścia nie dadzą efektu.

Przykład: Rozwiązywanie układu z macierzą rozrzedzoną w Octave

Najprostszy sposób na wykorzystanie metody bezpośredniego rozwiązywania układu z macierzą rzadką to zastosowanie znanego nam operatora do macierzy typu rozrzedzonego:

A = sparse(...);
x = A \ b;

Octave domyślnie przyłoży do $A$ reordering colamd i następnie skorzysta z biblioteki UMFPACK, by rozwiązać taki układ. Dodatkowo, badane jest wcześniej, czy macierz nie jest przypadkiem taśmowa o wąskiej wstędze: jeśli jest, to korzysta się z LAPACKa.

(Jak widzisz, operator rozwiązywania układu równań jest mocno przeciążony: działa i na macierzach kwadratowych, i na prostokątnych, i na rozrzedzonych --- na każdym rodzaju w inny sposób. Nie dziwi Cię chyba, dlaczego wygodnie było napisać Octave w C++?)

Przykład: Wypełnienie pewnej macierzy w zależności od użytego reorderingu

Rozważmy ponownie macierz pochodzącą z kolekcji Tima Davisa. Jest to macierz symetryczna, dodatnio określona, wymiaru 8032, o 355460 elementach niezerowych i, w konsekwencji, o wypełnieniu około 0.6 procent.

Struktura niezerowych elementów macierzy.

Zobaczymy, jak w zależności od użytego algorytmu permutacji kolumn i wierszy poradzi sobie algorytm rozkładu Cholesky'ego.

Czynnik rozkładu Cholesky'ego $L$ wykonanego standardowym algorytmem. Czas rozkładu: 0.892013

Czynnik rozkładu Cholesky'ego $L$ z reorderingiem COLAMD. Czas rozkładu: 0.813038

Czynnik rozkładu Cholesky'ego $L$ z reorderingiem SYMAMD. Czas rozkładu: 0.487683s. Prawie dwa razy szybciej niż bez reorderingu, chociaż i tak wskutek wzrostu wypełnienia macierzy w dolnym trójkącie mamy aż 2 procent niezerowych elementów.

Czynnik rozkładu Cholesky'ego $L$ z odwrotnym reorderingiem Cuthill-McKee. Czas rozkładu: 0.845928

Czynnik rozkładu Cholesky'ego $L$ z jeszcze innym (bardzo tanim, ale jak widać czasem zupełnie złym) reorderingiem. Czas rozkładu: 5.947936

Na zakończenie popatrzmy, jak ważne jest spostrzeżenie symetrii macierzy:

Rozkład LU, czynnik $L$ (bez reorderingu). Czas rozkładu LU: 1.696018s. Nieakceptowalny, podobnie jak drastyczne wypełnienie macierzy.

Rozkład LU, czynnik $U$ (bez reorderingu)

Jak widać, w naszym przypadku standardowe algorytmy (COLAMD i SYMAMD) poradziły sobie całkiem nieźle, chociaż wypełnienie i tak znacząco wzrosło. Zapewne, dla tej macierzy, która jest specyficznego typu --- pochodzi z dyskretyzacji równania różniczkowego --- algorytm ND mógłby tu jeszcze lepiej zadziałać.

Stacjonarne metody iteracyjne

Gdy macierz jest rozrzedzona, mnożenie takiej macierzy przez wektor jest bardzo tanie (koszt jest proporcjonalny do liczby niezerowych elementów macierzy). Dlatego, jeśli możemy zadowolić się rozwiązaniem przybliżonym układu, a w zamian osiągnąć je tanim kosztem, warto rozważyć metody iteracyjne.

Najprostsze metody iteracyjne (najprostsze w analizie i implementacji, ale --- jak można się domyślić --- w praktyce najmniej efektywne) polegają na rozkładzie macierzy na część "łatwo odwracalną", $M$, i "resztę", $Z$. Dokładniej, jeśli $M$ jest nieosobliwa, to równanie $Ax=b$ można zapisać jako zadanie punktu stałego

$Mx=Nx+b,$

gdzie $Z=M-A$. Inaczej:

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

i zastosować doń metodę iteracji prostej Banacha:

$x_{k+1}=M^{-1}(Z{x}_k+b).$

Takie metody nazywamy stacjonarnymi metodami iteracyjnymi. Aby przeanalizować zbieżność takiej metody, warto rozpatrzyć przypadek ogólniejszy

$x_k=B{x}_{k-1}+c,$

dla pewnej macierzy $B$ oraz wektora $c$. (Dla stacjonarnej metody iteracyjnej, $B=M^{-1}Z$ oraz $c=M^{-1}b$.)

W tym przypadku

$x_k-x^{*}=B^k(x_0-x^{*}),$

a stąd i z nierówności $\Vert B^k\Vert \le \Vert B{\Vert }^k$, mamy

$\Vert x_k-x^{*}\Vert \le \Vert B{\Vert }^k\Vert x_0-x^{*}\Vert .$

Warunkiem dostatecznym zbieżności iteracji prostych jest więc $\Vert B\Vert <1$. Okazuje się, że warunkiem koniecznym i dostatecznym zbieżności tej iteracji dla dowolnego wektora startowego $x_0$ jest

$\rho =\max \{|\lambda |:\lambda \text{ jest wartością własną }B\}<1.$

Tak więc, metody oparte na iteracji prostej będą zbieżne liniowo z ilorazem $\rho$.

Zaletą stacjonarnych metod iteracyjnych jest również ich prostota, przez co są one łatwe do zaprogramowania, co łatwo zobaczyć na przykładach metod: Jacobiego i Gaussa--Seidela, które teraz omówimy.

Metoda Jacobiego

Biorąc $M=\text{diag}(A)$, gdzie $\text{diag}(A)$ jest macierzą diagonalną składającą się z wyrazów stojących na głównej przekątnej macierzy $A$, układ $Ax=b$ jest równoważny układowi

$Mx=Zx+b,$

a stąd (o ile na przekątnej macierzy $A$ nie mamy zera) otrzymujemy metodę iteracyjną

$x_k=B{x}_{k-1}+c,$

gdzie $B=-M^{-1}Z$ i $c=M^{-1}b$, zwaną metodą Jacobiego.

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

$x_i^{(k)}=\frac{1}{a_{ii}}(b_i-\sum _{j\ne i}{a}_{ij}{x}_j^{(k-1)}),$

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

Twierdzenie O zbieżności metody Jacobiego

W metodzie Jacobiego warunek dostateczny zbieżności, $\Vert B\Vert <1$, jest spełniony np. wtedy, gdy macierz $A$ ma dominującą przekątną, tzn. gdy

$|a_{i,i}|>\sum _{j\ne i}|a_{i,j}|,\mathrm{\forall }i=1,\dots ,N.$

Dowód

Rzeczywiście, ponieważ wyraz $(i,j)$ macierzy $M^{-1}Z$ wynosi $0$ dla $i=j$ oraz $a_{i,j}/a_{i,i}$ dla $i\ne j$, a więc

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \lefteqn{\|M^{-1}Z\|_\infty \;=\; \max_{1\le i\le N} \sum_{j=1,j\ne i}^N{|a_{i,j}|}/{|a_{i,i}|} } \\ && =\;\max_{1\le i\le N} \sum_{j=1}^N{|a_{i,j}|}/{|a_{i,i}|}\,-\,1\;<\;1, \endaligned}

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

Przykład: Macierz laplasjanu

Macierz $N\times N$, zwana macierzą jednowymiarowego laplasjanu

$L=\left(\begin{array}{cccc}2& -1& & \\ -1& 2& \ddots & \\ & \ddots & \ddots & -1\\ & & -1& 2\end{array}\right)$

pojawia się w bardzo wielu zastosowaniach, także jako podzadanie w 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)$. Stosując do niej metodę Jacobiego mamy $M=2I$ oraz $Z=L-M$. Obliczając normę macierzy iteracji Jacobiego dostajemy $||M^{-1}Z|{|}_{\mathrm{\infty }}=1$, co nie rozstrzyga jeszcze o jej zbieżności lub niezbieżności.

Potrzebna będzie bardziej subtelna analiza. Okazuje się, że są znane wzory na wartości własne macierzy $L$:

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

dla $1\le j\le N.$ W konsekwencji, wartościami własnymi $M^{-1}Z=\frac{1}{2}Z=\frac{1}{2}L-I$ są liczby ${\mu }_i=\frac{1}{2}{\lambda }_i-1$. Ponieważ $0<{\mu }_i<1$, znaczy to, że metoda Jacobiego jest zbieżna dla macierzy $L$.

Z drugiej strony, nie dajmy się zwieść optymizmowi matematyka ("nie martw się, jest zbieżny..."): nietrudno sprawdzić, że $||Z|{|}_2=1-O(N^{-2})<1$, co oznacza, że metoda Jacobiego --- choć zbieżna --- dla dużych $N$ staje się zbieżna tak wolno, że w praktyce bezużyteczna!

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ż rozwiązując równanie metody Jacobiego:

$x_i^{(k)}=\frac{1}{a_{ii}}(b_i-\sum _{j\ne i}{a}_{ij}{x}_j^{(k-1)}),$

nietrudno zauważyć, że w części sumy moglibyśmy odwoływać się do "dokładniejszych" wartości $x_j^{(k)}$: dla $j<i$, tzn.

$x_i^{(k)}=\frac{1}{a_{ii}}(b_i-\sum _{j<i}{a}_{ij}{x}_j^{(k)}-\sum _{j>i}{a}_{ij}{x}_j^{(k-1)}).$

W języku rozkładu macierzy $A=M-Z$ i iteracji $x_{k+1}=M^{-1}(Z{x}_k+b)$ mamy więc $M=\text{tril}(A)$ (dolny trójkąt macierzy $A$).

Twierdzenie O zbieżności metody Gaussa-Seidela

Jeśli macierz $A$ jest diagonalnie dominująca, to metoda Gaussa--Seidela jest zbieżna dla dowolnego wektora startowego $x_0$.

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

Metoda Gaussa--Seidela jest w wielu przypadkach rzeczywiście szybciej zbieżna od metody Jacobiego, np. tak jest w przypadku macierzy jednowymiarowego Laplasjanu. Wciąż jednak, dodajmy, dla zadań bardzo źle uwarunkowanych jej zbieżność jest zbyt wolna by ją stosować jako samodzielną metodę.

Uwaga

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.

Złożoność stacjonarnych metod iteracyjnych

Zastanówmy się teraz nad złożonością metod iteracyjnych. Ponieważ możemy jedynie znaleźć pewne przybliżenie rozwiązania dokładnego $x^{*}$, przez złożoność metody będziemy rozumieli koszt kombinatoryczny obliczenia $x_k$ z zadaną dokładnością $ϵ>0$. Dla uproszczenia założymy, że medoda jest zbieżna liniowo z ilorazem $\rho$. Zauważmy, że aby zredukować błąd początkowy do $ϵ>0$, wystarczy wykonać $k=k(ϵ)$ iteracji, gdzie $k$ spełnia

${\rho }^k\Vert x_0-x^{*}\Vert \le ϵ,$

czyli

$k\ge \frac{\log (1/ϵ)-\log (1/\Vert x_0-x^{*}\Vert )}{\log (1/\rho )}.$

Liczba ta zależy więc w istotny sposób od błędu początkowego i (przede wszystkim) od współczynnika redukcji błędu $\rho$, natomiast zależność od dokładności $ϵ$ i wymiaru $N$ układu jest dużo mniej istotna (w zadaniach praktycznych -- takich jak jednowymiarowy laplasjan --- jednak często okazuje się, że... $\rho$ zależy od $N$!).

Zakładając, że koszt jednej iteracji wynosi $c=c(N)$ ($c(N)$ jest tym mniejszy, im mniejsza jest liczba niezerowych elementów macierzy $A$), złożoność metody jest proporcjonalna do

$c(N)\frac{\log (1/ϵ)}{\log (1/\rho )}.$

Stąd oczywisty wniosek (prawdziwy nie tylko dla metod stacjonarnych), że metody iteracyjne warto stosować zamiast metod bezpośrednich w przypadku, gdy

• wymiar $N$ układu $Ax=b$ jest "duży", oraz
• macierz $A$ układu jest "rozrzedzona", tzn. ma stosunkowo niewielką liczbę elementów niezerowych, np. proporcjonalną do $N$.

Układy o tych własnościach powstają często przy numerycznym rozwiązywaniu równań różniczkowych cząstkowych.

Metody przestrzeni Kryłowa

Zupełnie inny pomysł na realizację metody iteracyjnej przedstawiają metody przestrzeni Kryłowa, gdzie kolejne przybliżenie $x_k$ dobiera się w taki sposób, by minimalizowało pewną miarę błędu na podprzestrzeni Kryłowa

$K_k=\text{span}\{r,A_r,\dots ,A^{k-1}r\},$

gdzie $r=b-A{x}_0$ jest residuum na początku iteracji. (Zwróć uwagę, że przestrzeń Kryłowa jest rozpięta przez kolejne wektory metody potęgowej] --- to nie przypadek!). W zależności od wyboru sposobu miary błędu, dostajemy inną metodę iteracyjną, takie jak CG, GMRES, PCR, BiCG, i inne. Tutaj omówimy pokrótce tylko najpopularniejszą: CG.

CG

Metoda gradientów sprzężonych, w skrócie CG (conjugate gradients), działa przy założeniu, że $A$ jest symetryczna i dodatnio określona.

Kolejne przybliżenie $x_k$ ma minimalizować błąd w normie energetycznej indukowanej przez $A$,

$||x_k-x|{|}_A=\sqrt{(x_k-x{)}^{T}A(x_k-x)}$

na przestrzeni afinicznej $x_0+K_k$. Okazuje się (co nie jest oczywiste --- trzeba skorzystać z rozmaitych własności ortogonalności generowanych wektorów), że takie zadanie minimalizacji daje się bardzo efektywnie rozwiązać, skąd dostajemy bardzo zwarty algorytm:

Algorytm Metoda CG

---

r = b-A*x;
<math>\displaystyle \rho_0</math> = <math>\displaystyle ||r||_2^2</math>; <math>\displaystyle \beta</math> = 0; k = 1;
while (!stop)
{
	p = r + <math>\displaystyle \beta</math>*p;
	w = A*p;
	<math>\displaystyle \alpha</math> = <math>\displaystyle \rho_{k-1}</math>/<math>\displaystyle p^Tw</math>;
	x = x + <math>\displaystyle \alpha</math>*p;
	r = r - <math>\displaystyle \alpha</math>*w;
	<math>\displaystyle \rho_k</math> = <math>\displaystyle ||r||_2^2</math>;
	<math>\displaystyle \beta</math> = <math>\displaystyle \frac{\rho_{k}}{\rho_{k-1}}</math>;
	k++;	
}

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

Twierdzenie O zbieżności CG jako metody bezpośredniej

Niech $A\in R^{N\times N}$ będzie symetryczna i dodatnio określona. Algorytm CG znajdzie dokładne rozwiązanie po co najwyżej $N$ iteracjach.

Powyższe twierdzenie, choć teoretycznie interesujące, ma małą wartość praktyczną z dwóch powodów:

• dla bardzo dużych $N$, wykonanie $N$ iteracji może być wciąż zbyt kosztownym zadaniem;
• ponieważ w arytmetyce skończonej precyzji ortogonalność --- z której w bardzo istotny sposób korzysta się przy wyprowadzeniu algorytmu --- pogarsza się z iteracji na iterację i w konsekwencji, po wielu iteracjach, jakość $x_k$ przestaje się poprawiać.

Dlatego wygodniej potraktować CG jako metodę iteracyjną. Zachodzi bowiem

Twierdzenie Zbieżność CG jako metody iteracyjnej

Po $k$ iteracjach metody CG,

$||x_k-x|{|}_A\le 2{\left(\frac{\sqrt{\text{cond}(A)}-1}{\sqrt{\text{cond}(A)}+1}\right)}^k||x_0-x|{|}_A,$

gdzie $\text{cond}(A)=||A|{|}_2||A^{-1}|{|}_2$.

GMRES

Metoda GMRES (Generalized Minimum RESidual) nie wymaga ani symetrii, ani dodatniej określoności macierzy, jest więc bardziej uniwersalna, choć też bardziej kosztowna od CG. Jej szczegółowe omówienie, w tym --- oszacowania szybkości zbieżności --- wykracza niestety poza ramy niniejszego wykładu.

Ściskanie macierzy

Zbieżność wszystkich poznanych metod iteracyjnych zależy od własności spektralnych macierzy układu. Pojawiające się w zastosowaniach macierze często mają niekorzystne własności spektralne (np. bardzo duży wskaźnik uwarunkowania), przez co metody iteracyjne zbiegają na nich bardzo wolno.

Dlatego bardzo korzystne może być wstępne przetransformowanie układu

$Ax=b$

z macierzą o niekorzystnych własnościach, do układu

$MAx=Mb,$

gdzie macierz $MA$ ma znacznie korzystniejsze własności z punktu widzenia używanej metody iteracyjnej. W przypadku macierzy symetrycznych widzieliśmy, że kluczowe znaczenie dla zbieżności metody miał rozkład wartości własnych: jeśli wartości własne były bardzo rozrzucone po prostej, to uwarunkowanie było bardzo duże i w konsekwencji zbieżność powolna. Aby zbieżność była szybsza, kluczowe jest, by:

• wartości własne $MA$ były upakowane w klastrach
• najlepiej wszystkie w (małym) otoczeniu wartości 1

Jeśli więc chcielibyśmy przekształcić macierz tak, by metoda iteracyjna dla $MA$ zbiegała szybko, musimy w jakiś sposób "ścisnąć" spektrum macierzy $A$ w okolice jedności. Taką operację nazywamy ściskaniem (preconditioning), a macierz $M$ --- imadłem.

Aby całość miała sens, macierz ściskająca $M$ powinna:

• być łatwa w konstrukcji,
• być tania w mnożeniu przez wektor (głównym elementem każdej metody iteracyjnej jest mnożenie macierzy przez wektor: $M\cdot (A\cdot x)$),
• macierz $MA$ powinna mieć znacznie korzystniejsze własności z punktu widzenia używanej metody iteracyjnej.

Kilka ekstremalnych propozycji na macierz imadła to $M=I$ (łatwa w konstrukcji i tania w mnożeniu, ale niestety nic nie polepsza...) oraz $M=A^{-1}$ (rewelacyjnie poprawia zbieżność metody iteracyjnej, dając zbieżność w jednej iteracji, ale bardzo droga w konstrukcji i mnożeniu). Widać więc, że należy poszukiwać czegoś pośredniego, co niskim kosztem przybliża działanie macierzy odwrotnej.

Dlatego jednym z powszechniej stosowanych (aczkolwiek wciąż nie najbardziej skutecznych) sposobów ściskania są te oparte na zastosowaniu jednego kroku klasycznej metody iteracyjnej.

Zbieżność metody CG bez żadnego ściskania oraz ściśniętej imadłem opartym na jednej iteracji (blokowej) metody Jacobiego.

Inne sposoby ściśnięcia macierzy wykorzystują np. techniki tzw. niepełnego rozkładu macierzy, albo --- w specyficznych przypadkach --- tzw. metody wielosiatkowe.

Okazuje się, że zarówno CG jak i GMRES da się zaimplementować tak, by w jednej iteracji było konieczne tylko jedno mnożenie przez macierz imadła.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 4.6 i 4.7 w

• D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Tym razem wymienione rozdziały nie wystarczą do dogłębnego zrozumienia omawianych treści. Warto więc sięgnąć po jeden z podręczników

• C. T. Kelley, Iterative Solution of Systems of Linear and Nonlinear Equations, SIAM, 1995,
• Y. Saad, Iterative Methods for Sparse Linear Systems, PWS, 1996.