Zagadnienia

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=n^{d} węzłach, jej współczynnik uwarunkowania rośnie tak szybko z rozmiarem macierzy (jak n^{2}), ż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

Ax=b

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

\displaystyle MAx=Mb, (7.1)

którego macierz MA 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 MA były ciasno upakowane w kilku klastrach (najlepiej — w jednym).

Jeśli więc chcielibyśmy przekształcić macierz tak, by metoda iteracyjna dla MA 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 Mmacierzą ś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 MA 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”7W metodach operatorowych (matrix–free) nie jest nam potrzebna macierz M, tylko sposób obliczania M\cdot x, więc to ta operacja w praktyce wymaga ,,konstrukcji”. Doskonałym przykładem byłoby M=B^{{-1}}, gdzie B\approx A; wtedy M\cdot 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\cdot(A\cdot 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:

\displaystyle M_{L}AM_{R}y \displaystyle=M_{L}b, (7.2)
\displaystyle x \displaystyle=M_{R}y. (7.3)

(jeśli powyżej M_{L}=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 M_{L}) lub prawostronnej (jako M_{R}), operator M w przykładach poniżej, o ile nie zaznaczono inaczej, będzie mógł oznaczać zarówno M_{L} jak i M_{R}.

Ściskanie metodą Jacobiego

Ściskanie metodą Jacobiego możemy zastosować tylko wtedy, gdy a_{{ii}}\neq 0 dla wszystkich i. Wtedy — przez analogię do metody iteracyjnej Jacobiego — wybieramy M=\text{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\cdot y najlepiej wykonać, gdy zawczasu wyznaczymy czynniki rozkładu A=LU (co może nas drogo kosztować) i wtedy My=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(N^{3}) 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 l_{{ij}} dolnego trójkąta L oraz u_{{ij}} górnego trójkąta U, dla których a_{{ij}}\neq 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
    a_{{ik}}=a_{{ik}}/a_{{kk}}
  end
  for j=k+1 to N
    for i=k+1 to N
      a_{{ij}}=a_{{ij}}-a_{{ik}}*a_{{kj}}
    end
  end
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  a_{{ik}}\neq 0
      a_{{ik}}=a_{{ik}}/a_{{kk}}
    end
  end
  for j=k+1 to N
    for i=k+1 to N
      if a_{{ij}}\neq 0
        a_{{ij}}=a_{{ij}}-a_{{ik}}*a_{{kj}}
      end
    end
  end
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,

A_{h}x_{h}=b_{h},

indeksowana pewnym parametrem h. Rzeczywiście, w przypadku zadań takich jak równanie z macierzą dyskretnego laplasjanu (zob. przykłady 5.15.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 A_{h} są symetryczne i dodatnio określone, macierzy ściskających warto szukać w tej samej klasie.

Definicja 7.1

Rodzinę macierzy B_{h}=B_{h}^{T}>0 indeksowaną parametrem h nazywamy spektralnie równoważną rodzinie macierzy A_{h}=A_{h}^{T}>0, jeśli istnieją dwie stałe c,C niezależne od h takie, że

\displaystyle c\, x^{T}B_{h}x\leq x^{T}A_{h}x\leq C\, x^{T}B_{h}x\quad\forall x\in\mathbb{R}^{N}. (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 B_{h} jest rodziną macierzy spektralnie równoważnych macierzy A_{h}, to wartości własne \lambda(B_{h}^{{-1}}A_{h}) 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\leq\lambda(B_{h}^{{-1}}A_{h})\leq C.

Ponieważ dla dowolnej macierzy X, jej wartości własne są tożsame z wartościami macierzy B_{h}^{{1/2}}XB_{h}^{{-1/2}}, to w szczególności B_{h}^{{-1}}A_{h} ma spektrum identyczne z macierzą B_{h}^{{-1/2}}A_{h}B_{h}^{{-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 y^{T}y dostajemy

c\leq\dfrac{y^{T}B_{h}^{{-1/2}}A_{h}B_{h}^{{-1/2}}y}{y^{T}y}\leq C\quad\forall y\in\mathbb{R}^{N}.

Znaczy to, że wartości własne macierzy B_{h}^{{-1/2}}A_{h}B_{h}^{{-1/2}} leżą w przedziale [c,C].

Ćwiczenie 7.1

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

c\lambda _{{\min}}(B_{h}^{{-1}})\leq\lambda(A_{h})\leq C\lambda _{{\max}}(B_{h}^{{-1}})\quad\forall x\mathbb{R}^{N}.

(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

\dfrac{1}{C}\, x^{T}B_{h}^{{-1}}x\leq x^{T}A_{h}^{{-1}}x\leq\dfrac{1}{c}\, x^{T}B_{h}^{{-1}}x\quad\forall x\mathbb{R}^{N}.
Rozwiązanie: 

Jak wiemy z dowodu twierdzenia 7.1, \lambda(B_{h}^{{-1/2}}A_{h}B_{h}^{{-1/2}})\subset[c,C]. Ponieważ wartości własne macierzy odwrotnej są odwrotnościami wartości własnych macierzy danej, to w konsekwencji

\lambda(B_{h}^{{1/2}}A_{h}^{{-1}}B_{h}^{{1/2}})\subset[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

\dfrac{1}{C}\leq\dfrac{y^{T}B_{h}^{{1/2}}A_{h}^{{-1}}B_{h}^{{1/2}}y}{y^{T}y}\leq\dfrac{1}{c}\quad\forall y\mathbb{R}^{N}.

Podstawiając x=B_{h}^{{1/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

\displaystyle-\operatorname{div}(C(x)\nabla u(x)) \displaystyle=f(x), \displaystyle x\in\Omega=[0,1]^{d}, (7.5)
\displaystyle u(x) \displaystyle=0, \displaystyle x\in\partial\Omega, (7.6)

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

\exists\gamma>0\quad\forall x\in\Omega\qquad v^{T}C(x)v\geq\gamma v^{T}v\qquad\forall v\in\mathbb{R}^{d}.

Stąd między innymi wynika, że forma dwuliniowa określona na przestrzeni Sobolewa V=H^{1}_{0}(\Omega),

a(u,v)=\int _{\Omega}C(x)\nabla u(x)\cdot\nabla v(x)\, dx,

jest V–eliptyczna, to znaczy,

a(u,u)\geq\gamma\| u\| _{V}^{2}\qquad\forall u\in V,

gdzie \| u\| _{V}^{2}=\int _{\Omega}\nabla u(x)\cdot\nabla v(x)\, dx. Forma a(\cdot,\cdot) jest także ciągła,

a(u,v)\leq\Gamma\| u\| _{V}\| u\| _{V}\qquad\forall u,v\in V.

Stałe \gamma,\Gamma zależą oczywiście od C.

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

a(u,v)=f(v)\equiv\int _{\Omega}f(x)\, v(x)\, dx\qquad\forall v\in V

będziemy aproksymować metodą elementu skończonego w pewnej podprzestrzeni skończonego wymiaru V^{h}\subset V (por. wykład z Numerycznych równań różniczkowych). Będziemy więc szukać u^{h}\in V^{h} takiego, że

a(u^{h},v^{h})=f(v^{h})\qquad\forall v^{h}\in V^{h}.

Reprezentując8Oczywiście, N zależy od h, ale nie będziemy tego wyraźnie zaznaczać, by nie mnożyć indeksów. u^{h} w bazie \{\nu _{1}^{h},\ldots,\nu _{N}^{h}\} przestrzeni V^{h}: u^{h}=\sum _{{i=1}}^{N}U^{h}_{i}\nu _{i} mamy, że nasze zadanie dyskretne jest równoważne znalezieniu U_{h}\in\mathbb{R}^{N} spełniającego układ równań

\displaystyle A_{h}U_{h}=F_{h}, (7.7)

gdzie (A_{h})_{{ij}}=a(\nu _{i}^{h},\nu _{j}^{h}). (Macierz A_{h} 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 A_{h} jest bardzo źle uwarunkowana ze względu na h: \mbox{cond}_{2}(A_{h})=O(h^{{-2}}) i dlatego konieczne jest solidne jej ściśnięcie.

Niech B_{h} 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ą u_{h}\in V^{h} a jej reprezentacją U_{h}\in\mathbb{R}^{N} w wybranej bazie zachodzi następujący związek pomiędzy macierzą sztywności, a formą dwuliniową:

U_{h}^{T}A_{h}U_{h}=a(u_{h},u_{h}),

i w konsekwencji także

U_{h}^{T}B_{h}U_{h}=\| u_{h}\| _{V}^{2}.

Na mocy ciągłości i eliptyczności formy a w normie \|\cdot\| _{V}, macierze A_{h} i B_{h} są spektralnie równoważne.

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

Aby jednak trochę zbalansować nasz pogląd na tę sprawę warto zauważyć, że swój wkład w uwarunkowanie A_{h} ma także stosunek \Gamma/\gamma. Jeśli jest on bardzo duży (a może tak być, np. w materiałach silnie anizotropowych), wtedy uwarunkowanie B_{h}^{{-1}}A_{h} — 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, MA, w ogólnosci nie musi być dalej symetryczna — wbrew założeniom uczynionym w wyprowadzeniu metody CG.

Na szczęście, macierz MA jest symetryczna, ale w niestandardowym iloczynie skalarnym: \left\langle x,y\right\rangle=x^{T}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 M^{{1/2}}:

\displaystyle M^{{1/2}}AM^{{1/2}}\tilde{x} \displaystyle=M^{{1/2}}b, (7.8)
\displaystyle x \displaystyle=M^{{1/2}}\tilde{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 \tilde{A}=M^{{1/2}}AM^{{1/2}} jest symetryczna i dodatnio określona, można więc do niej zastosować bezpośrednio algorytm CG, wyznaczający kolejne przybliżenia \tilde{x}_{k}\in\tilde{x}_{0}+K_{k}(\tilde{r}_{0},\tilde{A}). Oznaczając \tilde{b}=M^{{1/2}}b, na k-tej iteracji wyznaczymy:

  • \tilde{r}_{k}=\tilde{b}-\tilde{A}\tilde{x}_{k},

  • \tilde{p}_{{k}}=\tilde{r}_{k}+\beta\tilde{p}_{{k-1}}

  • \rho _{k}=\|\tilde{r}_{k}\| _{2}^{2}

  • \alpha=\rho _{{k-1}}/\tilde{p}_{{k-1}}^{T}\tilde{A}\tilde{p}_{{k-1}}

  • \beta=\rho _{k}/\rho _{{k-1}}

  • \tilde{x}_{{k+1}}=\tilde{x}_{{k}}+\alpha\tilde{p}_{{k-1}}

Wprowadzając oznaczenia r_{k}=b-Ax_{k}, gdzie x_{k}=M^{{1/2}}\tilde{x}_{k}, i przyjmując p_{k}=M^{{1/2}}\tilde{p}_{{k}} możemy zapisać powyższe operacje wyłącznie w terminach wartości ,,bezfalkowych”:

  • \tilde{r}_{k}=M^{{1/2}}(b-Ax_{k})=M^{{1/2}}r_{k},

  • {p}_{{k}}={r_{k}}+\beta{p}_{{k-1}} (wynika z pomnożenia odpowiedniej równości z falkami przez M^{{1/2}})

  • \rho _{k}=\|\tilde{r}_{k}\| _{2}^{2}=r_{k}^{T}Mr_{k}

  • \alpha=\rho _{{k-1}}/\left((M^{{-1/2}}{p}_{{k-1}})^{T}M^{{1/2}}{A}M^{{1/2}}M^{{-1/2}}{p}_{{k-1}}\right)=\rho _{{k-1}}/{p}_{{k-1}}^{T}{A}{p}_{{k-1}}

  • \beta=\rho _{k}/\rho _{{k-1}}

  • {x}_{{k+1}}={x}_{{k}}+\alpha{p}_{{k-1}} (wynika z pomnożenia odpowiedniej równości z falkami przez M^{{1/2}}).

Znaczy to, że cały algorytm PCG możemy w praktyce zaimplementować tak, jak algorytm CG z tą różnicą, że do wyznaczenia \rho _{k} będziemy musieli obliczyć dodatkowy wektor Mr_{k}. Tym samym w implementacji nie pojawi się nigdzie mnożenie przez M^{{1/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=T_{N}+S,

gdzie S jest symetryczną, trójdiagonalną macierzą o losowo wybranych elementach z przedziału (0,\tau) (przy czym \tau jest małe), a T_{N} zadano wzorem (5.5).



Sprawdź, czy szybkość zbieżności PCG zależy od N i od \tau. 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 x_{k}\in x_{0}+K_{k}(Mr_{0},MA).

Rozwiązanie: 

Ponieważ \tilde{x}_{k}\in\tilde{x_{0}}+K_{k}(\tilde{r}_{0},\tilde{A}), to mnożąc stronami przez M^{{1/2}} mamy, że x_{k}\in x_{0}+M^{{1/2}}K_{k}(\tilde{r}_{0},\tilde{A}). Ponieważ

\displaystyle M^{{1/2}}\tilde{A}^{j}\tilde{r}_{0} \displaystyle=M^{{1/2}}(M^{{1/2}}AM^{{1/2}})^{j}M^{{1/2}}r_{0}
\displaystyle=M^{{1/2}}\cdot M^{{1/2}}AM^{{1/2}}\cdot M^{{1/2}}AM^{{1/2}}\cdots M^{{1/2}}AM^{{1/2}}\cdot M^{{1/2}}r_{0}
\displaystyle=(MA)^{j}Mr_{0},

to

M^{{1/2}}K_{k}(\tilde{r}_{0},\tilde{A})=K_{k}(Mr_{0},MA).

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ęciu9Jednostronne ściśnięcie dostaniemy, biorąc jedną z macierzy równą identyczności. macierzy A nieosobliwymi macierzami M_{L} i M_{R}:

M_{L}AM_{R}y=M_{L}b

wystarczy do macierzy M_{L}AM_{R} zastosować standardowy algorytm GMRES i pamiętać, by wyznaczone rozwiązanie przybliżone y_{k} przetransformować do rozwiązania oryginalnego układu:

x_{k}=M_{R}y_{k}.

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

M_{L}AM_{R}v=M_{L}(A(M_{R}v)).

Warto pamiętać, GMRES będzie obliczał residua równe r_{k}=M_{L}(b-Ax_{k}), 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
\| x_{k}-x^{*}\|\leq\epsilon _{{\text{abs}}},
Względne
\| x_{k}-x^{*}\|\leq\epsilon _{{\text{rel}}}\,\| x_{0}-x^{*}\|,
Mieszane
\| x_{k}-x^{*}\|\leq\epsilon _{{\text{abs}}}+\epsilon _{{\text{rel}}}\,\| x_{0}-x^{*}\|.

Sęk w tym, że zazwyczaj nie możemy wyznaczyć \| x_{k}-x^{*}\| z prostego powodu: nie znamy rozwiązania dokładnego x^{*}! Z drugiej strony, zawsze dysponujemy przecież inną, obliczalną wprost wartością: residuum r_{k}=b-Ax_{k}.

7.5.2. Kryterium residualne

Ponieważ r_{k}=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
\| r_{k}\|\leq\epsilon _{{\text{abs}}},
Względne
\| r_{k}\|\leq\epsilon _{{\text{rel}}}\,\| r_{0}\|,
Mieszane
\| r_{k}\|\leq\epsilon _{{\text{abs}}}+\epsilon _{{\text{rel}}}\,\| r_{0}\|.
Ćwiczenie 7.4

Podaj przykład takiej normy, w której możemy dokładnie wyznaczyć \| x_{k}-x^{*}\| na podstawie r_{k}.

Rozwiązanie: 

Pamiętajmy, że zakładamy zawsze, że A jest nieosobliwa. Ponieważ r_{k}=b-Ax_{k}=A(x^{*}-x_{k}), to

\| r_{k}\| _{2}^{2}=(x^{*}-x_{k})^{T}A^{T}A(x^{*}-x_{k})=\| x_{k}-x^{*}\| _{{A^{T}A}}^{2}.

Zatem standardowa norma euklidesowa residuum odpowiada normie energetycznej błędu liczonej w normie indukowanej przez macierz symetryczną i dodatnio określoną A^{T}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 Ax=b i niech r_{k}=b-Ax_{k} dla zadanego x_{k}. Wtedy

\dfrac{1}{\operatorname{cond}(A)}\dfrac{\| r_{k}\|}{\| b\|}\leq\dfrac{\| x_{k}-x^{*}\|}{\| x^{*}\|}\leq\operatorname{cond}(A)\dfrac{\| r_{k}\|}{\| b\|},

gdzie \operatorname{cond}(A)=\| A\|\cdot\| 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 \kappa=\lambda _{{\max}}(A)/\lambda _{{\min}}(A)=\| A\| _{2}\cdot\| A^{{-1}}\| _{2}.

7.5.3. Kryterium postępu

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

Stagnacji
\| x_{{k+1}}-x_{k}\|\leq\epsilon _{{\text{abs}}}+\epsilon _{{\text{rel}}}\,\| x_{k}\|,
Cierpliwości
k>k_{{\max}}.

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 k_{{\max}} iteracjach. Ponadto, jeżeli poprawka \| x_{{k+1}}-x_{k}\| 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ą

x_{{k+1}}=Mx_{k}+c,

przy czym \| M\|\leq\beta<1. Wykaż, że jeśli

\| x_{k}-x_{{k-1}}\|\leq\epsilon\dfrac{1-\beta}{\beta},

to

\| x_{k}-x^{*}\|\leq\epsilon.
Rozwiązanie: 

Ponieważ w metodzie stacjonarnej x^{*}=Mx^{*}+c, to

\displaystyle\| x_{k}-x^{*}\| \displaystyle=\| Mx_{{k-1}}+c-x^{*}\|=\| M(x_{{k-1}}-x^{*})\|\leq\| M\|\| x_{{k-1}}-x^{*}\|\leq\beta(\| x_{{k-1}}-x_{k}\|+\| x_{{k}}-x^{*}\|)
\displaystyle\leq\epsilon(1-\beta)+\beta\| x_{k}-x^{*}\|.

Zatem \| x_{k}-x^{*}\|(1-\beta)\leq\epsilon(1-\beta), skąd teza.

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

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.