MN08

Wielkie układy równań liniowych

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 $N$ ma tylko $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ą $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.).

Przykład: Macierz z kolekcji Boeinga

Macierz sztywności dla modelu silnika lotniczego wygenerowana swego czasu w zakładach Boeinga pochodzi z [http://www.cise.ufl.edu/research/sparse/matrices/Boeing kolekcji Tima Davisa] (autora bardzo dobrego solvera równań liniowych z macierzami rzadkimi). Jest to mała macierz, wymiaru 8032 (w kolekcji spotkasz równania z milionem i więcej niewiadomych).

Struktura niezerowych elementów macierzy bcsstk38.

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 MBEACXC

Dane pochodzą z serwisu MatrixMarket. Jest to niezbyt duża macierz, wymiaru $496\times 496$, niesymetryczna, o elementach rzeczywistych. Źródłem tej macierzy jest pewien model ekonomiczny.

Struktura niezerowych elementów macierzy MBEACXC

Tylko pozornie macierz ta wydaje się dość gęsta, w rzeczywistości jej współczynnik wypełnienia wynosi około 20 procent.

Formaty 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!).

Format współrzędnych (AIJ)

Do zapamiętania macierzy $A$ wymiaru $N\times M$ o $NNZ$ niezerowych elementów wykorzystujemy trzy wektory: I, J --- oba typu int --- oraz V, typu double, wszystkie o długości $NNZ$, przy czym

$A_{I[k],J[k]}=V[k],\mathrm{\forall }k=0,\dots ,NNZ-1.$

Format AIJ

W tym formacie wprowadzane są macierze rzadkie do Octave'a i MATLABa:

 
A = sparse(I,J,V,N,N);

Format spakowanych kolumn (wierszy)

Format współrzędnych nie narzucał żadnego uporządkowania elementów macierzy --- można było je umieszczać w dowolnej kolejności. Narzucenie sensownego porządku mogłoby wspomóc realizację wybranych istotnych operacji na macierzy, na przykład, aby wygodnie było realizować działanie (prawostronnego) mnożenia macierzy przez wektor, wygodnie byłoby przechowywać elementy macierzy wierszami. Tak właśnie jest zorganizowany format spakowanych wierszy, Compressed Sparse Row (CSR) . Analogicznie jest zdefiniowany format spakowanych kolumn, Compressed Sparse Column (CSC) , którym zajmiemy się bliżej.

Podobnie jak w przypadku formatu współrzędnych, macierz w formacie CSC jest przechowywana w postaci trzech wektorów: AV jest wektorem typu double o długości $NNZ$, zawierającym kolejne niezerowe elementy macierzy $A$ wpisywane kolumnami, AI jest wektorem typu int o długości $NNZ$, zawierającym numery wierszy macierzy $A$ odpowiadających elementom z AV. Natomiast zamiast tablicy J, jak to było w formacie współrzędnych, mamy krótszy wektor typu int, AP, o długości $M$ zawierający na $j$-tym miejscu indeks pozycji w AV, od którego rozpoczynają się w AV elementy $j$-tej kolumny macierzy $A$.

Format CSC

Mamy więc zależność, przy założeniu, że $AP[0]=0$,

$A_{AI[AP[j]+k],j+1}=AV[AP[j]+k],j=0,\dots ,M-1,k=0,\dots ,AP[j+1]-1.$

Taki (z drobnymi modyfikacjami) format macierzy wykorzystują np. pakiety ARPACK i UMFPACK.

Format diagonalny

Znacznie mniej uniwersalny niż poprzednie i dlatego rzadziej spotykany. Kolejne diagonale macierzy przechowujemy w kolejnych wierszach macierzy $P\times N$, gdzie $P$ jest liczbą niezerowych diagonal w $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

$Ax=b,$

ale w sytuacji, gdy macierz $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 $A$ taka, że dla $|i-j|\ge d$, $a_{ij}=0$, nazywamy macierzą taśmową z rozstępem $d$, o szerokości pasma $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 $\max \{2d,N\}$ --- tak więc, dla niezbyt dużych $d$ wciąż wynikowa macierz jest taśmowa.

W szczególności, gdy macierz jest taśmowa z pasmem o rozstępie $d$ i jednocześnie diagonalnie dominująca, wtedy rozkład LU takiej macierzy da się wykonać w miejscu kosztem $O(dN)$ czyli liniowym.

W LAPACKu zaimplementowano szybki solver równań z macierzami taśmowymi, DGBSV wymagający skądinąd specjalnego sposobu przechowywania macierzy wykorzystyjącego format diagonalny.

Macierze trójdiagonalne

Szczególnym przypadkiem macierzy taśmowych są macierze trójdiagonalne, tzn. taśmowe o rozstępie $d=1$:

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

Zadanie rozwiązywania równań z taką macierzą

$Ax=e$

jest tak często spotykane, że warto przytoczyć algorytm w całej okazałości --- popularnie zwany algorytmem przeganiania.

Zacznijmy jednak od stwierdzenia, kiedy taka macierz nie wymaga wyboru elementu głównego:

Stwierdzenie

Jeśli macierz $A$ ma słabo dominującą przekątną, tzn.

$|a_i|\ge |b_i|+|c_i|,1\le i\le N,$

($b_1=0=c_N$) i przynajmniej dla jednego indeksu "$i$" mamy powyżej ostrą nierówność "$>$", to algorytm przeganiania jest wykonalny bez przestawień wierszy. Ponadto wymaga on $7N$ operacji arytmetycznych, a więc jest prawie optymalny.

Algorytm Metoda przeganiania

---

<math>\displaystyle d_1</math> = <math>\displaystyle a_1</math>; 
<math>\displaystyle f_1</math> = <math>\displaystyle e_1</math>;
for (i = 2; i <= N; i++)
{
	<math>\displaystyle l</math> = <math>\displaystyle b_i</math>/<math>\displaystyle a_{i-1}</math>;
	<math>\displaystyle d_i</math> = <math>\displaystyle a_i</math> - <math>\displaystyle l</math> * <math>\displaystyle c_{i-1}</math>;
	<math>\displaystyle f_i</math> = <math>\displaystyle e_i</math> - <math>\displaystyle l</math> * <math>\displaystyle f_{i-1}</math>;
} 
<math>\displaystyle x_1</math> = <math>\displaystyle f_N</math>;
for (i = N-1; i >= 1; i--)
	<math>\displaystyle x_i</math> = <math>\displaystyle f_i</math> - <math>\displaystyle c_i</math> * <math>\displaystyle x_{i+1}</math>;

Metody bezpośrednie

Przykład: Strzałka Wilkinsona

Rozważmy układ równań z macierzą diagonalnie dominującą postaci

$A=\left(\begin{array}{ccccc}*& *& *& \cdots & *\\ *& *& & & \\ *& & *& & \\ ⋮& & & \ddots & \\ *& & & & *\end{array}\right)$

gdzie $*$ oznacza jakiś niezerowy element. Łatwo sprawdzić, że chociaż wyjściowa macierz jest rozrzedzona, to zastosowanie do niej eliminacji Gaussa powoduje, że w wyniku dostajemy gęste czynniki rozkładu.

Tymczasem wystarczy odwrócić kolejność równań i numerację niewiadomych (co dla macierzy jest równoznaczne z odwróceniem porządku wierszy i kolumn, korzystając z pewnej (jakiej?) macierzy permutacji $P$):

$\tilde{A}=PA{P}^T=\left(\begin{array}{ccccc}*& & & & *\\ & \ddots & & & ⋮\\ & & *& & *\\ & & & *& *\\ *& \cdots & *& *& *\end{array}\right)$

Wtedy okazuje się, że rozkład naszej macierzy nie powoduje już wypełnienia fill-in czynników rozkładu!

Właśnie na tym polega główny problem w rozwiązywaniu układów z macierzami rzadkimi metodami bezpośrednimi: maksymalnie wykorzystać rozrzedzenie macierzy tak, by czynniki rozkładu były możliwie mało wypełnione. Albowiem wiedząc to będziemy mogli ograniczyć się jedynie do fizycznego wyznaczenia wartości niezerowych elementów macierzy rozkładu. Ponadto wymagania pamięciowe algorytmu nie będą istotnie wykraczać ponad ilość pamięci potrzebnej na przechowanie danych (macierzy).

W ogólnym przypadku rozwiązanie takiego zadania jest trudne i większość algorytmów opiera się na pewnych heurystykach, które warto wspomóc wcześniejszą analizą konkretnego układu równań, który mamy rozwiązać. Najczęściej dąży się do takiego przenumerowania równań i niewiadomych, by w efekcie z góry przewidzieć, gdzie wystąpią zera w macierzach rozkładu --- i, by takich zer było jak najwięcej (by wypełnienie było jak najmniejsze)! Na architekturach z pamięcią hierarchiczną dąży się także do tego, by w trakcie rozkładu można było korzystać z BLAS3, a więc permutacje wierszy i kolumn macierzy muszą to także brać pod uwagę (tzw. metody wielofrontowe).

Stosuje się kilka strategii wyznaczania korzystnych permutacji reorderingu , z których warto wymienić

• przybliżone algorytmy minimalnego stopnia approximate minimum degree (AMD)
• techniki

podziału grafów na (prawie) rozłączne składowe (METIS).

Używają ich biblioteki takie jak UMFPACK, czy HSL.

W Octave mamy do dyspozycji także kilka procedur generujących takie permutacje, w tym: colamd (AMD dla macierzy niesymetrycznych) oraz symamd (AMD dla macierzy symetrycznych). Większy wybór oferuje MATLAB, jednak należy bezwzględnie pamiętać o jednym: nie ma uniwersalnej metody reorderingu i dla konkretnej macierzy może istnieć specjalna metoda, która da oszałamiające rezultaty, podczas gdy standardowe podejścia nie dadzą efektu.

Przykład: Wypełnienie pewnej macierzy w zależności od użytego algorytmu

Rozważmy macierz pochodzącą z kolekcji Tima Davisa, dotyczącą pewnego zadania mechaniki strukturalnej jakie pojawiło się u Boeinga. Jest to macierz symetryczna, dodatnio określona, wymiaru 8032.

Struktura niezerowych elementów macierzy.

Jest to macierz silnie rozrzedzona, ponieważ ma tylko 355460 niezerowych elementów (czyli wypełnienie to tylko 0.5 procent).

Zobaczymy, jak w zależności od użytego algorytmu permutacji kolumn i wierszy poradzi sobie algorytm rozkładu Cholesky'ego.

Czynnik rozkładu Cholesky'ego $L$ wykonanego standardowym algorytmem. Czas rozkładu: 0.892013

Czynnik rozkładu Cholesky'ego $L$ z reorderingiem COLAMD. Czas rozkładu: 0.813038

Czynnik rozkładu Cholesky'ego $L$ z reorderingiem SYMAMD. Czas rozkładu: 0.487683s. Prawie dwa razy szybciej niż bez reorderingu, chociaż i tak wskutek wzrostu wypełnienia macierzy w dolnym trójkącie mamy aż 2 procent niezerowych elementów.

Czynnik rozkładu Cholesky'ego $L$ z odwrotnym reorderingiem Cuthill-McKee. Czas rozkładu: 0.845928

Czynnik rozkładu Cholesky'ego $L$ z jeszcze innym (bardzo tanim, ale jak widać czasem zupełnie złym) reorderingiem. Czas rozkładu: 5.947936

Na zakończenie popatrzmy, jak ważne jest spostrzeżenie symetrii macierzy:

Rozkład LU, czynnik $L$ (bez reorderingu). Czas rozkładu LU: 1.696018s. Nieakceptowalny, podobnie jak drastyczne wypełnienie macierzy.

Rozkład LU, czynnik $U$ (bez reorderingu)

Jak widać, w naszym przypadku standardowe algorytmy (COLAMD i SYMAMD) poradziły sobie całkiem nieźle, chociaż wypełnienie i tak znacząco wzrosło. Zapewne algorytm oparty na podziale grafu na prawie rozłączne składowe mógłby tu jeszcze lepiej zadziałać.

Stacjonarne metody iteracyjne

Gdy macierz jest rozrzedzona, mnożenie takiej macierzy przez wektor jest tanie (koszt jest proporcjonalny do liczby niezerowych elementów macierzy). Dlatego, jeśli możemy zadowolić się rozwiązaniem przybliżonym układu, a w zamian osiągnąć je tanim kosztem, warto rozważyć metody iteracyjne.

Najprostsze metody iteracyjne (najprostsze w analizie i implementacji, ale --- jak można się domyślić --- w praktyce najmniej efektywne) polegają na rozkładzie macierzy na część "łatwo odwracalną", $M$, i "resztę", $Z$. Dokładniej, jeśli $M$ jest nieosobliwa, to równanie $Ax=b$ można zapisać jako zadanie punktu stałego

$Mx=Nx+b,$

gdzie $Z=M-A$. Inaczej:

$x=M^{-1}(Zx+b),$

i zastosować doń Dodaj link: metodę iteracji prostej Banacha:

$x_{k+1}=M^{-1}(Z{x}_k+b).$

Takie metody nazywamy stacjonarnymi metodami iteracyjnymi.

Aby przeanalizować zbieżność takiej metody, warto rozpatrzyć przypadek ogólniejszy

$x_k=B{x}_{k-1}+c,$

dla pewnej macierzy $B$ wektora $c$ (dla stacjonarnej metody iteracyjnej, $B=M^{-1}Z$ oraz $c=M^{-1}b$).

W tym przypadku

$x_k-x^{*}=B^k(x_0-x^{*}),$

a stąd i z nierówności $\Vert B^k\Vert \le \Vert B{\Vert }^k$, mamy

$\Vert x_k-x^{*}\Vert \le \Vert B{\Vert }^k\Vert x_0-x^{*}\Vert .$

Warunkiem dostatecznym zbieżności iteracji prostych jest więc $\Vert B\Vert <1$. Okazuje się, że warunkiem koniecznym i dostatecznym zbieżności tej iteracji dla dowolnego wektora startowego $x_0$ jest

$\rho =\max \{|\lambda |:\lambda \text{ jest wartością własną }A\}<1.$

Tak więc, metody oparte na iteracji prostej będą zbieżne liniowo z ilorazem $\rho$.

Metoda Jacobiego

Plik:Jacobi.jpg

Jacobi
Zobacz biografię

Biorąc $D=diag(A)$, gdzie $diag(A)$ jest macierzą diagonalną składającą się z wyrazów stojących na głównej przekątnej macierzy $A$, układ $Ax=b$ jest równoważny układowi

$Dx=Zx+b,$

a stąd (o ile na przekątnej macierzy $A$ nie mamy zera) otrzymujemy metodę iteracyjną

$x_k=B{x}_{k-1}+c,$

,

gdzie $B=-D^{-1}Z$ i $c=D^{-1}b$, zwaną metodą Jacobiego.

W metodzie Jacobiego warunek dostateczny zbieżności, $\Vert B\Vert <1$, jest spełniony np. wtedy, gdy macierz $A$ ma dominującą przekątną, tzn. gdy

$|a_{i,i}|>\sum _{j\ne i}|a_{i,j}|,\mathrm{\forall }i=1,\dots ,N.$

Rzeczywiście, ponieważ wyraz $(i,j)$ macierzy $D^{-1}Z$ wynosi $0$ dla $i=j$ oraz $a_{i,j}/a_{i,i}$ dla $i\ne j$, a więc

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \lefteqn{\|D^{-1}Z\|_\infty \;=\; \max_{1\le i\le N} \sum_{j=1,j\ne i}^N{|a_{i,j}|}/{|a_{i,i}|} } \\ && =\;\max_{1\le i\le N} \sum_{j=1}^N{|a_{i,j}|}/{|a_{i,i}|}\,-\,1\;<\;1, \endaligned}

przy czym ostatnia nierówność wynika z warunku diagonalnej dominacji.

Przykład: Macierz laplasjanu

Macierz $N\times N$, zwana macierzą jednowymiarowego laplasjanu

$L=\left(\begin{array}{cccc}2& -1& & \\ -1& 2& \ddots & \\ & \ddots & \ddots & -1\\ & & -1& 2\end{array}\right)$

pojawia się w bardzo wielu zastosowaniach, także jako podzadanie w algorytmach numerycznych. Ta macierz jest macierzą taśmową, symetryczną i dodatnio określoną, więc rozwiązanie układu równań z tą macierzą metodami bezpośrednimi jest łatwe, kosztem $O(N)$. Ciekawe będzie, jak poradzą sobie z nią metody iteracyjne.

Na przykład, w metodzie Jacobiego weźmiemy $D=2I$ oraz $Z=L-D$. Obliczając normę macierzy iteracji Jacobiego dostajemy $||D^{-1}Z|{|}_{\mathrm{\infty }}=1$, co nie rozstrzyga jeszcze o jej (nie)zbieżności.

Okazuje się, że są gotowe wzory na wartości własne macierzy $L$:

${\lambda }_j=4{\sin }^2\left(\frac{j\pi }{2(N+1)}\right),$

dla $1\le j\le N.$

i w konsekwencji, wartościami własnymi $D^{-1}Z=\frac{1}{2}Z=\frac{1}{2}L-I$ są liczby ${\mu }_i=\frac{1}{2}{\lambda }_i-1$. Ponieważ $0<{\mu }_i<1$, oraz $||Z|{|}_2=1-O(N^{-2})<0$, to metoda, choć zbieżna, dla dużych $N$ staje się zbieżna tak wolno, że w praktyce bezużyteczna.

Zaletą stacjonarnych metod iteracyjnych jest również ich prostota, przez co są one łatwe do zaprogramowania.

Metoda Gaussa--Seidela

Ciekawostką jest, że Gauss nie miał z nią nic wspólnego, a Seidel był wręcz jej przeciwnikiem...

Złożoność stacjonarnych metod iteracyjnych

Zastanówmy się teraz nad złożonością metod iteracyjnych. Ponieważ możemy jedynie znaleźć pewne przybliżenie rozwiązania dokładnego $x^{*}$, przez złożoność metody będziemy rozumieli koszt kombinatoryczny obliczenia $x_k$ z zadaną dokładnością $ϵ>0$. Dla uproszczenia założymy, że medoda jest zbieżna liniowo z ilorazem $\rho$. Zauważmy, że aby zredukować błąd początkowy do $ϵ>0$, wystarczy wykonać $k=k(ϵ)$ iteracji, gdzie $k$ spełnia

${\rho }^k\Vert x_0-x^{*}\Vert \le ϵ,$

czyli

$k\ge \frac{\log (1/ϵ)-\log (1/\Vert x_0-x^{*}\Vert )}{\log (1/\rho )}.$

Liczba ta zależy więc w istotny sposób od błędu początkowego i (przede wszystkim) od współczynnika redukcji błędu $\rho$, natomiast zależność od dokładności $ϵ$ i wymiaru $N$ układu jest dużo mniej istotna (w zadaniach praktycznych -- takich jak jednowymiarowy laplasjan --- jednak często okazuje się, że... $\rho$ zależy od $N$!).

Zakładając, że koszt jednej iteracji wynosi $c=c(N)$ ($c(N)$ jest tym mniejszy, im mniejsza jest liczba niezerowych elementów macierzy $A$), złożoność metody jest proporcjonalna do

$c(N)\frac{\log (1/ϵ)}{\log (1/\rho )}.$

Stąd oczywisty wniosek, że metody iteracyjne warto stosować zamiast metod bezpośrednich w przypadku, gdy

• wymiar $N$ układu $Ax=b$ jest "duży", oraz
• macierz $A$ układu jest "rozrzedzona", tzn. ma stosunkowo niewielką liczbę elementów niezerowych, np. proporcjonalną do $N$.

Układy o tych własnościach powstają często przy numerycznym rozwiązywaniu równań różniczkowych cząstkowych.

Metody przestrzeni Kryłowa

Zupełnie inny pomysł na realizację metody iteracyjnej przedstawiają metody przestrzeni Kryłowa, gdzie kolejne przybliżenie $x_k$ dobiera się w taki sposób, by minimalizowało pewną miarę błędu na podprzestrzeni Kryłowa

$K_k=\text{span}\{r,A_r,\dots ,A^{k-1}r\},$

gdzie $r=b-A{x}_0$ jest residuum na początku iteracji. W zależności od wyboru sposobu miary błędu, dostajemy inną metodę iteracyjną, takie jak CG, GMRES, PCR, BiCG, i inne. Tutaj omówimy pokrótce tylko dwie najpopularniejsze: CG i GMRES.

CG

Metoda gradientów sprzężonych conjugate gradients (CG) działa przy założeniu, że $A$ jest symetryczna i dodatnio określona.

Kolejne przybliżenie $x_k$ ma minimalizować błąd w normie energetycznej indukowanej przez $A$,

$||x_k-x|{|}_A=\sqrt{(x_k-x{)}^{T}A(x_k-x)}$

na przestrzeni afinicznej $x_0+K_k$. Okazuje się (co nie jest oczywiste --- trzeba skorzystać z rozmaitych własności ortogonalności generowanych wektorów), że takie zadanie minimalizacji daje się bardzo efektywnie rozwiązać, skąd dostajemy bardzo zwarty algorytm:

Algorytm

---

r = b-A*x;
<math>\displaystyle \rho_{-1}</math> = 0; <math>\displaystyle \rho_0</math> = <math>\displaystyle ||r||_2^2</math>; <math>\displaystyle \beta</math> = 0; k = 1;
while <math>\displaystyle \sim</math>STOP
{
	p = r + <math>\displaystyle \beta</math>*p;
	w = A*p;
	<math>\displaystyle \alpha</math> = <math>\displaystyle \frac{\rho_{k-1}}{p^Tw}</math>;
	x = x + <math>\displaystyle \alpha</math>*p;
	r = r - <math>\displaystyle \alpha</math>*w;
	<math>\displaystyle \rho_k</math> = <math>\displaystyle ||r||_2^2</math>;
	<math>\displaystyle \beta</math> = <math>\displaystyle \frac{\rho_{k}}{\rho_{k-1}}</math>;
	k++;	
}

Jak widać, całą iterację da się wykonać przechowując w pamięci tylko kilka wektorów (a nie, jak możnaby się obawiać, całą przestrzeń $K_k$), a najdroższym jej elementem jest mnożenie macierzy przez wektor.

Twierdzenie Zbieżność CG jako metody bezpośredniej

Niech $A\in R^{N\times N}$ będzie symetryczna i dodatnio określona. Algorytm CG znajdzie dokładne rozwiązanie po co najwyżej $N$ iteracjach.

Powyższe twierdzenie, choć teoretycznie interesujące, ma małą wartość praktyczną z dwóch powodów:

• dla bardzo dużych $N$, wykonanie $N$ iteracji może być wciąż zbyt kosztownym zadaniem;
• ponieważ w arytmetyce $f{l}_{\nu }$ skończonej precyzji ortogonalność, z której korzysta się przy wyprowadzeniu algorytmu nie jest zachowana, i w konsekwencji, po wielu iteracjach, jakość $x_k$ przestaje się poprawiać.

Dlatego wygodniej potraktować CG jako metodę iteracyjną. Zachodzi bowiem

Twierdzenie Zbieżność CG jako metody iteracyjnej

Po $k$ iteracjach metody CG,

$||x_k-x|{|}_A\le 2{\left(\frac{\sqrt{\text{cond}(A)}-1}{\sqrt{\text{cond}(A)}+1}\right)}^k||x_0-x|{|}_A,$

gdzie $\text{cond}(A)=||A|{|}_2||A^{-1}|{|}_2$.

Metoda GMRES Generalized Minimum RESidual nie wymaga ani symetrii, ani dodatniej określoności macierzy, jest więc bardziej uniwersalna, choć też bardziej kosztowna od CG.

Prekondycjoning

Zbieżność wszystkich poznanych metod iteracyjnych zależy od własności spektralnych macierzy układu. Często w zastosowaniach pojawiające się macierze mają niekorzystne własności spektralne (np. bardzo duży wskaźnik uwarunkowania), przez co metody iteracyjne zbiegają na nich bardzo wolno.

Dlatego bardzo korzystne może być wstępne przetransformowanie układu

$Ax=b$

z macierzą o niekorzystnych własnościach, do układu

$MAx=Mb,$

,

gdzie macierz $MA$ ma znacznie korzystniejsze własności z punktu widzenia używanej metody iteracyjnej. Taką operację nazywamy prekondycjoningiem, a macierz $M$ --- macierzą prekondycjonera.

Aby całość miała sens, macierz prekondycjonera $M$ powinna:

• być łatwa w konstrukcji,
• być tania w mnożeniu przez wektor (głównym elementem każdej metody

iteracyjnej jest mnożenie macierzy przez wektor: $M\cdot (A\cdot x)$),

• macierz $MA$ powinna mieć znacznie korzystniejsze własności z punktu widzenia

używanej metody iteracyjnej.

Kilka ekstremalnych propozycji na macierz prekondycjonera to $M=I$ (łatwa w konstrukcji i tania w mnożeniu, ale nic niestety nie polepsza) oraz $M=A^{-1}$ (rewelacyjnie poprawia zbieżność metody iteracyjnej, dając zbieżność w jednej iteracji, ale bardzo droga w konstrukcji i mnożeniu). Widać więc, że należy poszukiwać czegoś pośredniego, co niskim kosztem przybliża działanie macierzy odwrotnej.

Dlatego jednym z powszechniej stosowanych rodzajów prekondycjonerów są te oparte na zastosowaniu jednego kroku klasycznej metody iteracyjnej.

Zbieżność metody CG bez prekondycjonera oraz z prekondycjonerem opartym na jednej iteracji (blokowej) metody Jacobiego.

Inne prekondycjonery stosują np. techniki tzw. niepełnego rozkładu macierzy, albo --- w specyficznych przypadkach --- tzw. metody wielosiatkowe.

Okazuje się, że zarówno CG jak i GMRES da się zaimplementować tak, by w jednej iteracji było konieczne tylko jedno mnożenie przez macierz prekondycjonera.