MN08LAB

Ćwiczenia: układy liniowe z macierzami rzadkimi

Ć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$?

Ć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, bo w nich trzeba zacowywać uporządkowanie. Mnożenie przez wektor najwygodniejsze jest w CSR. Wyłuskanie wartości jest najmniej efektywne w formacie AIJ. Zob. Y. Saad, Iterative methods for sparse linear systems

Ć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ź!

Rozwiązanie

Gdyby pominąć ostatni wiersz i kolumnę, macierz byłaby trójdiagonalna, a z nią już wiemy, co 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ż dobry prekondycjoner mógłby pomóc:

$M{A}^{T}MAx=M{A}^{T}Mb,$

to jednak znacznie lepiej stosować metody opracowane specjalnie dla maicerzy niesymetrycznych, np. GMRES (też z dobrym prekondycjonerem...).

Implementacja metody to tylko decyzja, jak realizować mnożenie przez $A^{T}A$. Oczywiście niedobra metoda to

 
B = A^TA;
...
while ...
	y = B*x;
end

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

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

jest to koszt równy dwukrotnemu mnożeniu przez macierz $A$ (w formacie AIJ).