MN05LAB

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

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Oglądaj wskazówki i rozwiązania __SHOWALL__
Ukryj wskazówki i rozwiązania __HIDEALL__

Uwaga

Niektóre fragmenty zadań wymagają wykorzystania procedur LAPACKa; te części odłóż do momentu, gdy opanujesz następny wykład, dotyczący m.in. bibliotek algebry liniowej.

Ć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=LD{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}{\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: 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;

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ć jak poniżej. Duże macierze losowe są zazwyczaj dobrze uwarunkowane.

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ń

$\begin{array}{rl}10^{-18}{x}_1+x_2& =1,\\ x_1+2x_2& =4,\end{array}$

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$ jest obarczone bardzo dużym błędem! Prosta analiza pokazuje, że winę za to ponosi bardzo wielki współczynnik w macierzy $L$. Tak nie jest, gdy stosujemy wybór elementu: elementy $L$ są zawsze mniejsze od 1. Dlatego LAPACK i Octave dają w przypadku naszej macierzy rezultat będący bardzo dobrym przybliżeniem wyniku dokładnego.

Ć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? 
# Czy jest macierzą? Te linie pozostawiamy twojej inwencji

N = size(A,1);
for k=1:N-1
	if (A(k,k) == 0.0)
		return;
	end
	for i=k+1:N #wyznaczenie <math>k</math>-tej kolumny <math>L</math>
		A(i,k) = A(i,k)/A(k,k);
	end
	for j=k+1:N #modyfikacja podmacierzy <math>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>k</math>-tej kolumny <math>L</math>
		A(k+1:N,k) = A(k+1:N,k)/A(k,k);
#modyfikacja podmacierzy <math>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 losowych 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.

Jeśli chodzi o popełnione błędy dla macierzy Hilberta, wyjaśnienie znajdziesz w wykładzie o uwarunkowaniu układu równań liniowych.

Ć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 na przykład przy rozwiązywaniu 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>PA = LU</math>;
Utwórz macierz <math>B = [b_1,\ldots,b_k]</math>;
Rozwiąż <math>LY = B</math>;
Rozwiąż <math>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: 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.