Metody Numeryczne: 2012

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.

Podstawowe pojęcia

Weryfikacja hipotez parametrycznych

Przykład

Należy dokonać oceny partii pudelek zapalek liczącej 100 tys. sztuk. dostawca twierdzi, że w pudelku znajdują się przecietnie 54 zapalki. Zweryfikować hipotezę

H 0 (m = m 0 = 54)

. Ponieważ nie znamy rozkładu liczby zapałek w pudełkach w populacji generalnej, a mozemy łatwo pobrać próbę >= 30 możemy wieć w przybliżeniu skorzystać z rozkładu normalnego. Zkaldamy, że przy próbie o wielkości n = 100 odnotowano średnią arytmetyczną

m n = 51,21

natomiast

σ n' = 2,54

. Weryfikujemy przy poziomie istotności

α = 0,02

ponieważ obraliśmy dużą próbę wiec $\sigma_{M_n} \approx {{\sigma_n'} \over { \sqrt{n}}} = 0{,}245$ . Musimy zatem wyznaczyć t dla ktorego $P \left ( {{|M_n-m_0|} \over {\sigma_{M_n}}} \ge t \right ) = 0{,}02$
Definiujemy unormowaną zmienną Y:
$Y={{M_n-m_0} \over {\sigma_{M_n}}}$
podstawiamy do wzoru
$P(|Y| \ge t) = 0{,}02$
Z własności bezwzględnej wartości:
$P[(Y \ge t) \vee (Y \le -t)] = 0{,}02$
$P(Y \ge t)+P(Y \le -t) = 0{,}02$
Ponieważ funkcja gęstości jest dla rozkładu N(0,1) parzysta to zachodzi równość:
$P(Y \ge t) = P(Y \le -t)$
$2P(Y \ge t) = 0{,}02$
$P(Y \ge t) = 0{,}01$
Wiadomo, że P(A) to to samo co 1-P(A') więc:

1 - P (Y < t) = 0,01

P (Y < t) = 0,99

A P(Y<x) to dystrybuanta - czyli F(x):

F (t) = 0,99

Teraz w tablicy rozkładu normalnego znajdujemy najmniejszą wartość t dla której F(t) wynosi conajmniej 0,99. Jest to wartość 2,33.
Hipotezę

H 0

należy wiec odrzucić na poziomie istotności

α

jeżeli $|M_n - 54| \ge t \cdot 0{,}245$ , w przeciwnym wypadku przy zadanej istotności

α = 0,02

nie możemy ani potwierdzić hipotezy, ani jej odrzucić.
Zgodnie z naszymi danymi wychodzi:

| M n - 54 | = | 51,21 - 54 | = 2,79

$\sigma_{M_n} \cdot t = 0{,}245 \cdot 2{,}33 = 0{,}57085$
Więc:
$| M_n - 54 | \ge \sigma_{M_n} \cdot t$
$2{,}79 \ge 0{,}57085$
Zatem hipotezę możemy odrzucić (jeśli wyjdzie odwrotnie to piszemy że nie odrzucamy ani nie potwierdzamy - tak właśnie trzeba było zrobić na egzaminie, bo wychodziło <).

Metody Numeryczne

poniedziałek, 17 grudnia 2012

Wielkie układy równań liniowych