MN08LAB

Układy liniowe z macierzami rzadkimi

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

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

Ć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

• Y. Saad, Iterative methods for sparse linear systems, PWS, 1996.

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^{T}u+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^{T}Ax=A^{T}b,$

którego macierz $\displaystyle A^{T}A$ 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^{T}A$?

Rozwiązanie

Nietrudno sprawdzić, że dla normy spektralnej macierzy,

$\displaystyle {\mbox{cond}}(A^{T}A)=({\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^{T}MAx=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^{T}A$. 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.