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 $Ax=b$ jest metoda Richardsona, zadana wzorem

$x_{k+1}=x_k+\tau (b-A{x}_k)$,

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

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

Testuj na macierzy jednowymiarowego laplasjanu $L$ różnych wymiarów. Jak najefektywniej zaimplementować mnożenie przez $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 $A$ będą $0<{\lambda }_{\min }\le {\lambda }_{\max }$. Ponieważ macierz iteracji metody Richardsona $B=I-\tau A$, to jej wartości własne muszą leżeć w przedziale $[1-\tau {\lambda }_{\max },1-\tau {\lambda }_{\min }]$ i --- oczywiście --- aby iteracja miała sens, $\tau >0$ (dlaczego?).

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

$||x_k-x|{|}_2\le \frac{{\lambda }_{\max }-{\lambda }_{\min }}{{\lambda }_{\max }+{\lambda }_{\min }}||x_{k-1}-x|{|}_2$

(Przy okazji zauważ, że gdyby macierz $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 $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 $A$ przez wektor $x$,
• wyłuskania wartości elementu $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.

$A=\left(\begin{array}{ccccc}{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{array}\right)$

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ę:

$A=\left(\begin{array}{cc}T& v\\ w^T& a_N\end{array}\right)$,

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

$T=\left(\begin{array}{ccccc}{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{array}\right)$,

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

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

$\left(\begin{array}{cc}T& v\\ w^T& a_N\end{array}\right)=\left(\begin{array}{cc}L& \\ l^T& 1\end{array}\right)\cdot \left(\begin{array}{cc}U& u\\ & u_N\end{array}\right)$,

gdzie spełnione są zależności

• $T=LU$ (rozkład LU macierzy trójdiagonalnej $T$)
• $Ul=w$ (rozwiązanie układu równań z macierzą dwudiagonalną)
• $Lu=v$ (jw.)
• $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 $Ax=b$ z macierzą nieosobliwą można transformować do równoważnego mu układu równań normalnych

$A^{T}Ax=A^{T}b$,

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

Rozwiązanie

Nietrudno sprawdzić, że dla normy spektralnej macierzy,

$\text{cond}(A^{T}A)=(\text{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 $M$ mogłoby pomóc, np.

$M{A}^{T}MAx=M{A}^{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 $A^{T}A$. Niedobra metoda to

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

gdyż $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 $A$ (w formacie AIJ) i nie wymaga dodatkowej pamięci.