Tekst źródłowy strony MN07


<!-- 
Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki 
wprowadzi� przykry@mimuw.edu.pl
-->
 
=Uwarunkowanie układu równań liniowych=

{{powrot |Metody numeryczne | do strony głównej 
przedmiotu <strong>Metody numeryczne</strong>}}

Zajmiemy się wrażliwością układu równań na zaburzenia danych: prawej strony i
współczynników macierzy układu. Jak zobaczymy na poniższym przykładzie, bywają
równania, które są mało podatne na zaburzenia danych (a więc: dobrze
uwarunkowane) oraz równania, które są szalenie wrażliwe na zaburzenia, a więc
źle uwarunkowane.  Jak wkrótce się przekonamy, czułość danego układu równań na
zaburzenia da się precyzyjnie scharakteryzować, a cecha ta nie tylko będzie
miała wpływ na jakość rozwiązań możliwych do uzyskania w arytmetyce skończonej
precyzji, ale także na efektywność [[MN08|metod iteracyjnych]] rozwiązywania  układów równań liniowych, w których są tysięce (lub więcej) niewiadomych.

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład: Uwarunkowanie układu dwóch równań liniowych</span> 
<div class="solution" style="margin-left,margin-right:3em;">

Rozwiązanie układu dwóch równań liniowych można przedstawić w formie graficznej:
jest to punkt przecięcia się dwóch prostych wyznaczonych przez dane
współczynniki i wyrazy prawej strony.

[[Image:MNlinearcond.png|thumb|550px|center|Rozważmy pewien nieosobliwy układ dwóch równań
liniowych. Ma on dokładnie jedno rozwiązanie, na rysunku oznaczone kolorem czerwonym. Co się stanie, gdy trochę zaburzymy prawą stronę takiego układu?]]

[[Image:MNlinearcond1.png|thumb|550px|center|Wykresy zaburzonych prostych mogą zająć jedną z
zaznaczonych łososiowym kolorem pozycji.]]

[[Image:MNlinearcond2.png|thumb|550px|center|Obszar, gdzie mogą znaleźć się rozwiązania zaburzonego
układu, zaznaczyliśmy na czerwono. Jest on, kolokwialnie rzecz ujmując, z
grubsza tak
wielki jak wielkie były zaburzenia, co zgodne jest z typową intuicją "człowieka
z zewnątrz".]]

[[Image:MNlinearcond3.png|thumb|550px|center|Jednak bywają
równania, wrażliwe jak mimoza na nawet delikatne zaburzenia danych. Takie
równanie właśnie widzimy na rysunku: jego cechą szczególną jest to, że tym razem
proste, choć wciąż przecinają się dokładnie w jednym punkcie, są <strong>prawie</strong>
równoległe.]]

[[Image:MNlinearcond4.png|thumb|550px|center|Bierzemy zaburzenia takie same, jak poprzednio. Wykresy zaburzonych prostych mogą zająć jedną z zaznaczonych łososiowym kolorem pozycji.]]

[[Image:MNlinearcond5.png|thumb|550px|center|Tym razem obszar niepewności, gdzie mogą być
rozwiązania naszego zaburzonego układu, jest <strong>gigantyczny</strong>!]]

</div></div>

A więc równania liniowe mogą, choć nie muszą, być bardzo podatne na zaburzenia
danych. Gdy zamiast prawej strony, zaburzymy wyrazy macierzy układu, może nawet
okazać się, że dostaniemy układ równań sprzecznych (czy możesz podać przykład?)

Aby przedstawić ogólną teorię zaburzeń dla układów równań liniowych, musimy mieć
narzędzia do pomiaru błędu rozwiązań, a także zaburzeń danych zadania: czyli
macierzy i wektora prawej strony. Temu będą służyć normy.

==Normy wektorowe i macierzowe==

Aby badać odległość między rozwiązaniem dokładnym układu równań a jego
wartością przybliżoną uzyskaną np. algorytmem eliminacji Gaussa, będziemy
posługiwać się normami wektorów 
<math>x = (x_j)_{j=1}^n\in R^n</math> 
i macierzy <math>A = (a_{i,j})_{i,j=1}^n \in R^{n\times n}</math>. 
Najczęściej używanymi normami wektorowymi będą 
<strong>normy <math>p</math>-te</strong>,

<center><math>\| x\|\,=\,\| x\|_p\,=\,
   \left(\sum_{j=1}^n |x_j|^p\right)^{1/p},
   \qquad 1\le p< +\infty</math>,</center>

oraz 

<center><math>\| x\|_\infty\,=\,\lim_{p\to +\infty}\| x\|_p\,=\,
    \max_{1\le j\le n}|x_j|. 
</math></center>

W szczególności, norma <math>||\cdot||_2</math> jest dobrze nam znaną normą euklidesową wektora.

[[Image:MNball1.png|thumb|300px|center|Kula jednostkowa w normie <math>||\cdot||_1</math> w <math>R^2</math>]]
[[Image:MNball2.png|thumb|300px|center|Kula jednostkowa w normie <math>||\cdot||_2</math> w <math>R^2</math>]]
[[Image:MNballinf.png|thumb|300px|center|Kula jednostkowa w normie <math>||\cdot||_\infty</math> w <math>R^2</math>]]

Normą macierzową jest norma Frobeniusa 

<center><math>\|A\|_F\,=\,\sqrt{\sum_{i,j=1}^n |a_{i,j}|^2}</math>,</center>

a także <strong>normy indukowane</strong> przez normy wektorowe (np. przez 
normy <math>p</math>-te) 

<center><math>\|A\|\,=\,\sup_{x\ne 0}\frac{\|A x\|}{\| x\|}\,=\,
      \sup_{\|x\|=1}\|A x\|</math></center>

Jeśli norma macierzowa jest indukowana przez normę wektorową, 
to dla dowolnego wektora mamy 

<center><math>\|A x\|\,\le\,\|A\|\| x\|.  
</math></center>

Przypomnijmy, że w przestrzeniach liniowych skończenie wymiarowych 
(a więc także w <math>R^n</math> i w przestrzeni macierzy wymiaru <math>n\times n</math>) 
każde dwie normy są równoważne. To znaczy, że jeśli mamy dwie 
normy <math>\|\cdot\|</math> i <math>\|\cdot\|'</math> w przestrzeni skończenie wymiarowej 
<math>X</math>, to istnieją stałe <math>0<K_1\le K_2<\infty</math> takie, że 

<center><math>K_1\,\|x\|\,\le\,\|x\|'\,\le\,K_2\,\|x\|,\qquad\forall x\in X</math></center>

W szczególności dla <math>x\in R^n</math> mamy

<center><math>\begin{align} \| x\|_\infty &\le & \| x\|_1\,\le\,
        n\,\| x\|_\infty, \\
  \| x\|_\infty &\le & \| x\|_2\,\le\,
         \sqrt{n}\,\| x\|_\infty,\\
  \frac 1{\sqrt n}\,\| x\|_1 &\le & \| x\|_2\,\le\,
        \| x\|_1, 
\end{align}</math></center>

a dla <math>A=(a_{i,j})_{i,j=1}^n</math> mamy 

<center><math>
  \|A\|_2\, \le\, \|\,|A|\,\|_2 \,\le\, \|A\|_E 
     \,\le\, \sqrt n\, \|A\|_2</math>,</center>

gdzie <math>|A|=(|a_{i,j}|)_{i,j=1}^n</math>.  

Dla macierzy 
<math>A=(a_{i,j})_{i.j=1}^n\in R^{n\times n}</math> mamy 

<center><math>\|A\|_\infty \,=\, \max_{1\le i\le n}\sum_{j=1}^n |a_{i,j}| 
</math></center>

oraz 

<center><math>\|A\|_1 \,=\, \|A^T\|_\infty \,=\,
         \max_{1\le j\le n}\sum_{i=1}^n |a_{i,j}|</math></center>

Dowód tego faktu zostawiamy jako [[MN07LAB|ćwiczenie]].

==Uwarunkowanie układu równań liniowych==

Wyprowadzimy teraz wynik świadczący o tym, jak zaburzenie względne danych
przenosi się na błąd względny wyniku rozwiązania <math>x^*</math> układu równań liniowych <math>Ax=b</math>.

{{twierdzenie|O uwarunkowaniu układu równań|O uwarunkowaniu układu równań|

Niech <math>E</math> i <math>e</math> będą zaburzeniami 
odpowiednio macierzy <math>A</math> i wektora <math>b</math> na poziomie względnym <math>\epsilon</math>, tzn.

<center><math>\frac{\|E\|}{\|A\|} \,\le\,\epsilon \qquad  \mbox{i}  \qquad \frac{\|e\|}{\|b\|} \,\le\,\epsilon</math></center>

Jeśli 

<center><math>\epsilon\cdot \mbox{cond} (A)\,<\,1</math>,</center>

to układ zaburzony <math>(A+E) x=( b+ e)</math> ma jednoznaczne 
rozwiązanie <math>z^*</math> spełniające 

<center><math>\frac{\| z^*- x^*\|}{\| x^*\|} \; \le \; \frac{2\, \mbox{cond} (A)}{1-\epsilon \mbox{cond} (A)} \epsilon</math>,</center>

gdzie definiujemy <strong>współczynnik uwarunkowania układu</strong>

<center><math>\mbox{cond} (A) = ||A||\cdot ||A^{-1}||</math></center>

}}

Zauważmy najpierw, że zachodzi 

{{lemat|Neumanna o otwartości zbioru macierzy odwracalnych|Neumanna o otwartości zbioru macierzy odwracalnych|

Jeśli <math>F</math> jest macierzą 
taką, że <math>\|F\|<1</math>, to macierz <math>(I-F)</math> jest nieosobliwa oraz 

<center><math>\| (I-F)^{-1} \|\,\le\,\frac{1}{1-\|F\|}</math></center>

}}

{{dowod|||
Rzeczywiście, gdyby <math>(I-F)</math> była osobliwa, to istniałby niezerowy 
wektor <math>x</math> taki, że <math>(I-F) x=0</math>, co implikuje 
<math>\|F x\|/\| x\|=1</math> i w konsekwencji <math>\|F\|\ge 1</math>. Aby 
pokazać oszacowanie normy macierzy <math>(I-F)^{-1}</math> zauważmy, że 

<center><math>\begin{align} 1 &= \|I\|\,=\,\|(I-F)(I-F)^{-1}\| \\ &\ge & 
    \|(I-F)^{-1}\|\,-\,\|F\|\,\|(I-F)^{-1}\| \\
    &= (1-\|F\|)\,\|(I-F)^{-1}\|,
\end{align}</math></center>

skąd już wynika dowodzona nierówność. 
}}

{{dowod|twierdzenia o uwarunkowaniu|twierdzenia o uwarunkowaniu|

Po podstawieniu <math>F=-A^{-1}E</math> mamy teraz 

<center><math>\|F\|\,\le\,\|A^{-1}\|\,\|E\|\,\le\,\epsilon\|A\|\,\|A^{-1}\|\,<\,1</math>,</center>

co wobec równości <math>A+E=A(I+A^{-1}E)</math> daje, że macierz <math>(A+E)</math> 
jest nieosobliwa i układ zaburzony ma jednoznaczne rozwiązanie 
<math>z^*</math>. Przedstawmy to rozwiązanie w postaci 
<math>z^*= x^*+( z^*- x^*)</math>. Rozpisując układ 
zaburzony i wykorzystując równość <math>A x^*= b</math> otrzymujemy, 
że <math>(A+E)( z^*- x^*)= e\,-\,E x^*</math>, czyli  

<center><math>z^*- x^* \,=\, (I+A^{-1}E)^{-1}A^{-1}( e-E x^*)</math>,</center>

a stąd

<center><math>\begin{align} \| z^*- x^*\| &\le & \|(I+A^{-1}E)^{-1}\|\,\|A^{-1}\| \,
      (\| e\|+\|E\|\,\| x^*\| \\
   &\le & \frac{\|A^{-1}\|}{1-\epsilon\|A\|\,\|A^{-1}\|}
      \epsilon\left(\| b\|+\|A\|\,\| x^*\|\right) \\
   &\le & \frac{\|A\|\,\|A^{-1}\|}{1-\epsilon\|A\|\,\|A^{-1}\|}
      2\epsilon\cdot\| x^*\|, 
\end{align}</math></center>

co kończy dowód. 
}}

Gdy więc np. <math>\epsilon  \mbox{cond} (A) \leq \frac{1}{2}</math>, [[#O uwarunkowaniu układu równań|oszacowanie błędu rozwiązania układu zaburzonego]] możemy
zastąpić czytelniejszym (choć mniej precyzyjnym)

<center><math>\frac{\| z^*- x^*\|}{\| x^*\|} \leq 4 \,  \mbox{cond} (A) \, \epsilon</math></center>

Octave i MATLAB mają wbudowane funkcje wyznaczające normy wektorów i macierzy
 <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>N = 3;
x = [1:N]'
A = pascal(N)
norm(A,1)
norm(x,2)
norm(A,Inf)
</pre></div> 
 
a także funkcje wyznaczające uwarunkowanie macierzy, przy czym Octave liczy
tylko uwarunkowanie w normie <math>||\cdot||_2</math>:

 <div style="margin: 1em; padding:1em; color: #006; background-color:#fcfcfc;"><pre>cond(A)
</pre></div> 
 
W [http://www.netlib.org/lapack  LAPACKu] służy do tego funkcja <code style="color: #903">DGECON</code>. Zadanie wyznaczania uwarunkowania macierzy jest zadaniem bardzo intensywnym numerycznie. Problem, czy da się je wyznaczyć z dobrą dokładnością kosztem niższym niż wyznaczenie macierzy odwrotnej i jej normy, jest wciąż otwarty.

W praktyce obliczeniowej trafiają się zarówno układy dobrze uwarunkowane, jak i
macierze, których uwarunkowanie może być patologicznie duże (np. takie macierze
są chlebem powszednim osób rozwiązujących równania różniczkowe).

<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład: Macierz Hilberta</span> 
<div class="solution" style="margin-left,margin-right:3em;">

Przykładem macierzy o uwarunkowaniu wyjątkowo szybko rosnącym z wymiarem jest m.in. <strong>macierz Hilberta</strong> <math>H_N = (h_{ij})_{i,j=1}^N</math>, gdzie
<center><math>
h_{ij} = \frac{1}{i+j-1}</math></center>

[[Image:MNhilbertmatrix.png|thumb|550px|center|Macierz Hilberta wymiaru 25. Kolor odpowiada rzędowi wielkości elementu macierzy, dokładniej, <math>\log(h_{ij})</math>. Jak widzisz, większość elementów tej macierzy jest równa prawie zero, a więc w konsekwencji kolumny macierzy są prawie liniowo zależne. ]]

Taką macierz możemy wygenerować w Octave komendą <code style="color: #006">hilb(N)</code>. Jest to bardzo specyficzna macierz, co m.in. przejawia się tym, że uwarunkowanie macierzy Hilberta rośnie eksponencjalnie z <math>N</math>,  <math>\mbox{cond} (H_N) \approx  O(e^{3.5N})</math>:

<div style="font-family: monospace; white-space: pre; border-style: dashed; border-width: thin; border-color: black; margin: 1em; padding:1em; color: #444; background-color:#fdfdfd;"><realnowiki><realnowiki><nowiki>octave:2> cond(hilb(5))
ans =  4.7661e+05
octave:3> cond(hilb(10))
ans =  1.6025e+13
octave:4> cond(hilb(15))
ans =  3.7689e+17
octave:5> cond(hilb(20))
ans =  7.1209e+19
</nowiki></realnowiki></realnowiki></div>
 
</div></div>

==Numeryczna poprawność eliminacji Gaussa==

Przedstawimy bez dowodu klasyczne twierdzenie o "praktycznej numerycznej
poprawności" eliminacji Gaussa z wyborem w kolumnie.

{{twierdzenie|Wilkinsona|Wilkinsona|

Algorytm eliminacji Gaussa z wyborem elementu głównego w kolumnie, zrealizowany
w arytmetyce <math>fl_\nu</math>,

wyznacza <math>\widetilde{x}</math>  taki, że  <math>\widetilde{x}</math> jest <strong>dokładnym</strong> rozwiązaniem zadania zaburzonego

<center><math>\widetilde{A}\widetilde{x} = b</math>,</center>

przy czym 

<center><math>\frac{||A-\widetilde{A}||_\infty}{||A||_\infty} \leq K \, N^3 \, \rho_N \, \nu</math>,</center>

dla pewnej niedużej stałej <math>K = O(1)</math>. <strong>Wskaźnik wzrostu</strong> <math>\rho_N</math> definiujemy tutaj jako 
<center><math>\rho_N = \frac{\max_{i,j}|\widetilde{u}_{ij}|}{\max_{i,j} |a_{ij}|}</math>,</center>


gdzie <math>\widetilde{L}</math> i <math>\widetilde{U}</math> są numerycznie wyznaczonymi czynnikami rozkładu PA<math>=</math>LU.
}}

Jak widzimy, kluczowe dla numerycznej poprawności jest oszacowanie ''wskaźnika wzrostu'' <math>\rho_N</math>. Okazuje się, co wiedział już Wilkinson, że
* w ogólnym przypadku, zachodzi oszacowanie <math>\rho_N \leq 2^{N-1}</math>, które jest osiągane dla macierzy 

<center><math>W = \begin{pmatrix} 
1 & & & 1 \\
-1 & \ddots & & \vdots\\
\vdots &  & \ddots  & \vdots\\
-1 & \cdots & -1 & 1\\
\end{pmatrix} ;
</math></center>
* dla macierzy trójdiagonalnych lub diagonalnie dominujących, lub dla macierzy symetrycznych dodatnio określonych,  <math>\rho_N \leq 2</math>;
* w średnim przypadku, obserwuje się <math>\rho_N \leq N^{2/3}</math>, to znaczy macierze spotykane w praktyce obliczeniowej mają mały wskaźnik wzrostu.
 
Konkluzja jest więc taka, że <strong>algorytm eliminacji Gaussa z wyborem w kolumnie
jest ''praktycznie'' numerycznie poprawny</strong>. Z drugiej strony, dla bardzo dużych
<math>N</math> i niezbyt dobrze uwarunkowanych macierzy, może okazać się, że arytmetyka
pojedynczej precyzji może okazać się niewystarczająca dla uzyskania godnego
wyniku.

Algorytm eliminacji Gaussa z pełnym wyborem elementu głównego jest
numerycznie poprawny, ze wskaźnikiem wzrostu <math>\rho_N \leq \sqrt{n\cdot 2 \cdot
3^{1/2}\cdot 4^{1/3}\cdots N^{1/(N-1)}}</math>, a w praktyce grubo poniżej <math>\sqrt{N}</math>.

==Literatura==

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 4.4  i, nieobowiązkowo, 4.5</b> w
* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Bardzo dużo informacji na temat omawianych zagadnień znajduje się w monografiach 
* <span style="font-variant:small-caps">A.Kiełbasiński, H. Schwetlick</span>, <cite>Numeryczna algebra liniowa</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa, 1992,
* <span style="font-variant:small-caps">N. Higham</span>, <cite>Accuracy and Stability of Numerical Algorithms</cite>, SIAM, 2002.