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 zadań jest rozwiązywanie układów równań liniowych

$Ax=b,$

gdzie $A$ jest nieosobliwą macierzą $N\times N$, a dany wektor prawej strony $b\in R^N$.

W praktyce spotyka się zadania z $N=2,3,\dots 1000$. Zdarzają się także czaem specjalne macierze o wymiarach nawet rzędu $10^8$!

Rozwiązywanie układów równań liniowych jest sercem wielu innych algorytmów numerycznych, dlatego nie dziwi, że 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, takie jak:

• metoda wyznacznikowa (wzory Cramera)
• obliczenie macierzy $A^{-1}$ i następnie $x=A^{-1}b$

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

O tym, jak skutecznie rozwiązywać takie zadania, jakie jest ich uwarunkowanie --- o tym traktują następne dwa wykłady. Okaże się, że jedną z dobrych metod jest metoda eliminacji Gaussa.

Proste układy równań

Niektóre układy równań można bardzo łatwo rozwiązać. Zgodnie z zasadą

Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań

w dalszej kolejności pokażemy, jak dowolny układ równań sprowadzić do sekwencji dwóch (czasem trzech) łatwych do rozwiązania układów równań. Ale... jakie układy równań są "łatwe"?

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

Rozważmy układ z macierzą trójkątną $A$. Będą nas szczególnie interesować macierze trójkątne górne, dla których $a_{i,j}=0$ gdy $i>j$, oraz macierze trójkątne dolne z jedynkami na przekątnej, tzn. $a_{i,j}=0$, $i<j$, oraz $a_{i,i}=1$. Macierze pierwszego rodzaju będziemy oznaczać przez $U$, a drugiego rodzaju przez $L$.

$L=\left(\begin{array}{cccccc}1& & & & & \\ *& 1& & & & \\ *& *& 1& & & \\ *& *& *& 1& & \\ ⋮& ⋮& ⋮& \ddots & \ddots & \\ *& *& *& \cdots & *& 1\end{array}\right),U=\left(\begin{array}{cccccc}*& *& *& *& \cdots & *\\ & *& *& *& \cdots & *\\ & & *& *& \cdots & *\\ & & & *& \ddots & ⋮\\ & & & & \ddots & *\\ & & & & & *\end{array}\right)$

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

$Ux=c,$

$U=(u_{i,j})$, $c=(c_j)$, można rozwiązać stosując algorytm:

Algorytm Podstawienie w tył

---

<math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
for (i = N-1; i >= 1; i--)
	<math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
	u_{i,j}x_j^*\right)/u_{i,i}</math>;

(Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy implikuje, że $u_{i,i}\ne 0$, $\mathrm{\forall }i$.) Podobnie, układ $Lx=c$ rozwiązujemy algorytmem:

Algorytm Podstawienie w przód

---

<math>\displaystyle x_1 = c_1</math>;
for (i=2; i <= N; i++)
	<math>\displaystyle x_i = c_i\,-\,\sum_{j=1}^{i-1} l_{i,j} x_j^*</math>;

Oba algorytmy wymagają rzędu $N^2/2$ mnożeń lub dzieleń i $N^2/2$ dodawań lub odejmowań, a więc łącznie $O(N^2)$ działań arytmetycznych.

Układy z macierzą ortogonalną

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

$Qx=b,$

gdy $Q$ jest macierzą ortogonalną, to znaczy $Q^{T}Q=I$. Rzeczywiście, z ortogonalności mamy natychmiast, że

$x=Q^{T}b$

i w konsekwencji $x$ można wyznaczyć kosztem takim jak koszt mnożenia macierzy przez wektor, czyli $O(N^2)$ operacji.

Podobnie, gdy $Q\in C^{N\times N}$ jest unitarna, to znaczy $Q^{H}Q=I$, rozwiązaniem układu równań jest

$x=Q^{H}b.$

Metoda eliminacji Gaussa

Carl Friedrich Gauss
Zobacz biografię

W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym rozwiązywania układu równań

$Ax=b$

okazuje się popularna eliminacja Gaussa. Jednak z powodów, które będą dla nas później jasne, algorytm ten wyrazimy w języku tzw. rozkładu LU macierzy, to znaczy, sprowadzającego zadanie do znalezienia macierzy trójkątnej dolnej $L$ (z jedynkami na diagonali) oraz trójkątnej górnej $U$ takich, że

$A=LU,$

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

Algorytm Rozwiązanie układu z wykorzystaniem rozkładu LU

---

Znajdź rozkład <math>\displaystyle A=LU</math>;
Rozwiąż <math>\displaystyle Ly = b</math> przez podstawienie w przód;
Rozwiąż <math>\displaystyle Ux = y</math> przez podstawienie w tył;

Przypuśćmy, że taki rozkład $A=LU$ istnieje. Wówczas, zapisując macierze w postaci blokowej, eksponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy, mamy

$\left(\begin{array}{cc}{a}_{11}& a_{12}^T\\ a_{21}& A_{22}\end{array}\right)=\left(\begin{array}{cc}1& 0^T\\ l_{21}& L_{22}\end{array}\right)\left(\begin{array}{cc}{u}_{11}& u_{12}^T\\ 0& U_{22},\end{array}\right)$

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

• $u_{11}=a_{11}$ oraz $u_{12}=a_{12}$, więc pierwszy wiersz $U$ jest

kopią pierwszego wiersza $A$,

• $l_{21}=a_{21}/u_{11}$, więc pierwsza kolumna $L$ powstaje przez

podzielenie wszystkich elementów wektora $a_{21}$ przez element na diagonali $a_{11}$,

• $A_{22}-l_{21}{u}_{12}^T=L_{22}{U}_{22}$, a więc znalezienie podmacierzy

$L_{22}$ oraz $U_{22}$ sprowadza się do znalezienia rozkładu LU zmodyfikowanego bloku $A_{22}$ macierzy $A$, wymiaru $(N-1)\times (N-1)$.

Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z klasycznych pętli. 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ć w miejscu, nadpisując elementy $A$ elementami macierzy $U$ i $L$ (jedynek z diagonali $L$ nie musimy pamiętać, bo wiemy a priori, że tam są).

Algorytm Rozkład LU metodą eliminacji Gaussa

---

for k=1:N-1
	if <math>\displaystyle a_{kk}</math> == 0 
		STOP;
	end
	for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
		<math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
	end
	for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
		for i=k+1:N 
			<math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
		end
	end
end

Łatwo przekonać się, że $k$-ty obrót zewnętrznej pętli (tzn. $k$-ty krok algorytmu rozkładu LU) kosztuje rzędu $2(N-k{)}^2$ operacji arytmetycznych, skąd łączny koszt tego algorytmu rozkładu LU wynosi około $\frac{4}{3}{N}^3$.

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

• Koszt znalezienia rozkładu $A=LU$: $O(N^3)$;
• Koszt rozwiązania układu $Ly=b$: $O(N^2)$;
• Koszt rozwiązania układu $Ux=y$: $O(N^2)$.

Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania wynosi już tylko $O(N^2)$.

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

Z powyższego wynika, że łączny koszt rozwiązania równania liniowego wynosi $O(N^3)$. 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 $N$ równań liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch macierzy $N\times N$. Tymczasem znany jest całkiem prosty algorytm rekurencyjny, wyznaczający iloczyn dwóch macierzy kosztem $4.7\cdot N^{lo{g}_{2}7}\approx 4.7\cdot N^{2.807}$ (algorytm Strassena). Bardziej skomplikowany (i praktycznie nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt $O(N^{2.376})$. Tak więc równania liniowe daje się (teoretycznie) rozwiązać kosztem $O(N^{2.376})$.

Jednak w praktyce nawet prosty algorytm Strassena zazwyczaj nie jest stosowany. Wynika to stąd, że ma trochę gorsze własności numeryczne oraz, co istotniejsze, wymaga dużo dodatkowej pamięci na przechowywanie pośrednich wyników.

Wybór elementu głównego

Opisany powyżej algorytm rozkładu LU niestety czasem może się załamać: gdy napotka w czasie działania zerowy element na diagonali zmodyfikowanej podmacierzy, np. chociaż macierz

$A=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$

jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od razu zetknie się z dzieleniem przez $a_{11}=0$. Ale wystarczy zamienić ze sobą kolejnością wiersze macierzy $A$ (to znaczy, w układzie równań, zamienić ze sobą miejscami równania), a dostajemy macierz, dla której algorytm przejdzie bez problemu.

W praktyce obliczeniowej, aby uzyskać algorytm o [[|Dodaj link: możliwie dobrych własnościach numerycznych]], wykorzystujemy tzw. strategię wyboru elementu głównego w kolumnie. Polega to na tym, że zanim wykonamy $k$-ty krok algorytmu rozkładu LU,

• szukamy w pierwszej kolumnie podmacierzy $A(k:N,k:N)$ 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 $A(k,1:N)$ 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

$PA=LU,$

gdzie $P$ 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 kolumnie oraz tzw. wybór pełny, gdy elementu głównego szukamy w całej podmacierzy $A(k:N,k:N)$, 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 całą informację o wykonanych permutacjach przechowujemy nie w zerojedynkowej macierzy j.w. ale w jednym wektorze.

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 <math>\displaystyle a_{pk}</math>;
	zamień ze sobą wiersze A(k,1:N) i A(p,1:N);
	P(k) = p; P(p) = k;
	if <math>\displaystyle a_{kk}</math>
		STOP: macierz osobliwa!
	end
	
	/* kontunuuj tak jak w algorytmie bez wyboru */
	
	for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
		<math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
	end
	for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
		for i=k+1:N 
			<math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
		end
	end
end

<flash>file=Macierz.swf</flash><div.thumbcaption>Eliminacja 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 <math>\displaystyle PA = LU</math>;
rozwiąż względem <math>\displaystyle y</math> układ z macierzą górną trójkątną <math>\displaystyle Uy = Pb</math>;
rozwiąż względem <math>\displaystyle x</math> układ z macierzą dolną trójkątną <math>\displaystyle Lx = y</math>;

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: $A=A^T$ oraz $x^{T}Ax>0$, $\mathrm{\forall }x\ne 0$,
• silnie diagonalnie dominujących: macierz $A$ (lub $A^T$) spełnia

$|a_{ii}|>\sum _{j\ne i}|a_{ij}|,\mathrm{\forall }i.$