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 K_k(M r₀, M A) = span{M r₀, M AM r₀, . . . , (M A)^k−1M r₀}, a iteracja ma prostą

formu-łę x_k+1= x_k+ d_k, gdzie d_k∈ K_k(M r₀, M A) 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:

kxk_B= (x^TBx)^1/2.

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

GG^T = G^TG,

gdzie G = B^1/2M AB^−1/2.

W uogólnionej metodzie CG, którą będziemy oznaczali GCG(B, M , A), kolejna iteracja

x_k+1 ∈ x₀+ K_k(M r₀, M A) będzie określona tak, by

kx_k+1− x^∗k_B¬ kx − x^∗k_B ∀x ∈ x₀+ K_k(M r₀, M A) (8.1)

Ćwiczenia testowe 8.1. 1. Czy metoda CG to nic innego jak GCG(I, I, A)?

NIE. Komentarz do odpowiedzi prawidłowej: Oczywiście, CG to nic innego jak GCG(A, I, A). Komentarz do odpowiedzi błędnej: Sprawdź, w jakiej normie CG minimalizuje błąd.

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

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

TAK. Komentarz do odpowiedzi prawidłowej: Oczywiście. Komentarz do odpowiedzi błędnej: 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 ˜r0 = M r₀ oraz ˜A = M A, możemy zastosować twierdzenie 6.1 do przestrzeni Kryłowa K_k= K_k(˜r0, ˜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 kx_k− x^∗k_B¬ min p∈ ˜Pk max λ∈σ(M A) |p(λ)|kx₀− x^∗k_B.

Dowód. Na mocy wniosku 6.1,

kx_k− x^∗k_B= min p∈ ˜Pk

kp(M A)(x₀− x^∗)k_B.

Ponieważ dla dowolnego y

k(M A)^kyk_B = k(B^1/2M AB^−1/2)^kB^1/2yk₂, (8.2) to w konsekwencji kp(M A)yk_B = kp(B^1/2M AB^−1/2) B^1/2yk2 ¬ kp(B1/2M AB^−1/2)k₂kyk_B. Z założenia macierz G = B^1/2M AB^−1/2 jest normalna, jest więc diagonalizowalna i jej wektory własne są ortogonalne:

G = XΛX^T, gdzie X^TX = I.

Stąd zaś wynika, że kp(G)k₂ = kXp(Λ)XTk₂ = kp(Λ)k₂ = max_λ∈σ(G)|p(λ)|, co kończy dowód,

gdyż macierze G i M A są podobne, więc mają to samo spektrum.

Ćwiczenie 8.2. 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:

— x_k= x_k−1+α_kpk, gdzie p_k tworzą bazę B-ortogonalną przestrzeni Kryłowa K_k(M r₀, M A). — w^TB(x_k− x^∗) = 0 dla każdego w ∈ K_k(M r₀, M A).

Dowód. W oczywisty sposób x_k= x_k−1+ d_k, gdzie d_k∈ K_k(M r₀, M A). Jeśli przez V_k oznaczyć bazę w K_k(M r₀, M A), to x_k= x₀+ V_ka_k i na mocy (6.2) a_kspełnia układ równań normalnych,

V_k^TB(V_ka_k− (x^∗− x₀)) = 0. Ponieważ V_ka_k= x_k− x₀, dostajemy

V_k^TB(xk− x^∗) = 0, co dowodzi drugiego punktu.

Metoda B M Ograniczenia CG A I A = A^T > 0 CR A² I A = A^T PCG A M A = A^T > 0, M = M^T PCR AM A M A = A^T, M = M^T > 0 CGNR A^TA A^T CGNE I A^T D’yakonov B B⁻¹A^TB⁻¹ B = B^T > 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.

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

p₀= M r₀; k = 0;

while not stop begin

α_k =^p T kB(x^∗− x) pT kBpk xk= xk−1+ αkpk rk = rk−1− αkApk γk =^p T kBM Apk pT kBpk σk= ^p T kBM Apk−1 pT k−1Bp_k−1 pk+1= M Apk− γkpk− σkpk−1 k = k + 1 end

Algorytm ten jeszcze nie jest gotowy do użycia, albowiem wymaga obliczania B(x^∗− 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 zde-finiowania 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ę

p_k= 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!

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 sy-metryczna — 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⁻¹A, w wielu wypadkach skuteczniejsza (pod względem szybkości

98 8. Przegląd innych metod rozwiązywania wielkich układów równań liniowych Uwaga 8.1. Czasami wymaganie, by można było wykonywać mnożenie przez A^T 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 A · x.

Ćwiczenie 8.3. 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 ipcr.

8.1.2. Metody nie oparte na minimalizacji

Zamiast określać kolejną iterację metody, jak dotychczas, przez pewien warunek minima-lizacji, możemy inaczej położyć warunek na x_k. Na przykład, możemy wymagać by oprócz

x_k ∈ x₀+ K_k(r₀, A) zachodził warunek ortogonalności residuów: v^Trk= 0 ∀v ∈ K_k(r₀, A^T).

(Zwróć uwagę na to, że residua mają być ortogonalne nie do K_k(r₀, A), tylko do K_k(r₀, A^T).) Tego typu warunek prowadzi m.in. do metody BiCG (ang. Bi-Conjugate Gradient) która wymaga mnożenia przez A^T i nie jest też zbyt stabilna. Jej modyfikacje: CGS i BiCG-stab — nie wymagają już mnożenia przez A^T, 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 d(N + 1)/2e iteracji albo załamie się, albo osiągnie dokładne rozwiązanie.