Zagadnienia

8. Przegląd innych metod rozwiązywania wielkich układów równań liniowych

Na szczęście, nie ma jednej najlepszej metody rozwiązywania wielkich układów równań liniowych; czyni to oczywiście życie znacznie ciekawszym. Doprowadziło to do skonstruowania dużej liczby konkurencyjnych metod: zarówno bezpośrednich, jak i iteracyjnych. Celem nieniejszego rozdziału jest pobieżne zaznajomienie Czytelnika z niewielką ich reprezentacją. Dobór tej właściwej dla konkretnego zadania nie musi być oczywisty! Szerszy przegląd takich metod znajdziemy m.in. w monografii [13], [].

8.1. Inne metody Kryłowa

Na początek, na podstawie [1] przedstawimy w pewien systematyczny sposób uogólnienie metody CG na szerszą klasę problemów. Następnie zrobimy krótki przegląd wybranych metod Kryłowa, nie opartych na zasadzie minimalizacji.

8.1.1. Uogólnienia metody CG

Wielką zaletą metody PCG jest to, że nie wymaga ona pamiętania całej bazy przestrzeni Kryłowa KkMr0,MA=spanMr0,MAMr0,,MAk-1Mr0, a iteracja ma prostą formułę xk+1=xk+dk, gdzie dkKkMr0,MA jest tak dobrany by minimalizować pewną normę błędu. Można więc pójść krok dalej i określić całą klasę metod o podobnej strukturze.

Przyjmijmy więc, że

  • macierz A jest nieosobliwa (na tym poziomie ogólności nie zakładamy symetrii ani dodatniej określoności)

oraz są dane dwie dodatkowe macierze:

  • nieosobliwa macierz lewostronnie ściskająca M,

  • symetryczna, dodatnio określona macierz B, która będzie definiować normę, w której będziemy minimalizować błąd:

    xB=xTBx1/2.

Będziemy zakładali, że macierz ściśnięta MA jest normalna względem iloczynu skalarnego indukowanego przez macierz B, co możemy zapisać w formie warunku

GGT=GTG,

gdzie G=B1/2MAB-1/2.

W uogólnionej metodzie CG, którą będziemy oznaczali GCG(B, M, A), kolejna iteracja xk+1x0+KkMr0,MA będzie określona tak, by

xk+1-x*Bx-x*Bxx0+KkMr0,MA(8.1)
Quiz 8.1
  1. Czy metoda CG to nic innego jak GCG(I, I, A)?


    Dobrze Oczywiście, CG to nic innego jak GCG(A, I, A)

    Źle Sprawdź, w jakiej normie CG minimalizuje błąd

  2. Czy metoda PCG z macierzą ściskającą M to nic innego jak GCG(A, M, A)?


    Dobrze Oczywiście

    Źle Zob. [1]

Twierdzenie 8.1

Przy powyższych założeniach, metoda GCG(B, M, A) jest dobrze określona i w dokładnej arytmetyce osiąga rozwiązanie dokładne po co najwyżej N krokach.

Dowód

Kładąc r~0=Mr0 oraz A~=MA, możemy zastosować twierdzenie 6.1 do przestrzeni Kryłowa Kk=Kkr~0,A~, które zagwarantuje nam spełnienie tezy twierdzenia.

Twierdzenie 8.2

Przy powyższych założeniach, błąd na k-tym kroku metody GCG(B, M, A) spełnia oszacowanie

xk-x*BminpP~kmaxλσMApλx0-x*B.
Dowód

Na mocy wniosku 6.1,

xk-x*B=minpP~kpMAx0-x*B.

Ponieważ dla dowolnego y

MAkyB=B1/2MAB-1/2kB1/2y2,(8.2)

to w konsekwencji pMAyB=pB1/2MAB-1/2B1/2y2pB1/2MAB-1/22yB. Z założenia macierz G=B1/2MAB-1/2 jest normalna, jest więc diagonalizowalna i jej wektory własne są ortogonalne:

G=XΛXT, gdzie XTX=I.

Stąd zaś wynika, że pG2=XpΛXT2=pΛ2=maxλσGpλ, co kończy dowód, gdyż macierze G i MA są podobne, więc mają to samo spektrum.

Ćwiczenie 8.1

Udowodnij równość (8.2).

Można, postępując analogicznie jak w przypadku klasycznej metody CG wykazać, że GCG(B, M, A) ma własności podobne jak jej protoplastka:

Stwierdzenie 8.1

W metodzie GCG(B, M, A) zachodzi:

  • xk=xk-1+αkpk, gdzie pk tworzą bazę B-ortogonalną przestrzeni Kryłowa KkMr0,MA.

  • wTBxk-x*=0 dla każdego wKkMr0,MA.

Dowód

W oczywisty sposób xk=xk-1+dk, gdzie dkKkMr0,MA. Jeśli przez Vk oznaczyć bazę w KkMr0,MA, to xk=x0+Vkak i na mocy (6.2) ak spełnia układ równań normalnych,

VkTBVkak-x*-x0=0.

Ponieważ Vkak=xk-x0, dostajemy

VkTBxk-x*=0,

co dowodzi drugiego punktu.

Metodę GCG(B, M, A) możnaby zrealizować następującym uogólnieniem algorytmu CG, który nazwiemy za [1] algorytmem Odir(B, M, A):

Metoda Odir dla GCG, wersja bazowa
p0=Mr0k=0;
while not stop begin
  αk=pkTBx*-xpkTBpk
  xk=xk-1+αkpk
  rk=rk-1-αkApk
  γk=pkTBMApkpkTBpk
  σk=pkTBMApk-1pk-1TBpk-1
  pk+1=MApk-γkpk-σkpk-1
  k=k+1
end

Algorytm ten jeszcze nie jest gotowy do użycia, albowiem wymaga obliczania Bx*-x — a więc potencjalnie wymaga odwołania się do nieznanego a priori wektora błędu! Możliwość skutecznego obliczania współczynnika αk ogranicza więc zbiór B, których możemy użyć do zdefiniowania konkretnej metody. Z drugiej strony, wybierając konkretne B możemy dalej uprościć powyższą bazową implementację. Na tym etapie nie jest także oczywiste, czy może zdarzyć się pk=0 — a więc stagnacja i załamanie się algorytmu (dzielenie przez zero!).

Wśród metod, które dają się skutecznie zaimplementować na podstawie algorytmu Odir(B, M, A) — to znaczy: mają obliczalne αk — znajdują się metody wymienione w tabeli 8.1, którą przytaczamy za [1, Table 5.2]. Z tabeli wynika m.in., że PCG wcale nie wymaga dodatnio określonej macierzy ściskającej!

Metoda B M Ograniczenia
CG A I A=AT>0
CR A2 I A=AT
PCG A M A=AT>0, M=MT
PCR AMA M A=AT, M=MT>0
CGNR ATA AT
CGNE I AT
D'yakonov B B-1ATB-1 B=BT>0
Tabela 8.1. Wybrane metody Kryłowa dające się zinterpretować jako metoda GCG. W tabeli zamieszczono założenia na macierze A i M gwarantujące poprawność określenia metody.

Aby otrzymać dobrą implementację konkretnej metody opartej na algorytmie Odir, należy zawsze ją dopracować, wykorzystując zależności specyficzne dla konkretnej metody.

Warta uwagi jest metoda PCR, która działa w przypadku, gdy macierz A jest jedynie symetryczna — nie musi być dodatnio określona. Metoda D'aykonova [6] co prawda działa dla dowolnej nieosobliwej, niesymetrycznej macierzy A, ale ze względu na jej podobieństwo do równań normalnych dla B-1A, w wielu wypadkach skuteczniejsza (pod względem szybkości zbieżności, a nie pod względem oszczędności pamięci) będzie metoda typu GMRES.

Uwaga 8.1

Czasami wymaganie, by można było wykonywać mnożenie przez AT może być trudne do spełnienia — na przykład wtedy, gdy A jest zadana jako operator, tzn. wyłącznie przez procedurę mnożenia przez zadany wektor x, tzn. obliczania Ax.

Ćwiczenie 8.2

Porównaj zbieżność metod: PCG i PCR w przypadku, gdy A oraz M są macierzami symetrycznymi i dodatnio określonymi. Przeprowadź weryfikację eksperymentalną swoich przypuszczeń.

Wskazówka: 

W MATLABie dyponujesz gotową implementacją zarówno metody pcg, jak i pcr.

8.1.2. Metody nie oparte na minimalizacji

Zamiast określać kolejną iterację metody, jak dotychczas, przez pewien warunek minimalizacji, możemy inaczej położyć warunek na xk. Na przykład, możemy wymagać by oprócz xkx0+Kkr0,A zachodził warunek ortogonalności residuów:

vTrk=0vKkr0,AT.

(Zwróć uwagę na to, że residua mają być ortogonalne nie do Kkr0,A, tylko do Kkr0,AT.)

Tego typu warunek prowadzi m.in. do metody BiCG (ang. Bi-Conjugate Gradient) która wymaga mnożenia przez AT i nie jest też zbyt stabilna. Jej modyfikacje: CGS i BiCG-stab — nie wymagają już mnożenia przez AT, BiCG-stab jest też bardziej stabilna. Niestety, jak na razie dla BiCG-stab nie ma pełnej i satysfakcjonującej teorii zbieżności.

W przeciwieństwie do metod minimalizacji, kolejna iteracja może nie być realizowalna mimo, że nie osiągnięto jeszcze rozwiązania — mówimy wtedy potocznie o załamaniu się metody (ang. breakdown).

W innej metodzie, QMR (ang. Quasi-Minimal Residual) — i jej wariancie bez mnożenia przez macierz transponowaną: TFQMR (ang. Transpose-Free QMR) — kolejną iterację wybiera się tak, by zminimalizować (stosunkowo łatwą do obliczenia) wielkość, która ma tylko trochę wspólnego z normą residuum.

Można pokazać [9], że metoda TFQMR dla macierzy N×N (realizowana w dokładnej arytmetyce) w ciągu N+1/2 iteracji albo załamie się, albo osiągnie dokładne rozwiązanie.

8.2. Metody strukturalne

W bardzo wielu przypadkach, w macierzach układów równań jakie przychodzi nam rozwiązywać w praktycznych zastosowaniach, nie tylko niewiadome są powiązane jedynie z nielicznymi innymi niewiadomymi (co właśnie prowadzi do rozrzedzenia macierzy) ale w dość naturalny sposób można wskazać grupy niewiadomych i równań, które między sobą nie mają żadnych powiązań.

W najogólniejszym przypadku, mielibyśmy więc

A=A11A13A22A23A31A32A33,(8.3)

gdzie Aii są pewnymi macierzami kwadratowymi. Niewiadome i równania z grupy ,,1” nie mają bezpośredniego związku z grupą ,,2” (są one powiązane ze sobą jedynie poprzez równania i niewiadome grupy ,,3”). W dalszym ciągu będziemy zakładać, że macierze A11 i A22 są nieosobliwe. Strukturalizacja równania liniowego to nic innego jak wskazanie pewnej permutacji równań i niewiadomych, w wyniku której dostajemy układ postaci (8.3). W rzeczywistości, dla macierzy A11 i A22 mogłyby zachodzić podobne relacje10Taką permutację nazywamy wtedy zagnieżdżonym podziałem grafu macierzy (ang. nested dissection)., jak w poniższym przykładzie.

Przykład 8.1 (Strukturalizacja dyskretyzacji równania różniczkowego)

Rozważmy dyskretyzację równania różniczkowego

-Δu=f w Ω

z jednorodnymi warunkami brzegowymi Dirichleta. Jeśli obszar Ω podzielić na kilka mniejszych części, podobszarów, Ω¯=i=1KΩ¯i i przez Γ oznaczyć zestaw krawędzi (ang. interface) je spajających: Γ=(i=1KΩi)-Ω, to — ze względu na lokalny charakter działania pochodnej (i jej sensownych aproksymacji) — niewiadome z dwóch różnych podobszarów Ωi i Ωj nie będą ze sobą powiązane. Będą za to powiązane z niewiadomymi leżącymi na krawędziach spajających Γ, por. rysunek 8.1.

Rys. 8.1. Przykładowy podział obszaru na K=4 podobszary; węzły wspólne (należące do Γ) zaznaczono kolorem niebieskim.

Numerując niewiadome tak, by U=U1,U2,,UK,UΓT, gdzie wektor Ui zawiera niewiadome odpowiadające wewnętrznym węzłom w Ωi, a UΓ — niewiadome odpowiadające węzłom na interfejsie Γ, otrzymalibyśmy macierz o blokowej strukturze

Dh=D11D1ΓD22D2ΓDKKDKΓDΓ1DΓ2DΓKDΓΓN×N.(8.4)

W pustych miejscach Dh znajdują się bloki zerowe, a macierze Dii faktycznie są nieosobliwe.

Dzięki szczególnej postaci macierzy (8.3) możemy podać algorytm typu ,,dziel i rządź” rozwiązywania układu równań Ax=b z tą macierzą. Wprowadzając podział niewiadomych i wektora prawej strony zgodnie z blokowym podziałem macierzy A zadanym (8.3) mamy do rozwiązania układ

A11A13A22A23A31A32A33x1x2x3=b1b2b3.

Dokonując blokowej eliminacji Gaussa (najpierw na blokach A11 i A22) dostajemy następujący algorytm:

  1. Oblicz qi=Aii-1bi, i=1,2.

  2. Wyznacz x3 jako rozwiązanie układu

    Sx3=f,

    gdzie S=A33-A31A11-1A13-A32A22-1A23 oraz f=b3-A31q1-A32q2.

  3. Wyznacz xi=Aii-1bi-Ai3x3, i=1,2.

Oczywiście, wszędzie tam, gdzie występują odwrotności macierzy, w praktyce będziemy rozwiązywali układy równań z tymi macierzami.

Zwróćmy uwagę na kilka charakterystycznych cech tego algorytmu:

  • Zarówno w kroku 1. jak i w kroku 3. algorytmu, oba układy możemy rozwiązywać niezależnie od siebie (co daje możliwość równoległej implementacji).

  • Macierze Aij są rozrzedzone (bo A taka była), ale w ogólności S jest macierzą gęstszą.

Jeśli wymiar macierzy S jest niezbyt wielki, możemy ją wyznaczyć explicite i rozwiązać metodą bezpośrednią. Jednak jeśli S wciąż jest zbyt wielka (lub zbyt gęsta), by potraktować ją metodą bezpośrednią — wtedy możemy zastosować do niej metodę iteracyjną operatorową (matrix-free). Rzeczywiście, mnożenie macierzy S przez wektor z można zrealizować bez wyznaczania niej samej:

Sz=A33z-A31A11-1A13z-A32A22-1A23z.

(Nawiasy wskazują kolejność mnożenia.) Możemy więc wyznaczyć iloczyn Sz kosztem rozwiązania układów równań z macierzami A11 i A22 (co i tak musielibyśmy umieć zrobić!) oraz niewielkim dodatkowym kosztem pomnożenia pewnych wektorów przez rozrzedzone macierze trzeciej kolumny i trzeciego wiersza A.

Ćwiczenie 8.3

Przedyskutuj, jak rozwiązywać układ równań postaci (8.4).

8.3. Metody wielosiatkowe

Macierze rozrzedzone, które otrzymuje się z dyskretyzacji równań różniczkowych cząstkowych metodą różnic skończonych lub elementu skończonego, mają wiele specjalnych cech, które pozwalają na konstrukcję optymalnych metod rozwiązywania takich zadań. Jedną z takich klas metod są metody wielosiatkowe (ang. multigrid).

Omówimy jej ideę na przykładzie macierzy TN dyskretyzacji jednowymiarowego laplasjanu (5.5): wszystkie istotne wady metod iteracyjnych ujawniają się już w jej przypadku, podobnie jak zasadnicze cechy metody wielosiatkowej. Oczywiście, w praktyce metodę wielosiatkową stosowalibyśmy dopiero do dyskretyzacji dwu- lub trójwymiarowego laplasjanu.

Gdyby przyjrzeć się metodzie Jacobiego dla TN zobaczymy, że niektóre składowe błędu są w niej szybciej redukowane niż inne. Dokładniej, rozważmy metodę Jacobiego z parametrem relaksacyjnym ω,

xk+1=xk+ω2rk

(w macierzy TN diagonala jest stała równa 2I). Macierz iteracji tej metody ma postać

G=I-ω2TN.

Wyrażając błąd ek=xk-x* w bazie wektorów własnych macierzy G (są dane wzorem vij=sinijπN+1 — por. (5.18)) dostajemy, że

e0=iaivi

i w konsekwencji,

ek=Gke0=iλiGkaivi.

Jasne jest więc, że najszybszej redukcji ulegną te składowe błędu, które odpowiadają najbliższym zera wartościom λiG, najwolniej zaś będą znikać te, dla których λiG będzie duży. Rysunek 8.2 przedstawia wykres

λiG=1-ω1-cosiπN+1

w zależności od parametru relaksacyjnego ω.

Rys. 8.2. Wartości własne macierzy iteracji metody Jacobiego z parametrem ω1,34,12,14, zastosowanej do macierzy T100.

Jak możemy zauważyć, gdy ω=1 (tzn. gdy używamy standardowej metody Jacobiego), najsłabiej są redukowane składowe szybkozmienne i wolnozmienne (jedynie te ,,średniozmienne” są redukowane błyskawicznie).

Ale już dla ω=1/2, najmocniej są redukowane wszystkie szybkozmienne składowe błędu — zatem po kilku iteracjach takiej metody błąd zostanie wygładzony. Jeśliby więc wspomóc taką prostą iterację poprawką, która redukowałaby pozostałe — wolnozmienne — składowe błędu, moglibyśmy otrzymać szybko zbieżną iterację. Jak pamiętamy z wcześniejszego wykładu, metody projekcji zastępowały idealną poprawkę Aδ*=rk poprawką przybliżoną δ, wyznaczaną przez rozwiązanie zadania zmniejszonego wymiaru. W metodzie dwusiatkowej poprawkę wyznaczymy podobnie: rozwiązując metodą bezpośrednią (np. eliminacji Gaussa) równanie poprawki na rzadszej siatce — wszak będzie ono tylko gorszą aproksymacją (o mniejszej rozdzielczości!) tego samego rozwiązania równania Poissona, które ma aproksymować x*.

W metodzie wielosiatkowej zadanie na rzadszej siatce (potrzebne do wyznaczenia poprawki na siatce najdrobniejszej) potraktujemy rekurencyjnie w ten sam sposób, otrzymując sekwencję zadań coraz mniejszego rozmiaru, którą w pewnym momencie oczywiście będziemy musieli przerwać i rozwiązać zadanie bardzo małego wyniaru metodą bezpośrednią.

Oto więc kolejne etapy jednej iteracji metody wielosiatkowej, wyznaczającej kolejne przybliżenie xk+1 na podstawie xk:

  1. Wygładzenie (ang. pre-smoothing) — zaczynając od xk, wykonujemy kilka iteracji (zwykle jedną lub dwie) prostej metody typu Jacobiego lub Gaussa–Seidela z właściwie dobranym parametrem ω. Dostajemy przybliżenie pośrednie xk+11, ze zredukowaną szybkozmienną składową błędu.

  2. Restrykcja — sformułowanie na rzadszej siatce zadania na poprawkę δc dla xk+11.

  3. Gruba poprawka (ang. coarse grid correction) — wyznaczenie poprawki δc na rzadszej siatce — przez iterację taką, jak ta! (Jako przybliżenie początkowe dla δc sensownie wziąć zero.)

  4. Przedłużenie — wyrażenie poprawki obliczonej na rzadszej siatce w terminach wartości na gęstej siatki. Uwzględnienie poprawki i aktualizacja przybliżenia: xk+12=xk+11+δ.

  5. Dodatkowe wygładzenie (ang. post-smoothing) — ewentualnie. Znów kilka iteracji jak na pierwszym etapie, z przybliżeniem początkowym równym xk+12, zakończone wyznaczeniem xk+1. Gdy ten etap jest pominięty, za xk+1 bierze się xk+12.

Okazuje się, że w ten sposób uzyskujemy metodę, której szybkość zbieżności na modelowym zadaniu nie zależy od parametru dyskretyzacji h (a więc i od rozmiaru zadania). Przypomnijmy, że było to zmorą wszystkich dotychczas rozważanych metod iteracyjnych; aby to przezwyciężyć, musieliśmy tam użyć spektralnie równoważnej macierzy ściskającej. Metoda wielosiatkowa nie potrzebuje tego, co więcej, można pokazać, że koszt jednej iteracji metody jest rzędu ON — a więc (z dokładnością do stałej…) taki sam, jak koszt np. jednej iteracji metody CG.

Pełna metoda wielosiatkowa

Idąc dalej tym tropem możemy zauważyć, że jako dobre przybliżenie początkowe dla (iteracyjnej) metody wielosiatkowej można byłoby wybrać przybliżone rozwiązanie na rzadszej siatce: wtedy jedynymi istotnymi składowymi błędu na gęstej siatce byłyby składowe szybkozmienne — a te wyeliminowalibyśmy w kilku iteracjach! Oczywiście, wymagane przybliżenie na rzadszej siatce uzyskalibyśmy w ten sam sposób z jeszcze rzadszej siatki; na siatce o największych oczkach moglibyśmy użyć jakiejś metody bezpośredniej.

Pamiętając, że rozwiązywane przez nas zagadnienie — układ równań liniowych — pochodzi z dyskretyzacji pewnego równania różniczkowego musimy zdać sobie sprawę z tego, że nie ma sensu rozwiązywać go z dokładnością większą niż dokładność aproksymacji samej dyskretyzacji! Dlatego iterację można przerwać po osiągnięciu z góry zadanego poziomu błędu aproksymacji rozwiązania zadania różniczkowego. Gdy będziemy startować z przybliżenia wyznaczonego metodą opisaną powyżej, okazuje się, że często wystarczy tylko niewielka, niezależna od h, liczba iteracji by osiągnąć zadany poziom błędu.

Taka metoda wyznaczania rozwiązania zadania powstałego po dyskretyzacji równania różniczkowego — którą tutaj tylko zasygnalizowaliśmy w uproszczeniu — nazywa się pełną metodą wielosiatkową (ang. full multigrid). Choć wciąż jest to metoda iteracyjna, ma ona cechy metody bezpośredniej. Okazuje się, że koszt wyznaczenia rozwiązania tą metodą jest rzędu ON (kilku iteracji metody wielosiatkowej na finalnej, najdrobniejszej siatce) — a więc optymalny co do rzędu, gdyż samych niewiadomych w zadaniu jest N.

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.