MN05LAB

Ćwiczenia: rozwiązywanie układów równań liniowych

Ćwiczenie: Mnożyć przez odwrotność to nie zawsze to samo...

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

 
x = A \ b;

W tym celu Octave oczywiście wewnętrznie wykorzystuje funkcję LAPACKa DGESV. Ale w Octave jest także funkcja inv wyznaczająca macierz odwrotną, niektóre więc (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ć tak, jak poniżej. 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: Wybór elementu głównego jest naprawdę ważny!

Rozwiąż prosty układ równań

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned 10^{-18} x_1 + x_2 &= 1,\\ x_1 + 2x_2 &= 4, \endaligned}

czterema sposobami:

• na kartce papieru (dokładnie!)
• w komputerze, w arytmetyce podwójnej precyzji, korzystając z rozkładu LU bez wyboru elementu głównego
• jak poprzednio, ale z wyborem elementu głównego w kolumnie
• w LAPACKu lub Octave.

Porównaj uzyskane rozwiązania i wyciągnij wnioski.

Rozwiązanie

Rozwiązanie uzyskane bez wyboru elementu głównego jest bardzo mało precyzyjne $x_1$ i jest obarczone bardzo dużym błędem! Prosta analiza pokazuje, że winę za to ponosi bardzo wielki współczynnik w macierzy $L$. Nie jest tak, gdy stosujemy wybór elementu: elementy $L$ są zawsze mniejsze od 1. Dlatego inne rozwiązania numeryczne dają w przypadku naszej macierzy wynik będący bardzo dobrym przybliżeniem wyniku dokładnego. A propos: Octave korzysta z DGESV LAPACKa, a DGESV korzysta z wyboru w kolumnie...

Ć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 potró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

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: Czy Twoje programy działają naprawdę szybko?

Zrób zadanie poprzednie programując, niczym Zosia samosia, wszystko od początku do końca w czystym C. Porównaj czasy działania twojego programu i programu wywołującego po prostu procedurę biblioteczną LAPACKa DGESV, najlepiej wspartą dobrze podrasowanymi BLASami.

Wskazówka

To zadanie to wyzwanie!

Rozwiązanie

No tak, namawiamy cię do zmierzenia się z LAPACKiem i ATLASem....

Prościutki program --- po prostu tłumaczący algorytm z wykładu bezpośrednio na C --- u nas spisał się, dla macierzy wymiaru $N=1\dots 1024$, następująco:

Prosty solver w C kontra LAPACK i LAPACK z ATLASem. Zwróć uwagę na piki zwykłego kodu dla $N$ będących potęgami dwójki

Skala logarytmiczna pozwala lepiej ocenić różnice pomiędzy czasami wykonania. Prosty program w C jest około dziesięć razy wolniejszy od kodu korzystającego z LAPACKa i około 50 razy wolniejszy od kodu korzystającego z LAPACKa i ATLASa.

Ć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 algorytmów podstawienia w przód i w tył operujące od razu na wszystkich kolumnach macierzy $B$. 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 warunki:

$A=A^T\text{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), nie musimy więc 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. 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.