MN08LAB

From Studia Informatyczne


Układy liniowe z macierzami rzadkimi

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

Oglądaj wskazówki i rozwiązania
Ukryj wskazówki i rozwiązania

Ćwiczenie: Metoda Richardsona

Jedną z najprostszych klasycznych metod iteracyjnych dla równania \displaystyle Ax=b jest metoda Richardsona, zadana wzorem

\displaystyle x_{k+1} = x_k + \tau (b- Ax_k),

gdzie \displaystyle \tau jest pewnym parametrem. Gdy \displaystyle \tau=1, mamy do czynienia ze zwykłą metodą iteracji prostej, która najczęściej nie będzie zbieżna, dlatego wybór parametru \displaystyle \tau jest kluczowy dla skuteczności metody.

Dla \displaystyle A symetrycznej, dodatnio określonej sprawdź, przy jakich założeniach o \displaystyle \tau metoda będzie zbieżna do rozwiązania \displaystyle x^* z dowolnego wektora startowego \displaystyle x_0 i oceń szybkość tej zbieżności.

Testuj na macierzy jednowymiarowego laplasjanu \displaystyle L różnych wymiarów. Jak najefektywniej zaimplementować mnożenie przez \displaystyle L?

Wskazówka

Zadanie można rozwiązać przez analogię do opisywanego w wykładzie przypadku zastosowania metody Jacobiego do macierzy jednowymiarowego laplasjanu: wystarczy uogólnić!

Rozwiązanie

Niech ekstremalnymi wartościami własnymi macierzy \displaystyle A będą \displaystyle 0 < \lambda_{\min} \leq \lambda_{\max}. Ponieważ macierz iteracji metody Richardsona \displaystyle B= I - \tau A, to jej wartości własne muszą leżeć w przedziale \displaystyle [1-\tau \lambda_{\max}, 1 - \tau \lambda_{\min}] i --- oczywiście --- aby iteracja miała sens, \displaystyle \tau > 0 (dlaczego?).

Co więcej, \displaystyle ||B||_2 < 1 wtedy i tylko wtedy, gdy \displaystyle 0 < \tau < \frac{2}{\lambda_{\max}}, a najmniejszą normę spektralną macierzy \displaystyle B uzyskamy, gdy \displaystyle \tau = \frac{2}{\lambda_{\max} + \lambda_{\min}} i wówczas

\displaystyle ||x_k-x||_2 \leq \frac{\lambda_{\max} - \lambda_{\min}}{\lambda_{\max} + \lambda_{\min}}||x_{k-1}-x||_2.

(Przy okazji zauważ, że gdyby macierz \displaystyle A nie była określona, tzn. miałaby zarówno dodatnie, jak i ujemne wartości własne, to metoda Richardsona mogłaby w ogóle nie być zbieżna (dla pewnych wektorów startowych))

Jeśli chodzi o mnożenie przez macierz jednowymiarowego laplasjanu, to najprościej wcale nie używać struktury macierzy! Rzeczywiście, jedyne, czego potrzebujemy to operacja mnożenia wektora przez macierz \displaystyle L, a to realizuje pętla:

void LapMult(double *x, int N, double *y)
/* 
Mnożenie wektora przez macierz 1-wym laplasjanu L wymiaru N

	Wejście:
	x - wektor, który mnożymy przez L
	N - jego długość
	
	Wyjście:
	y - wynik: y = Lx
*/
{
	int i;

	for (i = 1; i < N-1; i++)
		y[i] = -x[i-1] + 2.0*x[i] - x[i+1];

	y[0] = 2.0*x[0] - x[1];
	y[N-1] = -x[N-2] + 2.0*x[N-1];
}

Ćwiczenie

Zaimplementuj operacje:

  • mnożenia macierzy \displaystyle A przez wektor \displaystyle x,
  • wyłuskania wartości elementu \displaystyle A_{ij},
  • zmiany wartości pewnego zerowego wyrazu macierzy na niezerową,

jeśli macierz jest zadana w formacie

  • AIJ,
  • CSC,
  • CSR.

Przetestuj dla kilku macierzy z kolekcji MatrixMarket.

Rozwiązanie

Dodawanie nowego elementu w formacie AIJ jest bardzo łatwe. W przeciwieństwie do pozostałych, w których trzeba zachowywać uporządkowanie. Mnożenie przez wektor najwygodniejsze jest w CSR, bo dodatkowo narzuca zasadę lokalności w przestrzeni. Wyłuskanie wartości jest najmniej efektywne w formacie AIJ. Szczegóły opisane są w rozdziale 3.5 książki

Zobacz także implementacje w Fortranie, w pakiecie SPARSKIT, będącym czymś w rodzaju odpowiednika BLAS dla macierzy rozrzedzonych.

Ćwiczenie: Konwersja formatu macierzy rzadkiej

Napisz procedurę aij2csr, konwertującą macierz w formacie AIJ do CSR i csr2aij, działającą w drugą stronę.

Rozwiązanie

Zobacz, jak to zrobiono w pakiecie SPARSKIT.


Ćwiczenie

Jak tanio rozwiązywać układ równań z macierzą cykliczną trójdiagonalną, tzn.

\displaystyle  A =  \begin{pmatrix}   a_1 & c_1 &  &  & b_1\\  b_2 & a_2 & c_2       &   & \\   &  b_3 & a_3  & \ddots &  \\   & & \ddots & \ddots  & c_{N-1}\\ c_N  &   &      & b_N  & a_N \end{pmatrix}

Dla uproszczenia załóżmy, że macierz jest dodatnio określona i symetryczna.

Zaimplementuj opracowaną metodę, korzystając z BLASów i LAPACKa.

Wskazówka

Dziel i rządź! (Tylko: jak?)

Rozwiązanie

Gdyby pominąć ostatni wiersz i kolumnę, macierz byłaby trójdiagonalna, a my już wiemy, co z nią zrobić... Jak więc sprytnie pozbyć się ostatniej niewiadomej i ostatniego równania? W naszej macierzy wyróżnijmy ostatni wiersz i kolumnę:

\displaystyle A =  \begin{pmatrix}  T & v\\ w^T & a_N \end{pmatrix} ,

gdzie \displaystyle T jest \displaystyle N-1 podmacierzą główną \displaystyle A,

\displaystyle T =  \begin{pmatrix}   a_1 & c_1 &  &  & \\  b_2 & a_2 & c_2       &   & \\   &  b_3 & a_3  & \ddots &  \\   & & \ddots & \ddots  & c_{N-2}\\   &   &      & b_{N-1}  & a_{N-1} \end{pmatrix} ,

natomiast \displaystyle w^T = [c_N, 0, \ldots, 0, b_N], \displaystyle v = [b_1, 0, \ldots, 0, c_{N-1}]^T.

Mając rozkład \displaystyle T=LU, łatwo stąd wygenerować rozkład \displaystyle A, gdyż

\displaystyle \begin{pmatrix}  T & v\\ w^T & a_N \end{pmatrix}  =  \begin{pmatrix}  L & \\ l^T & 1 \end{pmatrix}  \cdot \begin{pmatrix}  U & u\\   & u_N \end{pmatrix} ,

gdzie spełnione są zależności

  • \displaystyle T = LU (rozkład LU macierzy trójdiagonalnej \displaystyle T)
  • \displaystyle Ul = w (rozwiązanie układu równań z macierzą dwudiagonalną)
  • \displaystyle Lu = v (jw.)
  • \displaystyle l^Tu + u_N = a_N.

Ćwiczenie: CGNE

Ktoś mógłby sugerować, że skoro CG działa tylko dla macierzy symetrycznych, to dowolny układ \displaystyle Ax=b z macierzą nieosobliwą można transformować do równoważnego mu układu równań normalnych

\displaystyle A^TAx = A^T b,

którego macierz \displaystyle A^TA jest już oczywiście macierzą symetryczną i dodatnio określoną.

Wskaż potencjalne wady tej metody i podaj sposób jej implementacji.

Wskazówka

Jakie jest uwarunkowanie macierzy \displaystyle A^TA?

Rozwiązanie

Nietrudno sprawdzić, że dla normy spektralnej macierzy,

\displaystyle  \mbox{cond} (A^TA) = ( \mbox{cond} (A))^2,

a więc w przypadku macierzy źle uwarunkowanych należy spodziewać się patologicznie dużej liczby iteracji. Chociaż dobre (symetryczne, dodatnio określone) imadło \displaystyle M mogłoby pomóc, np.

\displaystyle MA^TMAx = MA^T Mb,

to jednak znacznie lepiej stosować metody opracowane specjalnie dla macierzy niesymetrycznych, np. GMRES (oczywiście z nieodzownym ściskaniem macierzy, gdy jest źle uwarunkowana...).

Implementacja metody iteracyjnej to tylko decyzja, jak realizować mnożenie przez \displaystyle A^TA. Niedobra metoda to

B = <math>\displaystyle A^TA</math>;
...
while ...
	y = B*x;
end

gdyż \displaystyle B będzie bardziej wypełniona niż A. Znacznie lepiej

...
while ...
	y = A*x;
	y = (y'*A)';
end

co realizuje się kosztem równym dwukrotnemu mnożeniu przez macierz \displaystyle A (w formacie AIJ) i nie wymaga dodatkowej pamięci.