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
ma tylko
niezerowych elementów. Wykorzystanie 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ą
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. 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.

Przykład: Macierz z kolekcji Boeinga

Struktura niezerowych elementów macierzy bcsstk38.

Spójrzmy na macierz sztywności dla modelu silnika lotniczego, wygenerowaną swego czasu w zakładach Boeinga i pochodzącą z dyskretyzacji pewnego równania różniczkowego cząstkowego. Pochodzi z kolekcji Tima Davisa. Jest to mała macierz, wymiaru 8032 (w kolekcji spotkasz równania z milionem i więcej niewiadomych). Jej współczynnik wypełnienia (to znaczy, stosunek liczby niezerowych do wszystkich elementów macierzy) wynosi jedynie
, a więc macierz jest bardzo rozrzedzona: w każdym wierszu są średnio tylko 44 niezerowe elementy.

Reprezentacja 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

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

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

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

Jednak wewnętrznie przechowywane są w innym formacie.

Przykład. Pokażemy jak w Octave wprowadzić macierz rozrzedzoną.

<realnowiki><realnowiki><realnowiki><realnowiki>octave:10> I = [1, 1, 1, 2, 3, 1, 4];

octave:11> J = [1, 3, 2, 2, 3, 4, 4];

octave:12> V = [1, 2, 3, 4, 5, 6, 7];

octave:13> N = 4;

octave:14> A = sparse(I,J,V,N,N)

A =

Compressed Column Sparse (rows = 4, cols = 4, nnz = 7)

(1, 1) -> 1

(1, 2) -> 3

(2, 2) -> 4

(1, 3) -> 2

(3, 3) -> 5

(1, 4) -> 6

(4, 4) -> 7

octave:15> spy(A);

</realnowiki></realnowiki></realnowiki></realnowiki>

Strukturę jej niezerowych elementów ukaże nam polecenie spy(A). Tak właśnie zostały wygenerowane obrazki macierzy w niniejszym wykładzie.

Format spakowanych kolumn (wierszy)

Format współrzędnych nie narzucał żadnego uporządkowania elementów macierzy oraz 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, dobrze 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
, zawierającym kolejne niezerowe elementy macierzy
wpisywane kolumnami, AI jest wektorem typu int o długości
, zawierającym numery wierszy macierzy
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
, zawierający na
-tym miejscu indeks pozycji w AV, od którego rozpoczynają się w AV elementy
-tej kolumny macierzy
.

Format CSC

Mamy więc zależność, przy założeniu, że
,

Taki (z drobnymi modyfikacjami) format macierzy wykorzystują np. pakiety ARPACK i UMFPACK (a także, wewnętrznie, Octave i MATLAB).

Format diagonalny

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

Uwagi praktyczne

Mnożenie macierzy w formacie AIJ przez wektor jest kilka razy wolniejsze w porównaniu do macierzy w formacie CSC/CSR (z tych dwóch oczywiście CSR jest najszybszy). Co więcej, okazuje się, że w typowych implementacjach, mnożenie macierzy rozrzedzonej (reprezentowanej np. w formacie CSC) przez wektor jest mało efektywne, mniej więcej na poziomie dziesięciu procent możliwości obliczeniowych procesora.

Jeśli już poważnie myśleć o przyspieszeniu mnożenia macierzy przez wektor (np. gdy chcemy stosować iteracyjną metodę rozwiązywania układu równań z tą macierzą), warto rozważyć inne formaty --- łączące w sobie podział macierzy na bloki (tak, jak w algorytmach BLAS) i przechowywanie (w zasadzie) tylko niezerowych elementów macierzy.

Macierze specjalne

Zajmiemy się teraz zadaniem rozwiązywania układu równań liniowych

ale w sytuacji, gdy macierz
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
taka, że dla
,
, nazywamy macierzą taśmową z rozstępem
, o szerokości pasma
.

Ł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
--- tak więc, dla niezbyt dużych
wciąż wynikowa macierz jest taśmowa.

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

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

Macierze trójdiagonalne

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

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

jest bardzo często spotykane, więc warto przytoczyć algorytm w całej okazałości --- popularnie zwany algorytmem przeganiania (niektóre źródła nazywają go algorytmem Thomasa).

Twierdzenie. Jeśli macierz
ma słabo dominującą diagonalę, tzn.

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

Algorytm Metoda przeganiania (w miejscu)

for (i = 2; i <= N; i++)

{

=
/
;

=
-
*
;

=
-
*
;

}

=
/
;

for (i = N-1; i >= 1; i--)

= (
-
*
)/
;

Metody bezpośrednie

Przykład: Strzałka Wilkinsona. Rozważmy układ równań z macierzą diagonalnie dominującą postaci

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 macierzy permutacji
):

Wtedy okazuje się, że rozkład naszej macierzy nie powoduje już wypełnienia 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: jak 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 (tzn. macierzy).

W ogólnym przypadku rozwiązanie takiego zadania jest trudne i większość algorytmów opiera się na pewnych heurystykach, które jednak w praktyce warto wspomóc wcześniejszą analizą konkretnego układu równań jaki 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 BLAS Level 3, a więc permutacje wierszy i kolumn macierzy muszą to także brać pod uwagę.

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

• przybliżone algorytmy minimalnego stopnia (approximate minimum degree, AMD), np. AMD;

• techniki podziału grafów na (prawie) rozłączne składowe nested dissection, ND, np. 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: Rozwiązywanie układu z macierzą rozrzedzoną w Octave. Najprostszy sposób na wykorzystanie metody bezpośredniego rozwiązywania układu z macierzą rzadką to zastosowanie znanego nam operatora do macierzy typu rozrzedzonego:

A = sparse(...);

x = A \ b;

Octave domyślnie przyłoży do
reordering colamd i następnie skorzysta z biblioteki UMFPACK, by rozwiązać taki układ. Dodatkowo, badane jest wcześniej, czy macierz nie jest przypadkiem taśmowa o wąskiej wstędze: jeśli jest, to korzysta się z LAPACKa.

Przykład: Wypełnienie pewnej macierzy w zależności od użytego reorderingu. Rozważmy ponownie macierz pochodzącą z kolekcji Tima Davisa. Jest to macierz symetryczna, dodatnio określona, wymiaru 8032, o 355460 elementach niezerowych i, w konsekwencji, o wypełnieniu około 0.6 procent.

Struktura niezerowych elementów macierzy.

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
wykonanego standardowym algorytmem. Czas rozkładu: 0.892013

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

Czynnik rozkładu Cholesky'ego
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
z odwrotnym reorderingiem Cuthill-McKee. Czas rozkładu: 0.845928

Czynnik rozkładu Cholesky'ego
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
(bez reorderingu). Czas rozkładu LU: 1.696018s. Nieakceptowany, podobnie jak drastyczne wypełnienie macierzy.

Rozkład LU, czynnik
(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, dla tej macierzy, która jest specyficznego typu --- pochodzi z dyskretyzacji równania różniczkowego --- algorytm ND mógłby tu jeszcze lepiej zadziałać.

Stacjonarne metody iteracyjne

Gdy macierz jest rozrzedzona, mnożenie takiej macierzy przez wektor jest bardzo 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 w praktyce najmniej efektywne) polegają na rozkładzie macierzy na część "łatwo odwracalną",
, i "resztę",
. Dokładniej, jeśli
jest nieosobliwa, to równanie
można zapisać jako zadanie punktu stałego

gdzie
. Inaczej:

i zastosować doń metodę iteracji prostej Banacha:

Takie metody nazywamy stacjonarnymi metodami iteracyjnymi. Aby przeanalizować zbieżność takiej metody, warto rozpatrzyć przypadek ogólniejszy

dla pewnej macierzy
oraz wektora
. (Dla stacjonarnej metody iteracyjnej,
oraz
). W tym przypadku

a stąd i z nierówności
, mamy

Warunkiem dostatecznym zbieżności iteracji prostych jest więc
. Okazuje się, że warunkiem koniecznym i dostatecznym zbieżności tej iteracji dla dowolnego wektora startowego
jest

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

Zaletą stacjonarnych metod iteracyjnych jest również ich prostota, przez co są one łatwe do zaprogramowania, co łatwo zobaczyć na przykładach metod: Jacobiego i Gaussa--Seidela, które teraz omówimy.

Metoda Jacobiego

Biorąc
, gdzie
jest macierzą diagonalną składającą się z wyrazów stojących na głównej przekątnej macierzy
, układ
jest równoważny układowi

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

gdzie
i
, zwaną metodą Jacobiego.

Rozpisując ją po współrzędnych dostajemy (numer iteracji wyjątkowo zaznaczamy w postaci górnego indeksu) układ rozszczepionych równań:

co znaczy dokładnie tyle, że w
-tym równaniu wyjściowego układu przyjmujemy za współrzędne
wartości z poprzedniej iteracji i na tej podstawie wyznaczamy wartość
.

Twierdzenie O zbieżności metody Jacobiego. W metodzie Jacobiego warunek dostateczny zbieżności,
, jest spełniony np. wtedy, gdy macierz
ma dominującą przekątną, tzn. gdy

Dowód. Rzeczywiście, ponieważ wyraz
macierzy
wynosi
dla
oraz
dla
, a więc

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

Przykład: Macierz laplasjanu. Macierz
, zwana macierzą jednowymiarowego laplasjanu

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 układ równań z tą macierzą można bez trudu rozwiązać metodami bezpośrednimi, kosztem
. Stosując do niej metodę Jacobiego mamy
oraz
. Obliczając normę macierzy iteracji Jacobiego dostajemy
, co nie rozstrzyga jeszcze o jej zbieżności lub niezbieżności. Potrzebna będzie bardziej subtelna analiza. Okazuje się, że są znane wzory na wartości własne macierzy
:

dla
W konsekwencji, wartościami własnymi
są liczby
. Ponieważ
, znaczy to, że metoda Jacobiego jest zbieżna dla macierzy
.

Z drugiej strony, nie dajmy się zwieść optymizmowi matematyka ("nie martw się, jest zbieżny..."): nietrudno sprawdzić, że
, co oznacza, że metoda Jacobiego --- choć zbieżna --- dla dużych
staje się zbieżna tak wolno, że w praktyce bezużyteczna.

Metoda Gaussa-Seidela

Heurystyka tej metody opiera się na zmodyfikowaniu metody Jacobiego tak, by w każdym momencie iteracji korzystać z najbardziej "aktualnych" współrzędnych przybliżenia rozwiązania
.

Rzeczywiście, przecież rozwiązując równanie metody Jacobiego:

nietrudno zauważyć, że w części sumy moglibyśmy odwoływać się do "dokładniejszych" wartości
: dla
, tzn.

W języku rozkładu macierzy
i iteracji
mamy więc
(dolny trójkąt macierzy
).

Twierdzenie O zbieżności metody Gaussa-Seidela. Jeśli macierz
jest diagonalnie dominująca, to metoda Gaussa--Seidela jest zbieżna dla dowolnego wektora startowego
.

Inny wariant tej metody dostalibyśmy, biorąc za
górny trójkąt macierzy
.

Metoda Gaussa--Seidela jest w wielu przypadkach rzeczywiście szybciej zbieżna od metody Jacobiego, np. tak jest w przypadku macierzy jednowymiarowego Laplasjanu. Wciąż jednak, dodajmy, dla zadań bardzo źle uwarunkowanych jej zbieżność jest zbyt wolna by ją stosować jako samodzielną metodę.

Obie metody, Jacobiego i (zwłaszcza) Gaussa--Seidela stosuje się także czasem w prostych algorytmach rozwiązywania układów równań nieliniowych: ich zaletą jest to, że głównym składnikiem iteracji jest rozwiązywanie skalarnego równania nieliniowego na każdym kroku metody.

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
, przez złożoność metody będziemy rozumieli koszt kombinatoryczny obliczenia
z zadaną dokładnością
. Dla uproszczenia założymy, że metoda jest zbieżna liniowo z ilorazem
. Zauważmy, że aby zredukować błąd początkowy do
, wystarczy wykonać
iteracji, gdzie
spełnia

czyli

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

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

Stąd oczywisty wniosek (prawdziwy nie tylko dla metod stacjonarnych), że metody iteracyjne warto stosować zamiast metod bezpośrednich w przypadku, gdy

• wymiar
  układu
  jest "duży", oraz

• macierz
  układu jest "rozrzedzona", tzn. ma stosunkowo niewielką liczbę elementów niezerowych, np. proporcjonalną do
  .

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
dobiera się w taki sposób, by minimalizowało pewną miarę błędu na podprzestrzeni Kryłowa

gdzie
jest residuum na początku iteracji. Przestrzeń Kryłowa jest rozpięta przez kolejne wektory metody potęgowej. 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 najpopularniejszą: CG.

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

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

na przestrzeni afinicznej
. 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 Metoda CG

r = b-A*x;

=
;
= 0; k = 1;

while (!stop)

{

p = r +
*p;

w = A*p;

=
/
;

x = x +
*p;

r = r -
*w;

=
;

=
;

k++;

}

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

Twierdzenie O zbieżności CG jako metody bezpośredniej. Niech
będzie symetryczna i dodatnio określona. Algorytm CG znajdzie dokładne rozwiązanie po co najwyżej
iteracjach.

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

• dla bardzo dużych
  , wykonanie
  iteracji może być wciąż zbyt kosztownym zadaniem;

• ponieważ w arytmetyce skończonej precyzji ortogonalność --- z której w bardzo istotny sposób korzysta się przy wyprowadzeniu algorytmu --- pogarsza się z iteracji na iterację i w konsekwencji, po wielu iteracjach, jakość
  przestaje się poprawiać.

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

Twierdzenie Zbieżność CG jako metody iteracyjnej. Po
iteracjach metody CG,

gdzie
.

GMRES

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. Jej szczegółowe omówienie, w tym --- oszacowania szybkości zbieżności --- wykracza niestety poza ramy niniejszego wykładu.

Ściskanie macierzy

Zbieżność wszystkich poznanych metod iteracyjnych zależy od własności spektralnych macierzy układu. Pojawiające się w zastosowaniach macierze często 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

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

gdzie macierz
ma znacznie korzystniejsze własności z punktu widzenia używanej metody iteracyjnej.

W przypadku macierzy symetrycznych widzieliśmy, że kluczowe znaczenie dla zbieżności metody miał rozkład wartości własnych: 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
  były upakowane w klastrach

• najlepiej wszystkie w (małym) otoczeniu wartości 1

Jeśli więc chcielibyśmy przekształcić macierz tak, by metoda iteracyjna dla
zbiegała szybko, musimy w jakiś sposób "ścisnąć" spektrum macierzy
w okolice jedności. Taką operację nazywamy ściskaniem (preconditioning), a macierz
--- macierzą ściskającą. Aby całość miała sens, macierz ściskająca
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:
  ),

• macierz
  powinna mieć znacznie korzystniejsze własności z punktu widzenia używanej metody iteracyjnej.

Kilka ekstremalnych (lecz wątpliwej jakości) "propozycji" na macierz ściskającą to
(łatwa w konstrukcji i tania w mnożeniu, ale niestety nic nie polepsza...) oraz
(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 (aczkolwiek wciąż nie najbardziej skutecznych) sposobów ściskania są te oparte na zastosowaniu jednego kroku klasycznej metody iteracyjnej.

Zbieżność metody CG bez żadnego ściskania oraz ściśniętej imadłem opartym na jednej iteracji (blokowej) metody Jacobiego.

Inne sposoby ściśnięcia macierzy wykorzystują 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 ściskającą.

Literatura

• D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

---