Metody Numeryczne: 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. 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ą $\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!
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 $\displaystyle 0.006$ , a więc macierz jest bardzo rozrzedzona: w każdym wierszu są średnio tylko 44 niezerowe elementy.

[Edytuj]

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

[Edytuj]

Format współrzędnych

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

$\displaystyle A_{I[k], J[k]} = V[k], \qquad\forall k = 0,\ldots, 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.

Format diagonalny

Znacznie mniej uniwersalny niż poprzednie i dlatego rzadziej spotykany. Kolejne diagonale macierzy przechowujemy w kolejnych wierszach macierzy $\displaystyle P\times N$ , gdzie $\displaystyle P$ jest liczbą niezerowych diagonal w $\displaystyle 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.

Macierze specjalne

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

$\displaystyle Ax = b,$

ale w sytuacji, gdy macierz $\displaystyle 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 $\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.

Macierze trójdiagonalne

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

$\displaystyle T = \begin{pmatrix} 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{pmatrix}$

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

$\displaystyle T x = 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 $\displaystyle T$ ma słabo dominującą diagonalę, tzn.

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

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

Metody Numeryczne

poniedziałek, 17 grudnia 2012

Wielkie układy równań liniowych