MN05

Układy równań liniowych

Czasem zadania obliczeniowe wymagają wykonania naprawdę wielkiej liczby obliczeń zmiennoprzecinkowych, przy czym wiele znanych, matematycznie równoważnych metod rozwiązywania takich zadań, ma diametralnie różne własności numeryczne. Bardzo ważną klasą takich intensywnych obliczeniowo zadań jest rozwiązywanie układów równań liniowych

,

gdzie
jest nieosobliwą macierzą
, a
jest dany wektor.

Szacuje się, że około 75 procent czasu obliczeniowego superkomputerów jest wykorzystywanych właśnie na rozwiązywanie takich zadań.

Okazuje się, że kilka znanych w matematyce sposobów rozwiązywania układów równań liniowych, takich jak:

• metoda wyznacznikowa (wzory Cramera);

• obliczenie macierzy
  i następnie
  ,

prawie nie nadaje się do numerycznego rozwiązywania takich zadań.

Proste układy równań

Układy z macierzą trójkątną

Rozważmy układ z macierzą trójkątną
. Będą nas szczególnie interesować macierze trójkątne górne
, dla których
, gdy
, oraz macierze trójkątne dolne z jedynkami na przekątnej
, tzn.
(
) i
, mianowicie

i
.

Układ z nieosobliwą macierzą trójkątną górną

,

gdzie
,
, można rozwiązać stosując algorytm:

Algorytm Podstawienie w tył

;

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

;

Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy implikuje, że
.

Podobnie, układ
rozwiązujemy algorytmem:

Algorytm Podstawienie w przód

;

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

;

Oba algorytmy wymagają rzędu
mnożeń lub dzieleń i
dodawań lub odejmowań, a więc łącznie
działań arytmetycznych.

Układy z macierzą ortogonalną

Równie prosto można rozwiązać układ równań

,

gdy
jest macierzą ortogonalną, to znaczy
. Z ortogonalności wynika wprost, że

i w konsekwencji
można wyznaczyć kosztem takim, jak koszt mnożenia macierzy przez wektor, czyli
operacji.

Podobnie, gdy
jest unitarna, to znaczy
(przypomnijmy:
oznacza macierz sprzężoną do
, tzn. taką, że
), rozwiązaniem układu równań jest

.

Metoda eliminacji Gaussa

Carl Friedrich Gauss

Algorytm numerycznego rozwiązywania układu równań metodą eliminacja Gaussa wyrazimy w terminach tzw. rozkładu LU macierzy, to znaczy sprowadzającego zadanie do znalezienia macierzy trójkątnej dolnej
(z jedynkami na diagonali) oraz trójkątnej górnej
takich, że

,

a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:

Algorytm Rozwiązanie układu równań z wykorzystaniem rozkładu LU

Znajdź rozkład
;

Rozwiąż
przez podstawienie w przód;

Rozwiąż
przez podstawienie w tył;

Przypuśćmy, że taki rozkład
istnieje (nie musi!). Wówczas, zapisując macierze w postaci blokowej, eksponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy, mamy

skąd (mnożąc blokowo macierz
przez
) wynika, że

• 
  oraz
  , więc pierwszy wiersz
  jest kopią pierwszego wiersza
  ,

• 
  , więc pierwsza kolumna
  powstaje przez podzielenie wszystkich elementów wektora
  przez element na diagonali
  ,

• 
  , a więc znalezienie podmacierzy
  oraz
  sprowadza się do znalezienia rozkładu LU zmodyfikowanego bloku
  macierzy
  , wymiaru
  . Macierz
  nazywamy uzupełnieniem Schura.

Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie rekurencyjne następuje na końcu, można je zastąpić pętlą. Jest to ważne w praktyce numerycznej, gdyż rekurencja kosztuje: zarówno pamięć, jak i czas.

Ponadto zauważmy, że opisany algorytm możemy wykonać in situ (w miejscu), nadpisując elementy
elementami macierzy
i
(jedynek z diagonali
nie musimy pamiętać, bo wiemy a priori, że tam są). Dzięki temu dodatkowo zaoszczędzimy pamięć.

Algorytm Rozkład LU metodą eliminacji Gaussa - wersja prymitywna

for k=1:N-1

if
== 0

STOP;

end

for i=k+1:N /* wyznaczenie
-tej kolumny
*/

=
;

end

for j=k+1:N /* modyfikacja podmacierzy
*/

for i=k+1:N

-=
;

end

end

end

Łatwo przekonać się, że
-ty obrót zewnętrznej pętli (tzn.
-ty krok algorytmu rozkładu LU) kosztuje rzędu
operacji arytmetycznych, skąd łączny koszt tego algorytmu rozkładu LU wynosi około
.

Jeśli więc do rozwiązywania układu równań
wykorzystamy rozkład LU, to mamy następujące zestawienie kosztów:

• Koszt znalezienia rozkładu
  :
  ;

• Koszt rozwiązania układu
  :
  ;

• Koszt rozwiązania układu
  :
  .

Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania wynosi już tylko
.

Wybór elementu głównego

Opisany powyżej algorytm rozkładu LU czasem może się niestety załamać: mianowicie wtedy, gdy napotka w czasie działania zerowy element w lewym górnym rogu zmodyfikowanej podmacierzy. Na przykład, macierz

jest ewidentnie nieosobliwa, ale nasz algorytm nawet nie ruszy z miejsca, bo od razu zetknie się z dzieleniem przez
... Ale wystarczy zamienić ze sobą wiersze macierzy
(to znaczy, w układzie równań, zamienić kolejność równań), a dostaniemy macierz, z którą nasz algorytm poradzi sobie bez problemu. Musimy więc --- aby stosować eliminację Gaussa do dowolnych macierzy nieosobliwych --- dokonywać takich permutacji równań, by elementem, przez który dzielimy, była zawsze liczba niezerowa (jest to możliwe, na mocy założenia nieosobliwości macierzy).

W praktyce obliczeniowej, aby uzyskać algorytm wykorzystujemy tzw. strategię wyboru elementu głównego w kolumnie. Polega to na tym, że zanim wykonamy
-ty krok algorytmu rozkładu LU,

• w pierwszej kolumnie podmacierzy
  szukamy elementu o największym module (taki element, na mocy założenia nieosobliwości macierzy, jest niezerowy) --- to jest właśnie element główny

• zamieniamy ze sobą wiersz
  z wierszem, w którym znajduje się element główny

• zapamiętujemy dokonaną permutację, bo potem --- gdy już przyjdzie do rozwiązywania układu równań --- będziemy musieli dokonać analogicznej permutacji wektora prawej strony

Wynikiem takiego algorytmu jest więc rozkład

gdzie
jest pewną (zerojedynkową) macierzą permutacji (tzn. macierzą identyczności z przepermutowanymi wierszami).

Oprócz wyboru elementu głównego w kolumnie, stosuje się czasem inne strategie, m.in. wybór w wierszu (analogicznie) oraz tzw. wybór pełny, gdy elementu głównego szukamy w całej podmacierzy
, co znacznie zwiększa liczbę porównań niezbędnych do wskazania elementu głównego, ale też trochę poprawia własności numeryczne takiego algorytmu.

W praktyce, do przechowywania całej informacji o wykonanych permutacjach wystarcza pojedynczy wektor.

Algorytm Rozkład LU z wyborem elementu głównego w kolumnie

P = 1:N; /* tu zapamiętujemy wykonane permutacje */

for k=1:N-1

w wektorze A(k:N,k) znajdź element główny
;

zamień ze sobą wiersze A(k,1:N) i A(p,1:N);

P(k) = p; P(p) = k;

if

STOP: macierz osobliwa!

end

/* kontunuuj tak jak w algorytmie bez wyboru */

for i=k+1:N /* wyznaczenie
-tej kolumny
*/

=
;

end

for j=k+1:N /* modyfikacja podmacierzy
*/

for i=k+1:N

-=
;

end

end

end

Przykład

Przebieg eliminacji Gaussa z wyborem elementu głównego

Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest już rozkład LU wykonany z wyborem w kolumnie.

Algorytm Rozwiązywanie układu równań metodą eliminacji Gaussa z wyborem elementu głównego w kolumnie

znajdź rozkład
;

rozwiąż względem
układ z macierzą górną trójkątną
;

rozwiąż względem
układ z macierzą dolną trójkątną
;

Dla niektórych ważnych klas macierzy wiadomo, że rozkład LU jest wykonalny bez wyznaczania elementu głównego, co istotnie może zmniejszyć całkowity czas działania algorytmu. Jest tak m.in. dla macierzy

• symetrycznych, dodatnio określonych:
  oraz jednocześnie
  ,
  ,

• silnie diagonalnie dominujących: macierz
  (lub
  ) spełnia

Złożoność obliczeniowa zadania rozwiązania układu równań liniowych

Z powyższego wynika, że łączny koszt rozwiązania równania liniowego poprzez rozkład LU wynosi
. Można zastanawiać się, jaka jest najmniejsza możliwa liczba operacji zmiennoprzecinkowych potrzebnych do rozwiązania układu równań liniowych.

Można pokazać, że minimalny koszt rozwiązania układu
równań liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch macierzy
. Tymczasem znany jest całkiem prosty algorytm rekurencyjny, wyznaczający iloczyn dwóch macierzy kosztem
(algorytm Strassena). Bardziej skomplikowany (i praktycznie nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt
. Równania liniowe daje się więc (teoretycznie) rozwiązać kosztem
.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 4 w

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

---