MN04LAB

Rozwiązywanie układów równań liniowych

Ćwiczenie

W Octave układ równań $Ax=b$ rozwiązujemy korzystając z "operatora rozwiązywania równania", tzn.

 
x = A \ b;

W tym celu Octave oczywiście wewenętrznie wykorzystuje funkcję LAPACKa DGESV. Ale w Octave jest także funkcja inv, wyznaczająca macierz odwrotną, więc niektóre (nie najlepsze, oględnie mówiąc) podręczniki zalecają

 
x = inv(A)*b;

Przedyskutuj, które podejście jest lepsze i dlaczego. Przeprowadź eksperymenty numeryczne weryfikujące Twoją tezę.

Rozwiązanie

Wyznaczenie $A^{-1}$ praktycznie nigdy nie jest konieczne w obliczeniach, a zwłaszcza do rozwiązania układu równań. Nie dość, że jest bardziej kosztowne, to dodatkowo nie jest numerycznie poprawne.

Przykładowy test mógłby wyglądać w poniższy sposób. Duże macierze losowe są zazwyczaj dobrze uwarunkowane i najczęściej nieosobliwe.

 
N = 400; A = rand(N,N);x = rand(N,1); b = A*x; 
tic; X = A \ b; tX = toc; 
tic; Y = inv(A)*b; tY = toc;
[norm(X-x)/norm(x),  norm(Y-x)/norm(x); tX, tY] 

 
ans =
   1.9861e-13   2.8687e-13
   2.3644e-01   4.0609e-01

Wyszło (z grubsza) dwa razy szybciej i dwa razy dokładniej.

Ćwiczenie

Zapisz w Octave algorytm rozkładu LU macierzy (bez wyboru elementu głównego) działający in situ.

Wskazówka

Pamiętaj, pętle w Octave wykonują się bardzo wolno!

Wykorzystaj go do napisania funkcji, która rozwiąże układ równań $Ax=b$.

Przetestuj tę funkcję na kilku macierzach i porównaj czas jego działania z czasem wykonania operacji x = A\b.

Spróbuj zastosować swój algorytm do kilku specjalnych macierzy:

• Hilberta dużego wymiaru
• diagonalnej z jednym elementem bardzo małym (a nawet równym zero)

Rozwiązanie

Najprostszy wariant to po prostu przepisanie algorytmu z wykładu, który już nawet trochę wykorzystywał notację dwukropkową MATLABa:

 
function [L,U] = lufa(A)
# Na początku warto dodać trochę sprawdzania danych: czy A jest kwadratowa? a w ogóle,
# czy jest macierzą? te linie pozostawiamy inwencji Czytelnika

N = size(A,1);
for k=1:N-1
	if (A(k,k) == 0.0)
		return;
	end
	for i=k+1:N #wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math>
		A(i,k) = A(i,k)/A(k,k);
	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 
			A(i,j) -= A(i,k)*A(k,j);
		end
	end
end
L = tril(A,-1) + eye(N);
U = triu(A);
end

Jednak w tej implementacji mamy do czynienia z nawet potójnie zagnieżdżonymi pętlami, co dla nawet średnich $N$ może być powodem sporego spowolnienia algorytmu. Warto więc uniknąć pętli, odwołując się do instrukcji wektorowych i blokowych.

 
function [L,U] = lufb(A)
# Na początku warto dodać trochę sprawdzania danych: czy A jest kwadratowa? a w ogóle,
# czy jest macierzą? te linie pozostawiamy inwencji Czytelnika

N = size(A,1);
for k=1:N-1
	if (A(k,k) == 0.0)
		return;
	end
#wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math>
		A(k+1:N,k) = A(k+1:N,k)/A(k,k);
#modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> 
		A(k+1:N,k+1:N) -= A(k+1:N,k)*A(k,k+1:N);
end
L = tril(A,-1) + eye(N);
U = triu(A);
end

 [Kod testujący lufa(0 i lufb()]
N = 100; 
A = rand(N) + 1000*eye(N);

format short e;

t = ta = tb = Inf;

K = 5;

for i = 0:K
	tic; [L, U] = lu(A); t = min(t,toc);
	tic; [La, Ua] = lufa(A); ta = min(ta,toc);
	tic; [Lb, Ub] = lufb(A); tb = min(tb,toc);
end

disp("Czasy: lu(), lufa(), lufb()"); [t, ta, tb]

disp("Błąd:");
[norm(L*U - A,inf), norm(La*Ua - A,inf), norm(Lb*Ub - A,inf)]

Dla macierzy wymiaru $N=100$ dostaliśmy na Pentium 1.5GHz

Uzupelnij tytul

Implementacja	lu (oryginalne)	lufa	lufb
Czas (s)	4.0436e-03	1.5877e+01	1.0767e-01
Błąd $||LU-A|{|}_{\mathrm{\infty }}$	3.4596e-13	8.0747e-13	8.0747e-13

Zwróć uwagę na charakterystyczne fakty:

• Usunięcie wewnętrznych pętli stukrotnie(!) przyspieszyło lufb w

stosunku do lufa

• Gotowa funkcja Octave wciąż jest stukrotnie(!) szybsza od mojej

najszybszej funkcji.

To jest typowe, gdy korzystamy ze środowisk typu Octave czy MATLAB. Czasem nawet lepiej napisać implementację wykonującą więcej obliczeń, ale korzystającą w większm stopniu z wbudowanych funkcji Octave, niż napisać samemu implementację liczącą mniej, ale korzystającą z licznych instrukcji for... i podobnych.

Innym sposobem przyspieszenia działania funkcji Octave jest ich prekompilacja. To w znacznym stopniu usuwa problemy związanie z korzystaniem z pętli, ale wykracza poza zakres naszego kursu.

Funkcja wykorzystująca gotowy rozkład LU do rozwiązania układu $Ax=b$ w MATLABie miałaby postać

 
function x = sol(L,U,b)
x = U \ ( L \ b);
end

gdyż MATLAB prawidłowo rozpoznaje układy z macierzą trójkątną i stosuje szybszy algorytm. Niestety, do tej pory w Octave nie jest to możliwe; nie ma też funkcji, która rozwiązywałaby układ z macierzą trójkątną (wykorzystując np. LAPACKa). Zrób więc to i... wyślij do twórców Octave'a!

Ćwiczenie

Zrób zadanie poprzednie w C i porównaj z procedurą biblioteczną LAPACKa DGESV, najlepiej wspartą dobrze podrasowanymi BLASami.

Ćwiczenie: Układy równań z wieloma prawymi stronami

Podaj sposób taniego wyznaczenia rozwiązania sekwencji $k<N$ układów równań z tą samą macierzą $A\in R^{N\times N}$ i różnymi prawymi stronami:

$A{x}_i=b_i,i=1,\dots ,k.$

Układy równań z tą samą macierzą, ale ze zmieniającą się prawą stroną równania powstają często przy rozwiązywaniu, np. równań różniczkowych cząstkowych, gdzie prawa strona układu odpowiada zmieniającym się warunkom brzegowym.

Wskazówka

Wystarczy tylko jeden raz wyznaczyć rozkład LU macierzy $A$

Rozwiązanie

No tak, oczywiście, wystarczy

Algorytm

---

Znajdź rozkład <math>\displaystyle PA = LU</math>;
Utwórz macierz <math>\displaystyle B = [b_1,\ldots,b_k]</math>;
Rozwiąż <math>\displaystyle LY = B</math>;
Rozwiąż <math>\displaystyle UX = Y</math>;

Rozwiązanie każdego z układów $LY=B$ i $UX=Y$ można przyspieszyć, tworząc warianty blokowe znanych algorytmów. Zaimplementowano je w LAPACKu, w tandemie funkcji DGETRF (rozkład) i DGETRS (rozwiązanie) oraz w DGESV, łączącej funkcjonalność obu.

Ćwiczenie: Metoda Cholesky'ego

Ważnym przykładem macierzy szczególnej postaci są macierze symetryczne i dodatnio określone. Są to macierze spełniające $A=A^T$ oraz

$x^{T}Ax>0,\mathrm{\forall }x\ne 0.$

Dla takich macierzy można nieco zmniejszyć koszt kombinatoryczny i zużycie pamięci przeprowadzając trochę inny rozkład na macierze trójkątne: tak, aby otrzymać rozkład

$A=L\cdot D\cdot L^T$

zamiast $PA=LU$, przy czym $L$ jest tu jak zwykle macierzą trójkątną dolną z jedynkami na przekątnej, a $D$ jest macierzą diagonalną z dodatnimi elementami na diagonali. Opracuj taki algorytm. W jego implementacji możesz porównywać się z procedurą LAPACKa DPOSV Inny wariant tego samego rozkładu to tak zwany rozkład Cholesky'ego--Banachiewicza, w którym przy tych samych założeniach na $A$, szukamy rozkładu wykorzystującego tylko jedną macierz trójkątną dolną:

$A=\tilde{L}\cdot {\tilde{L}}^T,$

(oczywiście tym razem nie żądamy, aby $\tilde{L}$ miała na diagonali jedynki). Jaka jest relacja między rozkładem $LD{L}^T$ a $\tilde{L}{\tilde{L}}^T$?

Rozwiązanie

Bez zmniejszenia ogólności rozpatrzymy tylko pierwszy krok rozkładu $LD{L}^T$. W tym celu, zauważmy najpierw, że $a_{1,1}=e_1^{T}A{e}_1>0$ (gdzie $e_1$ jest pierwszym wersorem), a więc nie musimy przestawiać wierszy, bo element na diagonali jest niezerowy. W pierwszym kroku mnożymy macierz $A$ z lewej strony przez odpowiednią macierz $L_1$, a potem z prawej przez $L_1^T$. Kluczem do zrozumienia algorytmu jest uwaga, że efektem mnożenia macierzy $L_{1}A$ z prawej strony przez $L_1^T$ jest wyzerowanie elementów pierwszego wiersza poza $a_{1,1}$ i pozostawienie niezmienionych pozostałych elementów. Ponadto macierz $A^{(1)}=L_{1}A{L}_1^T$ jest symetryczna i dodatnio określona. Rzeczywiście,

$(A^{(1)}{)}^T=(L_{1}A{L}_1^T{)}^T=(L_1^T{)}^T{A}^T{L}_1^T=L_{1}A{L}_1^T,$

oraz dla $x\ne 0$

$x^T{A}^{(1)}x=x^T{L}_{1}A{L}_1^{T}x=(L_1^{T}x{)}^{T}A(L_1^{T}x)>0,$

bo $x\ne 0$ implikuje $L_1^{T}x\ne 0$. Stąd $a_{2,2}^{(1)}=e_2^T{A}^{(1)}{e}_2>0$. Postępując tak dalej otrzymujemy

$L_{n-1}{L}_{n-2}\cdots L_2{L}_{1}A{L}_1^T{L}_2^T\cdots L_{n-2}^T{L}_{n-1}^T=D,$

przy czym macierz $D$ jest diagonalna i dodatnio określona, a więc wyrazy na diagonali są dodatnie. Oznaczając

$L=L_1^{-1}{L}_2^{-1}\cdots L_{n-1}^{-1}$

dostajemy żądany rozkład.

Zauważmy, że przy praktycznej realizacji rozkładu $A=LD{L}^T$ wystarczy modyfikować jedynie wyrazy pod i na głównej przekątnej macierzy wyjściowej, ponieważ, jak zauważyliśmy, wszystkie kolejne macierze $A^{(k)}$ są symetryczne. Pozwala to zmniejszyć --- w stosunku do rozkładu LU --- koszt obliczeniowy o połowę do $2n^3/3$ operacji arytmetycznych i zmieścić czynniki rozkładu w pamięci przeznaczonej na przechowywanie istotnych elementów $A$ (czyli w dolnym trójkącie macierzy).

Oczywiście, $\tilde{L}=L\sqrt{D}$.

Ćwiczenie: Obliczanie wyznacznika macierzy

Bardzo rzadko w praktyce numerycznej zdarza się potrzeba obliczenia wartości wyznacznika macierzy $A\in R^{N\times N}$. Zaproponuj metodę obliczania $\text{det}(A)$ oraz wskaż, jakiego rodzaju problemy numeryczne możesz napotkać.

Rozwiązanie

Wystarczy wykonać rozkład $PA=LU$. Wtedy

$\text{det}A=(-1{)}^p{u}_{11}{u}_{22}\cdots u_{nn},$

gdzie $p$ jest liczbą przestawień wierszy w eliminacji. Tak więc, koszt obliczenia wyznacznika to $O(N^3)$.

Kłopoty numeryczne możemy mieć z reprezentacją samej wartości wyznacznika w arytmetyce zmiennoprzecinkowej. Ponieważ

$\text{det}(cA)=c^{N}A,$

to dość łatwo będzie trafić (przy większych niż średnich wartościach $N$) na wartości wyznacznika, które będą albo potwornie duże, albo potwornie małe. Zresztą, zrób eksperyment:

Najpierw z losową macierzą $100\times 100$ i jej krotnościami...

 
octave:26> A=rand(100);
octave:27> det(A)
ans = -7.5240e+24
octave:28> det(10*A)
ans = -7.5240e+124
octave:29> det(100*A)
ans = -7.5240e+224
octave:30> det(1000*A)
ans = -Inf

oraz z coraz większymi macierzami losowymi:

 
octave:32> det(rand(100))
ans =  4.3751e+25
octave:33> det(rand(400))
ans = -3.0348e+218
octave:34> det(rand(800))
ans = -Inf

Można sobie z tym poradzić, szacując wpierw rząd wielkości wyznacznika, a następnie reprezentując jego wartość w formie cecha--mantysa, np. jako parę liczb $m,c$ takich, że $\text{det}(A)=m\cdot 10^c$, ale dopasowując obie tak, by nie nastąpił ani nadmiar, ani niedomiar.