Metody numeryczne i probabilstyczne: Wielkie układy równań liniowych

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 $\displaystyle N$ ma tylko $\displaystyle O(N)$ niezerowych elementów. Wykorzystanie 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ą $\displaystyle 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!

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

Macierz rozrzedzona (Octave):
<realnowiki><realnowiki><realnowiki><realnowiki>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); </realnowiki></realnowiki></realnowiki></realnowiki>

(http://osilek.mimuw.edu.pl/index.php?title=MN08)

Macierze taśmowe

Macierz $\displaystyle A$ taka, że dla $\displaystyle |i-j| \geq d$ , $\displaystyle a_{ij} = 0$ , nazywamy macierzą taśmową z rozstępem $\displaystyle d$ , o szerokości pasma $\displaystyle 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 $\displaystyle \max\{2d,N\}$ --- tak więc, dla niezbyt dużych $\displaystyle d$ wciąż wynikowa macierz jest taśmowa.

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

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

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ą", $\displaystyle M$ , i "resztę", $\displaystyle Z$ . Dokładniej, jeśli $\displaystyle M$ jest nieosobliwa, to równanie $\displaystyle Ax = b$ można zapisać jako zadanie punktu stałego

$\displaystyle Mx = Zx + b,$

gdzie $\displaystyle Z=M-A$ . Inaczej:

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

i zastosować doń metodę iteracji prostej Banacha:

$\displaystyle x_{k+1} = M^{-1}(Zx_k + b).$

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

$\displaystyle x_k\,=\,B x_{k-1}\,+\, c,$

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

W tym przypadku

$\displaystyle x_k- x^*\,=\,B^k( x_0- x^*),$

a stąd i z nierówności $\displaystyle \|B^k\|\le\|B\|^k$ , mamy

$\displaystyle \| x_k- x^*\|\,\le\, \|B\|^k\| x_0- x^*\|.$

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

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

Tak więc, metody oparte na iteracji prostej będą zbieżne liniowo z ilorazem $\displaystyle \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.

Metody numeryczne i probabilstyczne

poniedziałek, 17 grudnia 2012

Wielkie układy równań liniowych

Brak komentarzy:

Prześlij komentarz