Zagadnienia

7.1 Proste operatory ściskające
7.2 Macierze spektralnie równoważne
7.3 Metoda PCG
7.4 Ściskanie dla GMRES
7.5 Kryteria stopu metod iteracyjnych

7. Ściskanie macierzy (preconditioning)

Zbieżność wszystkich dotychczas przez nas poznanych metod iteracyjnych projekcji i Kryłowa zależy od własności spektralnych macierzy rozwiązywanego układu. Pojawiające się w zastosowaniach macierze często mają niestety niekorzystne własności spektralne (np. bardzo duży wskaźnik uwarunkowania), przez co metody iteracyjne zbiegają na nich bardzo wolno.

Przykładowo, dla wielokrotnie tu przywoływanej macierzy d-wymiarowego laplasjanu dyskretyzowanego na siatce o N=nd węzłach, jej współczynnik uwarunkowania rośnie tak szybko z rozmiarem macierzy (jak n2), że czyni bezsensownym iteracyjne rozwiązywanie takich układów gdy n jest bardzo duże. . Zgodnie z twierdzeniem o zbieżności CG, musielibyśmy wykonać rzędu O⁢n iteracji, by dostatecznie zredukować błąd. Jasne jest, że jeśli n=10000, to koszt iteracji staje się zaporowy, a w międzyczasie dodatkowo mogą odezwać się problemy związane z realizacją iteracji w arytmetyce skończonej precyzji.

Dlatego bardzo korzystna może być wstępna transformacja układu

A⁢x=b

z macierzą o niekorzystnych własnościach, do równoważnego układu

M⁢A⁢x=M⁢b,

(7.1)

którego macierz M⁢A ma znacznie korzystniejsze własności z punktu widzenia używanej metody iteracyjnej.

Okazuje się, że w m.in. przypadku macierzy laplasjanu i podobnych, można faktycznie wskazać taką macierz M, że iteracja nie jest wiele droższa, ale za to znacznie szybciej zbieżna. Z drugiej strony, w innych ważnych z punktu widzenia zastosowań problemach, takiej dobrej macierzy M wciąż nie umiemy wskazać.

W przypadku macierzy symetrycznych widzieliśmy, że kluczowe znaczenie dla zbieżności metody miało rozmieszczenie wartości własnych na osi liczbowej: jeśli wartości własne były bardzo rozrzucone po prostej, to uwarunkowanie było bardzo duże i w konsekwencji zbieżność powolna. Aby zbieżność była szybsza, kluczowe jest, by wartości własne M⁢A były ciasno upakowane w kilku klastrach (najlepiej — w jednym).

Jeśli więc chcielibyśmy przekształcić macierz tak, by metoda iteracyjna dla M⁢A zbiegała szybko, musimy w jakiś sposób ,,ścisnąć” spektrum macierzy A. Taką operację nazywamy więc ściskaniem (ang. preconditioning — bo zmiana uwarunkowania), a macierz M — macierzą ściskającą.

Aby całość miała sens, macierz ściskająca M powinna być tak dobrana, by jednocześnie zachodziło kilka nieco przeciwstawnych warunków:

macierz M⁢A powinna mieć znacznie korzystniejsze własności z punktu widzenia szybkości zbieżności (i ostatecznego kosztu) używanej metody iteracyjnej,
macierz M powinna być łatwa w ,,konstrukcji”⁷W metodach operatorowych (matrix–free) nie jest nam potrzebna macierz M, tylko sposób obliczania M⋅x, więc to ta operacja w praktyce wymaga ,,konstrukcji”. Doskonałym przykładem byłoby M=B-1, gdzie B≈A; wtedy M⋅x to obliczenie rozwiązania układu równań z macierzą B, więc ,,konstrukcją” M byłoby znalezienie rozkładu LU macierzy B.,
macierz M powinna być tania w mnożeniu przez wektor (głównym elementem każdej metody iteracyjnej jest mnożenie macierzy przez wektor: tutaj mamy M⋅A⋅x).

Kilka ekstremalnych propozycji na macierz ściskającą to M=I (łatwa w konstrukcji i tania w mnożeniu, ale niestety nic nie polepsza…) oraz M=A-1 (rewelacyjnie poprawia szybkość zbieżności metody iteracyjnej, dając zbieżność w jednej iteracji, ale jest droga w mnożeniu i bardzo droga w konstrukcji). Ponieważ obie skrajności są — ze zrozumiałych względów — nie do zaakceptowania, widać więc, że należy poszukiwać czegoś pośredniego: czegoś, co niskim kosztem przybliża (w jakimś sensie) działanie macierzy odwrotnej.

Uogólnieniem ściskania lewostronego (7.1) jest ściskanie obustronne:

	ML⁢A⁢MR⁢y	=ML⁢b,		(7.2)
	x	=MR⁢y.		(7.3)

(jeśli powyżej ML=I, to ściskanie jest oczywiście tylko prawostronne).

7.1. Proste operatory ściskające

Jednym z powszechniej stosowanych (aczkolwiek nie zawsze dostatecznie skutecznych) sposobów ściskania są te oparte na zastosowaniu jednego kroku klasycznej metody iteracyjnej (lub jej wersji blokowej). Ponieważ mogą one być stosowane zarówno w wersji lewostronnej (jako ML) lub prawostronnej (jako MR), operator M w przykładach poniżej, o ile nie zaznaczono inaczej, będzie mógł oznaczać zarówno ML jak i MR.

Ściskanie metodą Jacobiego

Ściskanie metodą Jacobiego możemy zastosować tylko wtedy, gdy ai⁢i≠0 dla wszystkich i. Wtedy — przez analogię do metody iteracyjnej Jacobiego — wybieramy M=diag⁢A-1. Oprócz zasadniczej wady, jaką jest najczęściej kiepska skuteczność w redukcji uwarunkowania, ściskanie metodą Jacobiego ma poza tym wiele zalet:

generowanie M jest praktycznie darmowe,
mnożenie przez M jest bardzo tanie (O⁢N flopów), a przy tym pięknie się zrównolegla.

Metody niepełnego rozkładu

Inne sposoby ściśnięcia macierzy obejmują np. techniki tzw. niepełnego rozkładu macierzy. Ich idea opiera się na obserwacji, że ,,idealną” macierzą ściskającą byłaby M=A-1. Przypomnijmy, że w takim przypadku mnożenie M⋅y najlepiej wykonać, gdy zawczasu wyznaczymy czynniki rozkładu A=L⁢U (co może nas drogo kosztować) i wtedy M⁢y=U-1⁢L-1⁢y (tu wykorzystamy fakt, że układy równań z macierzami L i U są ,,łatwiejsze” do rozwiązania, bo macierze są trójkątne). Oczywiście, taka ,,idealna” macierz ściskająca ma wadę, która w większości przypadków będzie ją całkowicie dyskwalifikować: mimo, że A jest rozrzedzona, macierze L i U mogą — i zazwyczaj faktycznie są — gęste. Oznacza to, że samo wyznaczenie rozkładu LU macierzy A będzie nas kosztować zabójcze O⁢N3 flopów.

Inżynierski pomysł na usunięcie tej wady jest taki, by w algorytmie wyznaczania czynników rozkładu LU wyznaczać jedynie te elementy li⁢j dolnego trójkąta L oraz ui⁢j górnego trójkąta U, dla których ai⁢j≠0. Zatem struktura niezerowych elementów macierzy L,U będzie powielać strukturę niezerowych elementów A — a więc dostaniemy w efekcie macierze rzadkie!

Jeśli więc ogólnie algorytm rozkładu LU (w miejscu, bez wyboru elementu głównego) jest postaci:

for k=1 to N-1

for i=k+1 to N

ai⁢k=ai⁢k/ak⁢k

end

for j=k+1 to N

for i=k+1 to N

ai⁢j=ai⁢j-ai⁢k*ak⁢j

end

to niepełny rozkład LU miałby (formalnie) postać:

Niepełny rozkład LU, wersja bazowa

for k=1 to N-1

for i=k+1 to N

if ai⁢k≠0

ai⁢k=ai⁢k/ak⁢k

end

for j=k+1 to N

for i=k+1 to N

if ai⁢j≠0

ai⁢j=ai⁢j-ai⁢k*ak⁢j

end

Opisana powyżej technika nosi nazwę niepełnego rozkładu LU (ang. incomplete LU factorization) o zerowym wypełnieniu (ang. zero fill-in), co oznaczamy ILU(0). Warto zwrócić uwagę, że może zdarzyć się, że czynnik U uzyskany niepełnym rozkładem macierzy nieosobliwej sam będzie macierzą osobliwą. Dla pewnych klas macierzy — np. dla tzw. M-macierzy [13] — wiadomo, że niepełny rozkład nie prowadzi do takich degeneracji.

Często rozkład ILU(0) produkuje czynniki LU, dla których U-1⁢L-1 jest tak bardzo odległe od A-1, że nie dają one dostatecznie silnego ściśnięcia macierzy A. Wtedy można zastosować któryś z wariantów niepełnego rozkładu

rozkład z progiem, ILUT (ang. treshold ILU), w którym zerujemy te elementy czynników rozkładu, które spełniają pewien warunek progowy — oględnie mówiąc, warunek ten dotyczy względnej wielkości danego elementu (,,pomijamy małe”);
rozkład według wzorca: w którym wprost wskazujemy, jakie elementy L i U mogą przyjąć niezerowe wartości;
rozkład według wypełnienia, ILU(p): w którym dopuszczamy pewne dodatkowe wypełnienie czynników rozkładu.

Szczegółowe omówienie różnych technik niepełnego rozkładu macierzy, razem z ich analizą w niektórych praktycznie użytecznych przypadkach, znajdziemy w monografii Saada [13].

Metody dopasowane do problemu

Najskuteczniejsze operatory ściskające dostajemy wykorzystując do maksimum całą specyfikę konkretnego rozwiązywanego zadania. Na przykład, dla macierzy pochodzących z dyskretyzacji niektórych równań różniczkowych cząstkowych opracowano kilka klas znakomitych sposobów ściskania takich, które gwarantują bardzo szybką zbieżność metod Kryłowa. Wśród nich znajdują się metody dekompozycji obszaru (por. rozdział 5.4, [14], [16], []) oraz metody wielosiatkowe. O tych ostatnich piszemy nieco więcej w rozdziale 8.3.

7.2. Macierze spektralnie równoważne

W wielu zastosowaniach, na przykład w wielokrotnie tu przywoływanych zadaniach dyskretyzacji eliptycznych równań różniczkowych cząstkowych, przedmiotem naszego zainteresowania jest nie tyle jeden układ równań, ale cała ich rodzina,

Ah⁢xh=bh,

indeksowana pewnym parametrem h. Rzeczywiście, w przypadku zadań takich jak równanie z macierzą dyskretnego laplasjanu (zob. przykłady 5.1 i 5.2), naszym parametrem h byłby bok oczka siatki dyskretyzacji.

Powstaje więc naturalna potrzeba wskazania klasy macierzy ściskających, których użycie pozwalałoby uzyskać szybkość zbieżności metody iteracyjnej niezależną od h. W przypadku, gdy macierze Ah są symetryczne i dodatnio określone, macierzy ściskających warto szukać w tej samej klasie.

Definicja 7.1

Rodzinę macierzy Bh=BhT>0 indeksowaną parametrem h nazywamy spektralnie równoważną rodzinie macierzy Ah=AhT>0, jeśli istnieją dwie stałe c,C niezależne od h takie, że

c⁢⁢xT⁢Bh⁢x≤xT⁢Ah⁢x≤C⁢⁢xT⁢Bh⁢x⁢⁢∀x∈RN.

(7.4)

Ważną własnością macierzy spektralnie równoważnych jest to, że na ich podstawie można skonstruować rodzinę macierzy ściskających, prowadzących do zadania o własnościach spektralnych niezależnych od h, por. [6].

Twierdzenie 7.1

Jeśli Bh jest rodziną macierzy spektralnie równoważnych macierzy Ah, to wartości własne λ⁢Bh-1⁢Ah są ograniczone niezależnie od h.

Dowód

Pokażemy nieco więcej: jeśli zachodzi warunek spektralnej równoważności (7.4), to dla dowolnego h

c≤λ⁢Bh-1⁢Ah≤C.

Ponieważ dla dowolnej macierzy X, jej wartości własne są tożsame z wartościami macierzy Bh1/2⁢X⁢Bh-1/2, to w szczególności Bh-1⁢Ah ma spektrum identyczne z macierzą Bh-1/2⁢Ah⁢Bh-1/2. Ta ostatnia jest symetryczna i dodatnio określona, więc wszystkie jej wartości własne są rzeczywiste i dodatnie.

Podstawiając x=B-1/2⁢y w (7.4) i dzieląc stronami przez yT⁢y dostajemy

c≤yT⁢Bh-1/2⁢Ah⁢Bh-1/2⁢yyT⁢y≤C⁢⁢∀y∈RN.

Znaczy to, że wartości własne macierzy Bh-1/2⁢Ah⁢Bh-1/2 leżą w przedziale c,C.

∎

Ćwiczenie 7.1

Wykaż, że jeśli zachodzi (7.4), to

c⁢λmin⁢Bh-1≤λ⁢Ah≤C⁢λmax⁢Bh-1⁢⁢∀x⁢RN.

(Por. [6, 3.10].)

Wskazówka:

Wynika z twierdzenia Couranta–Fishera.

Ćwiczenie 7.2

Wykaż, że jeśli zachodzi (7.4), to jednocześnie zachodzi

1C⁢⁢xT⁢Bh-1⁢x≤xT⁢Ah-1⁢x≤1c⁢⁢xT⁢Bh-1⁢x⁢⁢∀x⁢RN.

Rozwiązanie:

Jak wiemy z dowodu twierdzenia 7.1, λ⁢Bh-1/2⁢Ah⁢Bh-1/2⊂c,C. Ponieważ wartości własne macierzy odwrotnej są odwrotnościami wartości własnych macierzy danej, to w konsekwencji

λ⁢Bh1/2⁢Ah-1⁢Bh1/2⊂1/C,1/c.

Korzystając z faktu, że iloraz Rayleigh jest ograniczony z góry i z dołu przez ekstremalne wartości własne, dostajemy

1C≤yT⁢Bh1/2⁢Ah-1⁢Bh1/2⁢yyT⁢y≤1c⁢⁢∀y⁢RN.

Podstawiając x=Bh1/2⁢y otrzymujemy tezę.

Przykład 7.1 (Wykorzystanie szybkiego solvera do ściśnięcia macierzy dyskretyzacji równania różniczkowego o zmiennych współczynnikach)

Rozważmy eliptyczne równanie różniczkowe

	-div⁡C⁢x⁢∇⁡u⁢x	=f⁢x,		x∈Ω=0,1d,		(7.5)
	u⁢x	=0,		x∈∂⁡Ω,		(7.6)

gdzie C⁢x jest rzeczywistą i symetryczną macierzą rozmiaru d×d, o współczynnikach będących ciągłymi, ograniczonymi funkcjami na Ω, spełniającą warunek (gwarantujący eliptyczność w/w równania różniczkowego)

∃γ>0⁢⁢∀x∈Ω⁢⁢vT⁢C⁢x⁢v≥γ⁢vT⁢v⁢⁢∀v∈Rd.

Stąd między innymi wynika, że forma dwuliniowa określona na przestrzeni Sobolewa V=H01⁢Ω,

a⁢u,v=∫ΩC⁢x⁢∇⁡u⁢x⋅∇⁡v⁢x⁢⁢d⁢x,

jest V–eliptyczna, to znaczy,

a⁢u,u≥γ⁢uV2⁢⁢∀u∈V,

gdzie uV2=∫Ω∇⁡u⁢x⋅∇⁡v⁢x⁢⁢d⁢x. Forma a(⋅,⋅) jest także ciągła,

a⁢u,v≤Γ⁢uV⁢uV⁢⁢∀u,v∈V.

Stałe γ,Γ zależą oczywiście od C.

Nasze zadanie w postaci wariacyjnej: znaleźć u∈V takie, że

a⁢u,v=f⁢v≡∫Ωf⁢x⁢⁢v⁢x⁢⁢d⁢x⁢⁢∀v∈V

będziemy aproksymować metodą elementu skończonego w pewnej podprzestrzeni skończonego wymiaru Vh⊂V (por. wykład z Numerycznych równań różniczkowych). Będziemy więc szukać uh∈Vh takiego, że

a⁢uh,vh=f⁢vh⁢⁢∀vh∈Vh.

Reprezentując⁸Oczywiście, N zależy od h, ale nie będziemy tego wyraźnie zaznaczać, by nie mnożyć indeksów. uh w bazie ν1h,…,νNh przestrzeni Vh: uh=∑i=1NUih⁢νi mamy, że nasze zadanie dyskretne jest równoważne znalezieniu Uh∈RN spełniającego układ równań

Ah⁢Uh=Fh,

(7.7)

gdzie Ahi⁢j=a⁢νih,νjh. (Macierz Ah nazywana jest macierzą sztywności.) Gdy d=3, metody bezpośrednie mogą być mało efektywne wskutek problemów z wypełnieniem czynników rozkładu, należy więc układ (7.7) rozwiązywać iteracyjnie. Niestety, w przypadku, gdy bazę wybierzemy jako standardową bazę elementu skończonego, macierz Ah jest bardzo źle uwarunkowana ze względu na h: cond2⁢Ah=O⁢h-2 i dlatego konieczne jest solidne jej ściśnięcie.

Niech Bh będzie macierzą sztywności uzyskaną w przypadku, gdy C⁢x=I, czyli — ,,zwykłą” macierzą dyskretnego laplasjanu.

Zauważmy, że przy zachowaniu wyżej wspomnianej odpowiedniości pomiędzy funkcją uh∈Vh a jej reprezentacją Uh∈RN w wybranej bazie zachodzi następujący związek pomiędzy macierzą sztywności, a formą dwuliniową:

UhT⁢Ah⁢Uh=a⁢uh,uh,

i w konsekwencji także

UhT⁢Bh⁢Uh=uhV2.

Na mocy ciągłości i eliptyczności formy a w normie ∥⋅∥V, macierze Ah i Bh są spektralnie równoważne.

A więc, znając szybką metodę rozwiązywania zadania z macierzą Bh — lub, ogólniej, dobrą macierz ściskającą dla Bh — dysponujemy jednocześnie dobrym (czytaj: niezależnym od h) sposobem ściśnięcia macierzy Ah.

Aby jednak trochę zbalansować nasz pogląd na tę sprawę warto zauważyć, że swój wkład w uwarunkowanie Ah ma także stosunek Γ/γ. Jeśli jest on bardzo duży (a może tak być, np. w materiałach silnie anizotropowych), wtedy uwarunkowanie Bh-1⁢Ah — choć już nie będzie zależało od h — jednak wciąż będzie wielkie: nad tym, jak najlepiej ściskać takie macierze są obecnie prowadzone badania.

7.3. Metoda PCG

Naturalnym kandydatem na macierz ściskającą w metodzie CG (która jest określona dla macierzy A symetrycznej i dodatnio określonej) jest inna macierz symetryczna i dodatnio określona, M. Jak widzieliśmy w twierdzeniu 7.1 i ćwiczeniu 7.2, dobrym kandydatem na M mogłaby być macierz spektralnie równoważna macierzy A-1. Wydawać by się jednak mogło, że jednostronne (na przykład, lewostronne) ściskanie taką macierzą nie powiedzie się, bo otrzymana przez nas macierz, M⁢A, w ogólnosci nie musi być dalej symetryczna — wbrew założeniom uczynionym w wyprowadzeniu metody CG.

Na szczęście, macierz M⁢A jest symetryczna, ale w niestandardowym iloczynie skalarnym: x,y=xT⁢M-1⁢y. Aby uniknąć przeformułowania wyników uzyskanych we wcześniejszych rozdziałach, użyjemy pewnego triku: teoretycznie użyjemy zadania symetrycznie ściśniętego obustronnie macierzami M1/2:

	M1/2⁢A⁢M1/2⁢x~	=M1/2⁢b,		(7.8)
	x	=M1/2⁢x~,		(7.9)

ale implementację doprowadzimy do takiego stanu, że będzie w niej występować tylko mnożenie wektora przez macierz M — i to tylko jeden raz w każdej iteracji!

Rzeczywiście, ściśnięta macierz A~=M1/2⁢A⁢M1/2 jest symetryczna i dodatnio określona, można więc do niej zastosować bezpośrednio algorytm CG, wyznaczający kolejne przybliżenia x~k∈x~0+Kk⁢r~0,A~. Oznaczając b~=M1/2⁢b, na k-tej iteracji wyznaczymy:

r~k=b~-A~⁢x~k,
p ~ k = r ~ k + β ⁢ p ~ k - 1
ρ k = r ~ k 2 2
α = ρ k - 1 / p ~ k - 1 T ⁢ A ~ ⁢ p ~ k - 1
β = ρ k / ρ k - 1
x ~ k + 1 = x ~ k + α ⁢ p ~ k - 1

Wprowadzając oznaczenia rk=b-A⁢xk, gdzie xk=M1/2⁢x~k, i przyjmując pk=M1/2⁢p~k możemy zapisać powyższe operacje wyłącznie w terminach wartości ,,bezfalkowych”:

r~k=M1/2⁢b-A⁢xk=M1/2⁢rk,
pk=rk+β⁢pk-1 (wynika z pomnożenia odpowiedniej równości z falkami przez M1/2)
ρ k = r ~ k 2 2 = r k T ⁢ M ⁢ r k
α = ρ k - 1 / M - 1 / 2 ⁢ p k - 1 T ⁢ M 1 / 2 ⁢ A ⁢ M 1 / 2 ⁢ M - 1 / 2 ⁢ p k - 1 = ρ k - 1 / p k - 1 T ⁢ A ⁢ p k - 1
β = ρ k / ρ k - 1
xk+1=xk+α⁢pk-1 (wynika z pomnożenia odpowiedniej równości z falkami przez M1/2).

Znaczy to, że cały algorytm PCG możemy w praktyce zaimplementować tak, jak algorytm CG z tą różnicą, że do wyznaczenia ρk będziemy musieli obliczyć dodatkowy wektor M⁢rk. Tym samym w implementacji nie pojawi się nigdzie mnożenie przez M1/2 — w każdej iteracji będzie potrzebne tylko jedno mnożenie przez M!

Przykład 7.2

Działanie PCG z dobrą macierzą ściskającą prześledzimy na przykładzie macierzy jednowymiarowego laplasjanu, o losowo zaburzonych elementach:

A=TN+S,

gdzie S jest symetryczną, trójdiagonalną macierzą o losowo wybranych elementach z przedziału 0,τ (przy czym τ jest małe), a TN zadano wzorem (5.5).

%
% rozwiazania zadania z macierza rozrzedzona
% 
disp('Matematyka obliczeniowa II');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = jacobi(A, b, tol=1e-10, maxit=1000, x0)
%
% metoda Jacobiego na macierzy A
%
if (nargin < 5)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy: 
% A = M - Z
M = diag(A); Z = tril(A,-1) + triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = (b - Z*x)./M;
	iter = iter+1;
	resid(iter) = norm(b-A*x);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = sor(A, b, omega=1, tol=1e-10, maxit=1000, x0)
%
% metoda SOR (gdy omega=1, to Gaussa-Seidela) na macierzy A
%
if (nargin < 6)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0;
resid(iter) = norm(b - A*x); 
% rozklad macierzy: 
% A = M - Z;
M = tril(A,-1) + spdiag(diag(A)/omega);
Z = spdiag((1/omega - 1)*diag(A)) - triu(A,1);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	x = M \ (Z*x + b);
	iter = iter+1;
	resid(iter) = norm(b-A*x);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, info, iter, resid] = steepestdescent(A, b, tol=1e-8, maxit=1000, x0)
%
% metoda najszybszego spadku na macierzy A
% wymaga A = A' > 0
%
if (nargin < 5)
	x0 = zeros(size(b));
end
resid = zeros(maxit,1);
iter = 1; 
x = x0; r = b - A*x;
resid(iter) = norm(r);

while( (resid(iter) > tol*resid(1)) && (info = (iter < maxit)))
	w = A*r;
	a = (r'*r)/(r'*w);
	x = x + a*r;
	r = r - a*w;
	iter = iter+1;
	resid(iter) = norm(r);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nrange = [10, 100, 200, 400]; % zestaw wartosci testowych parametru p
for k = 1:length(Nrange)
	N = Nrange(k); % rozmiar macierzy
	b = rand(N,1);
	% losowe zaburzenie macierza trojdiagonalna
	g = 0.1; S = spdiags(g*rand(size(e,3)), [-1,0,1], N,N);
	% losowe zaburzenie macierza nieregularna
	% g = 0.01; S = g*sprand(N,N,2/N); 
	S = 0.5*(S+S');
	e = ones(N,1); M = spdiags([-e, 2*e, -e], [-1,0,1], N,N);
	A = M + S;
	spy(A); pause(0.5);
	
	omega = 2/(1+sin(pi/(N+1))); % optymalny parametr dla SOR, gd macierz jest 1 wym laplasjanu
	% omega = 1.3;
	tic; [x, info, its, resid] = jacobi(A, b); 
	iter(k) = its; time(k) = toc;
	tic; [x, info, its, residg] = sor(A, b, omega); 
	iterg(k) = its; timeg(k) = toc; omegag(k) = omega;
	tic; [x, info, relres, its, resids] = pcg(A, b, 1e-8, 1000); 
	iters(k) = its; times(k) = toc; 
	tic; [x, info, relres, its, residc] = pcg(A, b, 1e-8, 1000, M); 
	iterc(k) = its; timec(k) = toc; 
	tmp = zeros(size(resid)); tmp(1:length(residc)) = residc; residc=tmp; 
	mi = min([iter(k),iterg(k),iters(k),iterc(k)]);
	semilogy(1:mi,[resid(1:mi), residg(1:mi), resids(1:mi), residc(1:mi)],[['b-;N=',num2str(N),' Jacobi;']; ['r-;N=',num2str(N),' (SOR(',num2str(omega),');']; ['g-;N=',num2str(N),' (CG);']; ['k-;N=',num2str(N),' (PCG);']]); pause(1.5);
end

% dla porownania - solver bezposredni

tic; x = A\b; timed = toc;

Sprawdź, czy szybkość zbieżności PCG zależy od N i od τ. Zobacz, jak zmienią się rezultaty, gdy zaburzenie nie będzie zachowywać symetrii lub dodatniej określoności A. A co stanie się, gdy S nie będzie trójdiagonalna?

Ćwiczenie 7.3

Wykaż, że w metodzie PCG xk∈x0+Kk⁢M⁢r0,M⁢A.

Rozwiązanie:

Ponieważ x~k∈x0~+Kk⁢r~0,A~, to mnożąc stronami przez M1/2 mamy, że xk∈x0+M1/2⁢Kk⁢r~0,A~. Ponieważ

	M1/2⁢A~j⁢r~0	=M1/2⁢M1/2⁢A⁢M1/2j⁢M1/2⁢r0
		=M1/2⋅M1/2⁢A⁢M1/2⋅M1/2⁢A⁢M1/2⁢⋯⁢M1/2⁢A⁢M1/2⋅M1/2⁢r0
		=M⁢Aj⁢M⁢r0,

M1/2⁢Kk⁢r~0,A~=Kk⁢M⁢r0,M⁢A.

7.4. Ściskanie dla GMRES

Ponieważ metoda GMRES nie wymaga żadnych specjalnych własności macierzy (z wyjątkiem nieosobliwości), po obustronnym ściśnięciu⁹Jednostronne ściśnięcie dostaniemy, biorąc jedną z macierzy równą identyczności. macierzy A nieosobliwymi macierzami ML i MR:

ML⁢A⁢MR⁢y=ML⁢b

wystarczy do macierzy ML⁢A⁢MR zastosować standardowy algorytm GMRES i pamiętać, by wyznaczone rozwiązanie przybliżone yk przetransformować do rozwiązania oryginalnego układu:

xk=MR⁢yk.

Wszelkie mnożenia przez macierz ściśniętą będziemy oczywiście realizować zgodnie z nawiasami:

ML⁢A⁢MR⁢v=ML⁢A⁢MR⁢v.

Warto pamiętać, GMRES będzie obliczał residua równe rk=ML⁢b-A⁢xk, zatem tylko gdy stosujemy prawostronne ściskanie, residua wyznaczone algorytmem GMRES będą identyczne z residuami układu wyjściowego (co, generalnie, samo w sobie niekoniecznie musi być wartością na której powinno nam zależeć).

7.5. Kryteria stopu metod iteracyjnych

Do tej pory zajmowaliśmy się samymi metodami przybliżonego rozwiązywania równań, zwracając uwagę na koszt jednej iteracji, sposób implementacji czy też szybkość zbieżności. Jednak każdą iterację musimy kiedyś w końcu przerwać; musimy więc zadać sobie pytanie: kiedy i na jakiej podstawie powinniśmy zakończyć działanie metody iteracyjnej?

7.5.1. Kryterium błędu

Idealnym kryterium stopu metody iteracyjnej byłaby redukcja błędu poniżej zadanego poziomu, mierzonego w zadanej przez użytkownika normie. Kryterium błędu zatrzymania metody iteracyjnej może więc mieć jedną z trzech postaci:

Bezwzględne

xk-x*≤ϵabs,

Względne

xk-x*≤ϵrel⁢⁢x0-x*,

Mieszane

xk-x*≤ϵabs+ϵrel⁢⁢x0-x*.

Sęk w tym, że zazwyczaj nie możemy wyznaczyć xk-x* z prostego powodu: nie znamy rozwiązania dokładnego x*! Z drugiej strony, zawsze dysponujemy przecież inną, obliczalną wprost wartością: residuum rk=b-A⁢xk.

7.5.2. Kryterium residualne

Ponieważ rk=0 odpowiada znalezieniu rozwiązania dokładnego, to kryterium stopu możemy oprzeć na redukcji residuum poniżej zadanego poziomu. Kryterium residualne zatrzymania metody iteracyjnej może więc mieć jedną z trzech postaci:

Bezwzględne

rk≤ϵabs,

Względne

rk≤ϵrel⁢⁢r0,

Mieszane

rk≤ϵabs+ϵrel⁢⁢r0.

Ćwiczenie 7.4

Podaj przykład takiej normy, w której możemy dokładnie wyznaczyć xk-x* na podstawie rk.

Rozwiązanie:

Pamiętajmy, że zakładamy zawsze, że A jest nieosobliwa. Ponieważ rk=b-A⁢xk=A⁢x*-xk, to

rk22=x*-xkT⁢AT⁢A⁢x*-xk=xk-x*AT⁢A2.

Zatem standardowa norma euklidesowa residuum odpowiada normie energetycznej błędu liczonej w normie indukowanej przez macierz symetryczną i dodatnio określoną AT⁢A.

Jednak zależność między normą residuum a tą samą normą błędu może nie być bardzo wyrazista: po raz kolejny może tu dać znać o sobie złe uwarunkowanie macierzy.

Stwierdzenie 7.1 (o oszacowaniu błędu przez residuum)

Niech x* będzie rozwiązaniem układu A⁢x=b i niech rk=b-A⁢xk dla zadanego xk. Wtedy

1cond⁡A⁢rkb≤xk-x*x*≤cond⁡A⁡rkb,

gdzie cond⁡A=A⋅A-1.

Powyższe stwierdzenie należy interpretować tak, że jedynie w przypadku, gdy macierz A jest dobrze uwarunkowana w danej normie, mała norma residuum jest dobrym indykatorem małej (tej samej) normy błędu. Warto wiedzieć, że w metodzie CG można niskim kosztem (O⁢k, gdzie k jest liczbą iteracji) wyznaczyć (coraz lepsze, wraz z postępem iteracji) dolne oszacowanie na spektralny współczynnik uwarunkowania κ=λmax⁢A/λmin⁢A=A2⋅A-12.

7.5.3. Kryterium postępu

W praktyce obliczeniowej stosuje się dwa dodatkowe kryteria stopu metody iteracyjnej

Stagnacji

xk+1-xk≤ϵabs+ϵrel⁢⁢xk,

Cierpliwości

k>kmax.

Rzeczywiście, musimy przecież określić moment, kiedy należałoby zakończyć obliczenia, jeśli do tej pory nie przyniosły sukcesu: zrobimy to więc po kmax iteracjach. Ponadto, jeżeli poprawka xk+1-xk jest niewielka, to metoda być może na dobre ugrzęzła w jakimś przybliżeniu i wtedy, jeśli nadal jesteśmy daleko od rozwiązania, należałoby spróbować czegoś innego. To kryterium należy stosować z dużym wyczuciem: być może stagnacja metody była tylko chwilowa i po kilku następnych iteracjach zbieżność przyspieszy?

Ćwiczenie 7.5

Rozważmy stacjonarną metodę iteracyjną

xk+1=M⁢xk+c,

przy czym M≤β<1. Wykaż, że jeśli

xk-xk-1≤ϵ⁢1-ββ,

xk-x*≤ϵ.

Rozwiązanie:

Ponieważ w metodzie stacjonarnej x*=M⁢x*+c, to

	xk-x*	=∥Mxk-1+c-x∥=∥M(xk-1-x)∥≤∥M∥∥xk-1-x∥≤β(∥xk-1-xk∥+∥xk-x∥)
		≤ϵ⁢1-β+β⁢xk-x*.

Zatem xk-x*⁢1-β≤ϵ⁢1-β, skąd teza.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.