MN04: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
Przykry (dyskusja | edycje)
mNie podano opisu zmian
m Zastępowanie tekstu – „,↵</math>” na „</math>,”
 
(Nie pokazano 19 wersji utworzonych przez 3 użytkowników)
Linia 1: Linia 1:


=Układy równań liniowych=
<!--
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
-->
=Własności zadania obliczeniowego i algorytmu numerycznego=


Czasem zadania obliczeniowe wymagają wykonania naprawdę wielkiej liczby obliczeń
{{powrot |Metody numeryczne | do strony głównej
zmiennoprzecinkowych, przy czym wiele znanych, <strong>matematycznie równoważnych</strong> metod
przedmiotu <strong>Metody numeryczne</strong>}}
rozwiązywania takich zadań, ma <strong>diametralnie różne własności numeryczne</strong>.
Bardzo ważną klasą takich zadań jest rozwiązywanie układów równań liniowych


<center><math>\displaystyle
==Uwarunkowanie zadania obliczeniowego==
Ax = b,
</math></center>


gdzie <math>\displaystyle A</math> jest nieosobliwą macierzą <math>\displaystyle N\times N</math>, a dany wektor prawej strony <math>\displaystyle b\in R^N</math>.  
Jak zobaczyliśmy w poprzednich przykładach, dane, jakimi dysponujemy wykonując
zadanie obliczeniowe, są z natury rzeczy wartościami zaburzonymi. Źródłem tych
zaburzeń są:
* błąd reprezentacji danych w arytmetyce zmiennoprzecinkowej (np. 0.1 nie jest równe dokładnie <math>1/10</math>)
* błędy w parametrach będących rezultatem poprzednich obliczeń (np. chcemy rozwiązać równanie <math>f(x) = a</math>, ale <math>a</math> jest rezultatem innej symulacji), a także
* błędy pomiarowe w zadaniach praktycznych (np. chcemy policzyć numeryczną prognozę pogody, ale temperaturę wyjściową znamy tylko z dokładnością do 0.1 stopnia, co gorsza --- z niektórych stacji w ogóle nie mamy danych)
Okazuje się, że powszechna intuicja, że małe zaburzenia danych powinny dawać
małe zaburzenia wyniku, nie znajduje potwierdzenia nawet w bardzo prostych
przypadkach. Z drugiej strony, umiejętność oceny jakościowego <strong>wpływu
zaburzenia danych na wynik</strong> jest kapitalna w świecie obliczeń numerycznych w
ogólności, a w szególności --- inżynierskich.


W
Wprowadza się pojęcie <strong>uwarunkowania</strong> zadania, to znaczy jego podatności na
praktyce spotyka się zadania z <math>\displaystyle N = 2, 3, \ldots 1000</math>. Zdarzają się także czaem
zaburzenia danych. Dla przejrzystości przypuśćmy, że nasze zadanie obliczeniowe
specjalne macierze o wymiarach nawet rzędu <math>\displaystyle 10^8</math>!
polega na wyznaczeniu <math>f(x)</math> dla danego <math>x</math>.


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ń
Jak bardzo będzie odległe
liniowych, takie jak:
<math>f(\widetilde{x})</math>, gdy <math>\widetilde{x}\approx x</math>? Rozważa się dwa przypadki:
* metoda wyznacznikowa (wzory Cramera)
* <strong>uwarunkowanie względne</strong>: jak względne zaburzenie danych wpływa na błąd względny wyniku: <center><math>\frac{||f(x) - f(\widetilde{x})||}{||f(x)||} \leq  \mbox{cond} _{rel}(f,x) \cdot \frac{||x - \widetilde{x}||}{||x||}</math></center>
* obliczenie macierzy <math>\displaystyle A^{-1}</math> i następnie <math>\displaystyle x = A^{-1}b</math>
Najmniejszy mnożnik <math>\mbox{cond} _{rel}(f,x)</math> spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia <math>f(x)</math> dla danego <math>x</math>.
* <strong>uwarunkowanie bezwzględne</strong>: jak bezwzględne zaburzenie danych wpływa na błąd bezwzględny wyniku: <center><math>||f(x) - f(\widetilde{x})|| \leq  \mbox{cond} _{abs}(f,x) \cdot ||x - \widetilde{x}||</math></center>
Najmniejszy mnożnik <math>\mbox{cond} _{abs}(f,x)</math>  spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia <math>f(x)</math> dla danego <math>x</math>.
   
   
<strong>nie nadają się</strong> do numerycznego rozwiązywania takich zadań.
Powiemy, że zadanie <math>f(x)</math> jest
* <strong>dobrze uwarunkowane</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) \approx 1</math>,
* <strong>źle uwarunkowane</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) \gg 1</math>,
* <strong>źle postawione</strong> w punkcie <math>x</math>, gdy <math>\mbox{cond} (f,x) = +\infty</math>.
Teraz już rozumiemy, że redukcja cyfr przy odejmowaniu jest tylko
odzwierciedleniem faktu, że zadanie odejmowania dwóch bliskich liczb jest po
prostu zadaniem źle uwarunkowanym!
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład: Uwarunkowanie zadania obliczenia sumy</span>
<div class="solution" style="margin-left,margin-right:3em;">
 
Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia <math>s(x,y) = x + y</math> ma
<center><math>
\mbox{cond} _{abs}(s, (a,b)) = 1, \qquad  \mbox{cond} _{rel}(s, (a,b)) = \frac{|a|+|b|}{|a+b|}
</math></center>
 
Tak więc, gdy <math>a\approx -b</math>, to <math>\mbox{cond} _{rel}(s, (a,b)) \approx +\infty</math> i zadanie
jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może
skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego,
najczęściej rzeczywiście tak będzie...
</div></div>
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład</span>
<div class="solution" style="margin-left,margin-right:3em;">


O tym, jak <strong>skutecznie</strong> rozwiązywać takie zadania, jakie jest ich
Dla różniczkowalnej funkcji skalarnej <math>f : R \rightarrow R</math> mamy
uwarunkowanie --- o tym traktują następne dwa wykłady. Okaże się, że jedną z
<center><math>
dobrych metod jest metoda eliminacji Gaussa.
|f(x) - f(\widetilde{x})| \approx |f'(x) | | x - \widetilde{x} |
</math></center>


==Proste układy równań==
i w konsekwencji dla zadania obliczenia <math>f(x)</math> dla danego <math>x</math> mamy, przy
założeniu małych zaburzeń,
<center><math>
\mbox{cond} _{abs}( f, x) = |f'(x)|, \qquad  \mbox{cond} _{rel}( f, x) =
\frac{|f'(x)|\cdot|x|}{|f(x)|}</math></center>


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


<blockquote style="background-color:#fefeee">
Jest zrozumiałe, że złe uwarunkowanie jest szkodliwe w praktyce numerycznej: zaburzenia danych są nieuniknione, a źle uwarunkowane zadanie tylko je wzmocni na wyjściu. Jednak za jakiś czas zobaczysz wyjątkową sytuację, gdy [[MN13#Odwrotna metoda potęgowa|złe uwarunkowanie pewnego podzadania nie tylko nie przeszkadza, ale wręcz <strong>pomaga</strong>]] szybciej rozwiązać zadanie główne!
Trudne zadania rozwiązujemy sprowadzając je do sekwencji łatwych zadań
</blockquote>  


w dalszej kolejności pokażemy, jak dowolny układ równań sprowadzić do sekwencji
==Rozkład algorytmu względem informacji==
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ą====
<strong>Algorytm</strong> to dokładnie określona i dozwolona w danym modelu
obliczeniowym sekwencja akcji, pozwalająca na rozwiązanie danego
zadania (w sposób dokładny lub przybliżony).


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


<center><math>\displaystyle L = \begin{pmatrix}  
<center><math>{\bf ALG}:\,F\longrightarrow G,  
1 &  &  &        &  &  \\
* & 1 &  &        &  &  \\
* & * & 1 &        &  &  \\
* & * & * & 1 &  &        \\
\vdots  & \vdots  & \vdots  & \ddots  & \ddots &  \\
*  &  *  & * &  \cdots  &  *    & 1
\end{pmatrix} ,  
\qquad
U = \begin{pmatrix}
* & * & * & *      & \cdots & * \\
  & * & * & *      & \cdots & * \\
  &  & * & *      & \cdots & * \\
  &  &        & * & \ddots &  \vdots \\
  &  &  &        & \ddots & * \\
  &  &  &        &        & * \end{pmatrix}
</math></center>
</math></center>


Układ z nieosobliwą macierzą trójkątną górną
taki że <math>{\bf ALG}(f)</math> jest wynikiem działania algorytmu
w arytmetyce idealnej dla danej <math>f</math>.


<center><math>\displaystyle
Zauważmy, że wynik <math>{\bf ALG}(f)</math> działania algorytmu nie
  U\, x\;=\; c,
zależy bezpośrednio od <math>f</math>, ale raczej od <strong>informacji</strong>
o <math>f</math> (uzyskanej dzięki poleceniu <math>{\cal IN}</math>). Informacja
ta może być <strong>pełna</strong> albo tylko <strong>częściowa</strong>.
Informacja jest pełna gdy, np.
<math>f=(f_1,\ldots,f_n)\in R^n</math> i wczytamy wszystkie
współrzędne <math>f_i</math>. Informacja może być częściowa, gdy
<math>f</math> jest funkcją. Wtedy wiele danych może posiadać tę
samą informację, co łatwo zaobserwować na przykładzie
zadania całkowania.
 
Niech <math>N:F\to \cup_{n=0}^\infty R^n</math> będzie
<strong>operatorem informacji</strong>, tzn.
 
<center><math>N(f)\,=\,(y_1,y_2,\ldots,y_n)
</math></center>
</math></center>


<math>\displaystyle U=(u_{i,j})</math>, <math>\displaystyle  c=(c_j)</math>, można rozwiązać stosując algorytm:
jest informacją o <math>f</math> zebraną przy idealnej realizacji
algorytmu. Zauważmy, że informacja jest pełna, gdy <math>N</math> jest
przekształceniem różnowartościowym, tzn. jeśli
<math>f_1\ne f_2</math> implikuje <math>N(f_1)\ne N(f_2)</math>.
W przeciwnym przypadku mamy do czynienia z informacją
częściową.


{{algorytm|Podstawienie w tył||
Każdy algorytm <math>{\bf ALG}</math> może być przedstawiony jako złożenie
<pre>
operatora informacji i pewnego operatora
<math>\varphi:N(F)\to G</math> zdefiniowanego równością


<math>\displaystyle x_N^*\, =\, c_N / u_{N,N}</math>;
<center><math>\varphi\left(N(f)\right)\,=\,{\bf ALG}(f).
for (i = N-1; i >= 1; i--)
</math></center>
<math>\displaystyle x_i^*\,:=\,\left( c_i\,-\, \sum_{j=i+1}^N
 
u_{i,j}x_j^*\right)/u_{i,i}</math>;
Zauważmy, że w przypadku informacji częściowej zwykle nie
</pre>}}
istnieje algorytm dający dokładne rozwiązanie zadania dla
każdej danej <math>f\in F</math>, ponieważ dla danych o tej samej
informacji mogą istnieć różne rozwiązania.


(Algorytm ten jest wykonalny, ponieważ nieosobliwość macierzy
==Problem wyboru algorytmu==
implikuje, że <math>\displaystyle u_{i,i}\ne 0</math>, <math>\displaystyle \forall i</math>.) Podobnie, układ
<math>\displaystyle L x= c</math> rozwiązujemy algorytmem:


{{algorytm|Podstawienie w przód||
Wybór algorytmu jest najistotniejszą częścią całego procesu
<pre>
numerycznego rozwiązywania zadania. Kierujemy się przy tym przede
wszystkim następującymi kryteriami:
* dokładnością algorytmu,
* złożonością algorytmu,
* własnościami numerycznymi algorytmu.
Przez dokładność algorytmu rozumiemy różnicę między
rozwiązaniem dokładnym <math>S(f)</math> a rozwiązaniem
<math>{\bf ALG}(f)</math> dawanym przez algorytm w arytmetyce idealnej.
Jeśli <math>{\bf ALG}(f) = S(f)</math>,
<math>\forall f \in F</math>,
to algorytm nazywamy <strong>dokładnym</strong>.


<math>\displaystyle x_1 = c_1</math>;
Mówiąc o złożoności, mamy na myśli złożoność pamięciową
for (i=2; i <= N; i++)
(zwykle jest to liczba stałych i zmiennych używanych przez
<math>\displaystyle x_i = c_i\,-\,\sum_{j=1}^{i-1} l_{i,j} x_j^*</math>;
algorytm), jak również złożoność obliczeniową.
</pre>}}
Na złożoność obliczeniową algorytmu dla danej <math>f</math> składa
się koszt uzyskania infomacji <math>y=N(f)</math> (zwykle jest on
proporcjonalny do liczby wywołań polecenia <math>{\cal IN}</math>), oraz
koszt <strong>kombinatoryczny</strong> przetworzenia tej informacji, aż do
uzyskania wyniku <math>\varphi(y)</math>. Koszt kombinatoryczny zwykle
mierzymy liczbą operacji arytmetycznych wykonywanych przez
algorytm.


Oba algorytmy wymagają rzędu <math>\displaystyle N^2/2</math> mnożeń lub dzieleń i
Przez własności numeryczne algorytmu rozumiemy jego
<math>\displaystyle N^2/2</math> dodawań lub odejmowań, a więc łącznie <math>\displaystyle O(N^2)</math>
własności przy realizacji w arytmetyce <math>fl_\nu</math>. Temu
działań arytmetycznych.  
ważnemu tematowi poświęcimy teraz osobny paragraf.  


====Układy z macierzą ortogonalną====
==Numeryczna poprawność algorytmu==


Równie tanio można rozwiązać układ równań
Pożądane jest, aby algorytm dawał "dobry" wynik zarówno
w arytmetyce idealnej, jak i w arytmetyce <math>fl_\nu</math>. Niestety,
jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm
jest dokładny, to w wyniku jego realizacji w <math>fl_\nu</math> możemy
otrzymać wynik <math>fl_\nu({\bf ALG}(f))</math> daleko odbiegający od
<math>S(f)</math>. W szczególności, prawie zawsze mamy


<center><math>\displaystyle
<center><math>S(f)\,\ne\,fl_\nu\left({\bf ALG}(f)\right).
Q x =  b,
</math></center>
</math></center>


gdy <math>\displaystyle Q</math> jest macierzą ortogonalną, to znaczy <math>\displaystyle Q^TQ = I</math>. Rzeczywiście, z
Zauważmy również, że o ile z reguły znamy dokładne zachowanie
ortogonalności mamy natychmiast, że
się algorytmu w arytmetyce idealnej dla danej informacji, to nie
można tego samego powiedzieć o jego zachowaniu się w arytmetyce
<math>fl_\nu</math>. W związku z tym powstaje pytanie, jak kontrolować błąd 
algorytmów wynikający z błędów zaokrągleń i jakie algorytmy
uznamy za te o najwyższej jakości numerycznej.
 
Istnienie błędów reprezentacji liczb rzeczywistych powoduje,
że informacja <math>y=N(f)</math> o danej <math>f</math> nie jest w
ogólności reprezentowana dokładnie. Znaczy to, że zamiast na
informacji dokładnej, dowolny algorytm będzie operować na
informacji <strong>nieco zaburzonej</strong> <math>y_\nu</math>, tzn. zaburzonej na
poziomie błędu reprezentacji. Tak samo wynik dawany przez algorytm
będzie w ogólności zaburzony na poziomie błędu reprezentacji.
W najlepszym więc wypadku wynikiem działania algorytmu w <math>fl_\nu</math>
będzie <math>(\varphi(y_\nu))_\nu</math> zamiast <math>\varphi(y)</math>. Algorytmy
dające tego rodzaju wyniki uznamy za posiadające najlepsze
własności numeryczne w arytmetyce <math>fl_\nu</math> i nazwiemy numerycznie
poprawnymi.
 
Powiemy, że ciąg rzeczywisty
<math>a_\nu=(a_{\nu,1},\ldots,a_{\nu,n})</math>
(a właściwie rodzina ciągów <math>\{a_\nu\}_\nu</math>) jest
<strong>nieco zaburzonym</strong> ciągiem <math>a=(a_1,\ldots,a_n)</math>, jeśli
istnieje stała <math>K</math> taka, że dla wszystkich dostatecznie
małych <math>\nu</math> zachodzi


<center><math>\displaystyle
<center><math>
x = Q^T  b
  |a_{\nu,j} - a_j|\,\le\,K\,\nu\,|a_j|,\qquad 1\le j\le n,
</math></center>
</math></center>


i w konsekwencji <math>\displaystyle x</math> można wyznaczyć kosztem takim jak koszt mnożenia macierzy
albo ogólniej
przez wektor, czyli <math>\displaystyle O(N^2)</math> operacji.  
 
<center><math>
  \|a_\nu - a\|\,\le\,K\,\nu\,\|a\|</math>,</center>
 
gdzie <math>\|\cdot\|</math> jest pewną normą w <math>R^n</math>. W pierwszym
przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim
o zaburzeniu w normie <math>\|\cdot\|</math>.
 
Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają
za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas
 
<center><math>\|a_\nu - a\|_\infty \,=\, \max_{1\le j\le n} |a_{\nu,j} - a_j|
  \,\le\,K\,\nu\,\max_{1\le j\le n} |a_j|\,=\,K\,\nu\,\|a\|_\infty</math>,</center>
 
i korzystając z faktu, że w przestrzeni skończenie wymiarowej
wszystkie normy są równoważne, otrzymujemy dla pewnych stałych
<math>K_1</math> i <math>K_2</math>
 
<center><math>\|a_\nu - a\|\,\le\,K_1\|a_\nu-a\|_\infty\,\le\,
    K_1 K\,\nu\,\|a\|_\infty\,\le\,K_2 K_1 K\,\nu\,\|a\|</math>,</center>
 
czyli nierówność dla zaburzenia w normie, ze stałą <math>K = K_2 K_1 K</math>.
 
{{definicja|Algorytm numerycznie poprawny|Algorytm numerycznie poprawny|
 
Algorytm <math>{\bf ALG}</math> rozwiązywania zadania
nazywamy <strong>numerycznie poprawnym</strong> w zbiorze danych
<math>F_0\subset F</math>, jeśli dla każdej danej <math>f\in F_0</math>
wynik <math>fl_\nu({\bf ALG}(f))</math> działania algorytmu w arytmetyce
<math>fl_\nu</math> można zinterpretować jako nieco zaburzony wynik
algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji
<math>y_\nu=(N(f))_\nu\in N(F)</math> o <math>f</math>, przy czym
poziom zaburzeń nie zależy od <math>f</math>.  


Podobnie, gdy <math>\displaystyle Q\in C^{N\times N}</math> jest unitarna, to znaczy <math>\displaystyle Q^HQ = I</math>,
Formalnie znaczy to, że istnieją stałe <math>K_1</math>, <math>K_2</math>  oraz
rozwiązaniem układu równań jest
<math>\nu_0>0</math> takie, że spełniony jest następujący warunek.
Dla dowolnej <math>\nu\le\nu_0</math> oraz informacji <math>y\in N(F_0)</math>
można dobrać <math>y_\nu\in N(F)</math> oraz
<math>\left(\varphi(y_\nu)\right)_\nu</math> takie, że


<center><math>\displaystyle
<center><math>\|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|,
x = Q^H  b.
</math></center>
</math></center>


==Metoda eliminacji Gaussa==
<center><math>\|\left(\varphi(y_\nu)\right)_\nu - \varphi(y_\nu)\|\,\le\,
    K_2\,\nu\,\|\varphi(y_\nu)\|</math>,</center>


[[grafika:Gauss.jpg|thumb|right||Carl Friedrich Gauss<br>  [[Biografia Gauss|Zobacz biografię]]]]
oraz


W przypadku dowolnych macierzy, bardzo dobrym algorytmem numerycznym
<center><math>fl_\nu\left({\bf ALG}(f)\right)\,=\,
rozwiązywania układu równań
      fl_\nu\left(\varphi(N(f))\right)\,=\,
      \left(\varphi(y_\nu)\right)_\nu.
</math></center>


<center><math>\displaystyle Ax=b</math></center>
}}


okazuje się popularna
[[Image:MNcondition7.png|thumb|550px|center|Numerycznie poprawny algorytm daje w arytmetyce <math>fl_\nu</math> wynik <math>ALG(N(x))</math>, który daje
eliminacja Gaussa. Jednak z powodów, które będą dla nas później jasne, algorytm
się zinterpretować jako mało zaburzony wynik <math>f(y)</math> zadania na mało zaburzonych
ten wyrazimy w języku tzw. <strong>rozkładu LU</strong> macierzy, to znaczy,
danych <math>x</math>.]]
sprowadzającego zadanie do znalezienia
 
macierzy trójkątnej dolnej <math>\displaystyle L</math> (z jedynkami na diagonali) oraz trójkątnej górnej
Zauważmy,że jeśli <math>f\in R^n</math>,
<math>\displaystyle U</math> takich, że
<math>N(f)=(f_1,\ldots,f_n)</math>, oraz algorytm jest
<center><math>\displaystyle
dokładny, <math>{\bf ALG}\equiv\varphi\equiv S</math>, to numeryczną
A = LU,
poprawność algorytmu można równoważnie zapisać jako
 
<center><math>fl_\nu\left({\bf ALG}(f)\right)\,=\,
  \left(S(f_\nu)\right)_\nu.
</math></center>
</math></center>


a następnie rozwiązania sekwencji dwóch układów równań z macierzami trójkątnymi:
Numeryczna poprawność algorytmu jest wielce pożądaną cechą.


{{algorytm|Rozwiązanie układu z wykorzystaniem rozkładu LU||
<blockquote  style="background-color: #fefeee; padding:1em;  margin-left,margin-right:2em;  margin-top,margin-bottom: 1em;"> 
<pre>
Algorytm numerycznie poprawny działa w arytmetyce zmiennoprzecinkowej (niemal) tak, jakby wszystkie obliczenia były wykonywane w arytmetyce dokładnej, a tylko dane i wyniki podlegały reprezentacji w skończonej precyzji.
</blockquote>


Znajdź rozkład <math>\displaystyle A=LU</math>;
==Rola uwarunkowania zadania==
Rozwiąż <math>\displaystyle Ly = b</math> przez podstawienie w przód;
 
Rozwiąż <math>\displaystyle Ux = y</math> przez podstawienie w tył;
Niech <math>{\bf ALG}(\cdot)=\varphi(N(\cdot))</math> będzie algorytmem numerycznie
</pre>}}
poprawnym dla danych <math>F_0\subset F</math>. Wtedy jego błąd w <math>fl_\nu</math>
można oszacować następująco:
 
<center><math>\begin{align} \|S(f)-fl_\nu({\bf ALG}(f))\| \;=\;
    \|S(f)-\left(\varphi(y_\nu)\right)_\nu\|  \\
  &\le \|S(f)-\varphi(y)\|\,+\,
          \|\varphi(y)-\varphi(y_\nu)\|\,+\,
          \|\varphi(y_\nu)-\left(\varphi(y_\nu)\right)_\nu\| \\
  &\le \|S(f)-{\bf ALG}(f)\|\,+\,
          \|\varphi(y)-\varphi(y_\nu)\|\,+\,
          K_2\,\nu\,\|\varphi(y_\nu)\| \\
  &\le \|S(f)-{\bf ALG}(f)\|\,+\,
          (1 + K_2 \nu) \|\varphi(y)-\varphi(y_\nu)\|\,+\,
          K_2\,\nu\,\|\varphi(y)\|,
\end{align}</math></center>


Przypuśćmy, że taki rozkład <math>\displaystyle A=LU</math> istnieje. Wówczas, zapisując macierze w postaci
przy czym <math>\|y_\nu-y\|\,\le\,K_1\,\nu\,\|y\|</math>. Stąd
blokowej, eksponując pierwszy wiersz i pierwszą kolumnę zaangażowanych macierzy,
w szczególności wynika, że jeśli algorytm jest numerycznie
mamy
poprawny i ciągły ze względu na informację <math>y</math>, to


<center><math>\displaystyle
<center><math>\lim_{\nu\to 0}\,\|S(f)-fl_\nu({\bf ALG}(f))\|\,=\,
\begin{pmatrix}
      \|S(f)-{\bf ALG}(f)\|.
a_{11} & a_{12}^T\\
a_{21} & A_{22}
\end{pmatrix}
=
\begin{pmatrix}
1 & 0^T\\
l_{21} & L_{22}
\end{pmatrix}  
\begin{pmatrix}
u_{11} & u_{12}^T\\
0 & U_{22},
\end{pmatrix}  
</math></center>
</math></center>


skąd (mnożąc blokowo macierz <math>\displaystyle L</math> przez <math>\displaystyle U</math>) wynika, że
To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie
* <math>\displaystyle u_{11} = a_{11}</math> oraz <math>\displaystyle u_{12} = a_{12}</math>, więc pierwszy wiersz <math>\displaystyle U</math> jest
się zachowywał w <math>fl_\nu</math> prawie tak, jak w arytmetyce idealnej.
kopią pierwszego wiersza <math>\displaystyle A</math>,
 
* <math>\displaystyle l_{21} = a_{21}/u_{11}</math>, więc pierwsza kolumna <math>\displaystyle L</math> powstaje przez
Z powyższych wzorów wynika, że błąd w <math>fl_\nu</math> algorytmu
podzielenie wszystkich elementów wektora <math>\displaystyle a_{21}</math> przez element na diagonali
numerycznie poprawnego zależy w dużym stopniu od:
<math>\displaystyle a_{11}</math>,
* dokładności algorytmu w arytmetyce idealnej,  
* <math>\displaystyle A_{22} - l_{21}u_{12}^T = L_{22}U_{22}</math>, a więc znalezienie podmacierzy
* dokładności <math>\nu</math> arytmetyki <math>fl_\nu</math>,  
<math>\displaystyle L_{22}</math> oraz <math>\displaystyle U_{22}</math> sprowadza się do znalezienia rozkładu LU zmodyfikowanego
* wrażliwości algorytmu na małe względne zaburzenia informacji <math>y</math>.  
bloku <math>\displaystyle A_{22}</math> macierzy <math>\displaystyle A</math>,
wymiaru <math>\displaystyle (N-1)\times (N-1)</math>.
   
   
Dostajemy więc algorytm rekurencyjny, jednak ze względu na to, że wywołanie
Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy
rekurencji następuje na końcu każdej iteracji, można rozwinąć korzystająć z
trochę więcej uwagi jedynie trzeciemu.  
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
Jeśli <math>\varphi</math> spełnia warunek Lipschitza ze stałą <math>L</math>,
elementy <math>\displaystyle A</math> elementami macierzy <math>\displaystyle U</math> i <math>\displaystyle L</math> (jedynek z diagonali <math>\displaystyle L</math> nie musimy
a dokładniej
pamiętać, bo wiemy ''a priori'', że tam są).


{{algorytm|Rozkład LU metodą eliminacji Gaussa||
<center><math>\|\varphi(y_\nu)-\varphi(y)\|\,\le\,L\,\|y_\nu-y\|</math>,</center>
<pre>


for k=1:N-1
to
if <math>\displaystyle a_{kk}</math> == 0  
 
STOP;
<center><math>\begin{align} \|S(f)-fl_\nu({\bf ALG}(f))\| \\
end
    &\le \|S(f)-{\bf ALG}(f)\|\,+\,
for i=k+1:N /* wyznaczenie <math>\displaystyle k</math>-tej kolumny <math>\displaystyle L</math> */
      (1+K_2\nu)L\|y_\nu-y\|\,+\,K_2\nu\|\varphi(y)\|  \\
<math>\displaystyle a_{ik}</math> = <math>\displaystyle a_{ik}/a_{ii}</math>;
    &\le \|S(f)-{\bf ALG}(f)\|\,+\,
end
        (1+K_2\nu)LK_1\nu\|y\|\,+\,K_2\nu\|\varphi(y)\|.
for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
\end{align}</math></center>
for i=k+1:N
 
<math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny
end
algorytmu proporcjonalnie do <math>\nu</math>.
end
 
end
Bardziej jednak interesuje nas błąd <strong>względny</strong>. Wybierzmy
</pre>}}
"małe" <math>\eta\ge 0</math> i przypuśćmy, że
 
<center><math>\|\varphi(y_\nu)-\varphi(y)\|\;\le\;
    M\,K_1\,\nu\,\max(\eta,\|\varphi(y)\|)</math>,</center>
 
dla pewnej <math>M</math> niezależnej od <math>y</math>, tzn. błąd względny informacji,
<math>\|y_\nu-y\|\le K_1\nu\|y\|</math>, przenosi się na błąd względny
wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia"
<math>M</math>, albo na błąd bezwzględny ze współczynnikiem <math>M\eta</math>.
(Zauważmy, że gdybyśmy wzięli <math>\eta=0</math>, to dla <math>y</math> takiej, że
<math>\varphi(y)=0</math>, musiałoby być <math>\varphi(y_\nu)=0</math> --- co zwykle, choć
nie zawsze, nie jest prawdą.) Wtedy
 
<center><math>\begin{align} \|S(f)-fl_\nu({\bf ALG}(f))\|
  & \le \|S(f)-{\bf ALG}(f)\|+
    (1 + K_2 \nu) M K_1 \nu \max (\eta, \|\varphi(y)\|)+
      K_2 \nu \|\varphi(y)\| \\
    &= \|S(f)-{\bf ALG}(f)\|\,+\,\nu\,
        \Big(\,MK_1(1+K_2\nu)+K_2\Big)\max(\eta,\|\varphi(y)\|).
\end{align}</math></center>
 
W szczególności, gdy algorytm jest dokładny i korzysta z pełnej
informacji o <math>f</math>, tzn. <math>S\equiv{\bf ALG}\equiv\varphi</math>, to
błąd
 
<center><math>\frac{\|S(f)-fl_\nu({\bf ALG}(f))\|}
      {\max (\eta, \|S(f)\|)} \;\le\;
        \Big( M K_1 (1+K_2\nu) + K_2\Big)\,\nu
        \,\approx\,(M\,K_1\,+\,K_2)\,\nu. 
</math></center>
 
Stąd wynika, że jeśli <math>(MK_1+K_2)\nu\ll 1</math>, to błąd względny
algorytmu w <math>fl_\nu</math> jest mały, o ile <math>\|S(f)\|\ge\eta</math>.
Błąd względny jest proporcjonalny do dokładności <math>\nu</math>,
arytmetyki <math>fl_\nu</math>, współczynników proporcjonalności <math>K_i</math>
algorytmu numerycznie poprawnego, oraz do wrażliwości <math>M</math>
zadania <math>S</math> na małe względne zaburzenia danych.
 
Zwróćmy uwagę na istotny fakt, że interesują nas właściwie
tylko te zaburzenia danych (informacji), które powstają przy
analizie algorytmu numerycznie poprawnego. I tak, jeśli algorytm
jest numerycznie poprawny z pozornymi zaburzeniami danych w normie,
to trzeba zbadać wrażliwość zadania ze względu na zaburzenia
danych w normie. Jeśli zaś mamy pozorne zaburzenia
"po współrzędnych" (co oczywiście implikuje pozorne zaburzenia
w normie) to wystarczy zbadać uwarunkowanie zadania ze względu na
zaburzenia "po współrzędnych", itd.
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
<span  style="font-variant:small-caps;">Przykład: Iloczyn skalarny</span>
<div class="solution" style="margin-left,margin-right:3em;">
 
Załóżmy, że dla danych ciągów rzeczywistych o ustalonej
długości <math>n</math>, <math>a_j</math>, <math>b_j</math>, <math>1\le j\le n</math>, chcemy obliczyć


Łatwo przekonać się, że <math>\displaystyle k</math>-ty obrót zewnętrznej pętli (tzn. <math>\displaystyle k</math>-ty krok
<center><math>S(a,b)\,=\,\sum_{j=1}^n a_j b_j</math></center>
algorytmu rozkładu LU) kosztuje rzędu <math>\displaystyle 2(N-k)^2</math> operacji arytmetycznych, skąd łączny koszt tego
algorytmu rozkładu LU wynosi około <math>\displaystyle \frac{4}{3}N^3</math>.


Jeśli więc wykorzystać rozkład LU do rozwiązywania układu równań <math>\displaystyle Ax=b</math>, to mamy
Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem
następujące zestawienie kosztów:
i korzystający z pełnej informacji o kolejnych współrzednych.  
* Koszt znalezienia rozkładu <math>\displaystyle A=LU</math>: <math>\displaystyle O(N^3)</math>;
* Koszt rozwiązania układu <math>\displaystyle Ly=b</math>: <math>\displaystyle O(N^2)</math>;
* Koszt rozwiązania układu <math>\displaystyle Ux=y</math>: <math>\displaystyle O(N^2)</math>.
Tak więc, gdy znany już jest rozkład LU macierzy, koszt rozwiązania równania
wynosi już tylko <math>\displaystyle O(N^2)</math>.


{{uwaga|Złożoność obliczeniowa zadania rozwiązania układu równań liniowych||
Oznaczmy przez <math>\tilde a_j</math> i <math>\tilde b_j</math> reprezentacje liczb
<math>a_j</math> i <math>b_j</math> w <math>fl_\nu</math>, <math>\tilde a_j=a_j(1+\alpha_j)</math>,
<math>\tilde b_j=b_j(1+\beta_j)</math>, oraz przez <math>\gamma_j</math> i <math>\delta_j</math>
błędy względne powstałe przy kolejnych mnożeniach i dodawaniach.
Oczywiście <math>|\alpha_j|,|\beta_j|, |\gamma_j|, |\delta_j|\le\nu</math>.
Otrzymujemy


Z powyższego wynika, że łączny koszt rozwiązania równania liniowego wynosi
<center><math>\begin{align} fl_\nu\left(\sum_{j=1}^n a_jb_j\right) &=
<math>\displaystyle O(N^3)</math>. Można zastanawiać się, jaka jest <strong>najmniejsza możliwa</strong> liczba
    \Big(\,(fl_\nu(\sum_{j=1}^{n-1}a_jb_j)\,+\,\tilde a_n\tilde b_n
operacji zmiennoprzecinkowych potrzebnych do rozwiązania układu równań
      (1+\gamma_n)\,\Big)(1+\delta_n)\,=\,\ldots \\
liniowych.
  &= \bigg(\cdots\Big(
      \tilde a_1\tilde b_1(1+\gamma_1)+\tilde a_2\tilde b_2
        (1+\gamma_2)\Big)(1+\delta_2) +\cdots+
      \tilde a_n\tilde b_n(1+\gamma_n)\bigg)(1+\delta_n) \\
  &= \tilde a_1\tilde b_1(1+\gamma_1)(1+\delta_2)
                \cdots(1+\delta_n) +\cdots+\tilde a_j
      \tilde b_j(1+\gamma_j)(1+\delta_j)\cdots(1+\delta_n) \\
  &= \sum_{j=1}^n a_jb_j(1+e_j),
\end{align}</math></center>


Można pokazać,  że minimalny koszt rozwiązania układu <math>\displaystyle N</math> równań
gdzie w przybliżeniu (tzn. gdy <math>\nu\to 0</math>) mamy <math>|e_1|\leq (n+2)\nu</math>  
liniowych nie może być wyższego rzędu niż minimalny koszt mnożenia dwóch
i <math>|e_j|\leq (n-j+4)\nu</math>, <math>2\le j\le n</math>. Algorytm naturalny jest więc
macierzy <math>\displaystyle N\times N</math>. Tymczasem znany jest całkiem prosty algorytm rekurencyjny,
numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany
wyznaczający iloczyn dwóch macierzy kosztem <math>\displaystyle 4.7\cdot N^{log_27} \approx 4.7
w <math>fl_\nu</math> można zinterpretować jako dokładny wynik dla danych
\cdot N^{2.807}</math> (algorytm Strassena). Bardziej skomplikowany (i praktycznie
<math>a_{\nu,j}=a_j</math> i <math>b_{\nu,j}=b_j(1+e_j)</math>, przy czym
nieimplementowalny) algorytm Coppersmitha i Winograda daje nawet koszt
<math>\|b_\nu-b\|_p\leq (n+2)\nu\|b\|_p</math>.  
<math>\displaystyle O(N^{2.376})</math>.  Tak więc równania liniowe daje się (teoretycznie) rozwiązać
kosztem <math>\displaystyle O(N^{2.376})</math>.


Jednak w praktyce nawet prosty algorytm Strassena zazwyczaj nie jest stosowany.
Zobaczmy teraz, jak błąd we współrzędnych <math>b_j</math> wpływa
Wynika to stąd, że ma trochę gorsze własności numeryczne oraz, co istotniejsze, wymaga dużo
na błąd wyniku. Mamy
dodatkowej pamięci na przechowywanie pośrednich wyników.
}}


==Wybór elementu głównego==
<center><math>\begin{align} \Big|\sum_{j=1}^n a_jb_j-fl_\nu\Big(\sum_{j=1}^n a_jb_j\Big)\Big|
    &= \Big|\sum_{j=1}^n a_jb_j-\sum_{j=1}^n a_jb_j(1+e_j)\Big| \\
    &= \Big|\sum_{j=1}^n e_ja_jb_j\Big|
      \,\le\, \sum_{j=1}^n |e_j||a_jb_j| \\
    &\leq  (n+2)\nu\sum_{j=1}^n |a_jb_j|.
\end{align}</math></center>


Opisany powyżej algorytm rozkładu LU niestety czasem może się załamać: gdy
Stąd dla <math>\eta\ge 0</math>
napotka w czasie działania zerowy element na diagonali zmodyfikowanej
podmacierzy, np. chociaż macierz


<center><math>\displaystyle
<center><math>\frac{|\sum_{j=1}^n a_jb_j-fl_\nu(\sum_{j=1}^n a_jb_j)|}
A = \begin{pmatrix} 0 & 1\\ 1 & 0
      {\max(\eta,|\sum_{j=1}^n a_jb_j|)} \,\leq\,
\end{pmatrix}  
        K_{\eta}\,(n+2)\,\nu,
</math></center>
</math></center>


jest ewidentnie nieosobliwa, to nasz algorytm nawet nie ruszy z miejsca, bo od
gdzie
razu zetknie się z dzieleniem przez <math>\displaystyle a_{11}=0</math>. Ale wystarczy zamienić ze sobą
 
kolejnością wiersze macierzy <math>\displaystyle A</math> (to znaczy, w układzie równań, zamienić ze sobą
<center><math>K_\eta\,=\,K_\eta(a,b)\,=\,\frac{\sum_{j=1}^n |a_jb_j|}
miejscami równania), a dostajemy macierz, dla której algorytm przejdzie bez
            {\max(\eta,|\sum_{j=1}^n a_jb_j|) }.  
problemu.
</math></center>


W praktyce obliczeniowej, aby uzyskać algorytm o [[|Dodaj link: możliwie dobrych własnościach
Zauważmy, że jeśli iloczyny <math>a_jb_j</math> są wszystkie dodatnie
numerycznych]], wykorzystujemy tzw. strategię <strong>wyboru elementu głównego w
albo wszystkie ujemne, to <math>K_\eta=1</math>, tzn. zadanie jest dobrze
kolumnie</strong>.
uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej
Polega to na tym, że zanim wykonamy <math>\displaystyle k</math>-ty krok algorytmu rozkładu LU,
<math>n\nu</math>. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile
* szukamy w pierwszej kolumnie podmacierzy <math>\displaystyle A(k:N,k:N)</math> szukamy elementu o
liczba <math>n</math> składników nie jest horrendalnie duża. W ogólności
największym module (taki element, na mocy założenia nieosobliwości macierzy,
jednak <math>K_\eta</math> może być liczbą dowolnie dużą i wtedy nie możemy
jest niezerowy) --- to jest właśnie element główny
być pewni uzyskania dobrego wyniku w <math>fl_\nu</math>.
* zamieniamy ze sobą wiersz <math>\displaystyle A(k,1:N)</math> 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


<center><math>\displaystyle PA = LU,
</div></div>
</math></center>


gdzie <math>\displaystyle P</math> 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,
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
m.in. wybór w kolumnie oraz tzw. <strong>wybór pełny</strong>, gdy elementu głównego szukamy w
<span  style="font-variant:small-caps;">Przykład: Pierwiastki trójmianu</span>  
<strong>całej</strong> podmacierzy <math>\displaystyle A(k:N,k:N)</math>, co znacznie zwiększa liczbę porównań
<div class="solution" style="margin-left,margin-right:3em;">
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
Rozpatrzymy teraz zadanie obliczenia wszystkich pierwiastków
zerojedynkowej macierzy j.w. ale w jednym wektorze.
rzeczywistych równania kwadratowego.  
Będziemy zakładać, że model obliczeniowy dopuszcza obliczanie
pierwiastków kwadratowych z liczb nieujemnych oraz
<math>fl_\nu(\sqrt{x})=rd_\nu(\sqrt{rd_\nu(x)})</math>.  


{{algorytm|Rozkład LU z wyborem elementu głównego w kolumnie||
Okazuje się, że nie umiemy pokazać numerycznej poprawności
<pre>
"szkolnego" algorytmu obliczającego pierwiastki równania
bezpośrednio ze wzorów omawianych powyżej. Można jednak pokazać
numeryczną poprawność drobnej jego modyfikacji wykorzystującej
wzory Viete'a.


P = 1:N; /* tu zapamiętujemy wykonane permutacje */
{{algorytm|||
for k=1:N-1
<pre>\EATWSDelta = p*p - q;  
if  (Delta == 0)  
w wektorze A(k:N,k) znajdź element główny <math>\displaystyle a_{pk}</math>;
      OUT(p);
zamień ze sobą wiersze A(k,1:N) i A(p,1:N);
else
P(k) = p; P(p) = k;
if  (Delta > 0)
if <math>\displaystyle a_{kk}</math>
{
STOP: macierz osobliwa!
Delta1 = sqrt(d);  
end
if (p >= 0)
{
/* kontunuuj tak jak w algorytmie bez wyboru */
x1 = p + Delta1;
x2 = q/z1;
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>;
else
end
{  
for j=k+1:N /* modyfikacja podmacierzy <math>\displaystyle A_{22} = A_{22} - l_{21}u_{12}^T</math> */
x2 = p - Delta1;
for i=k+1:N
x1 = q/ź2;  
<math>\displaystyle a_{ij}</math> -= <math>\displaystyle a_{ik}a_{kj}</math>;
}
end
OUT(x1);  OUT(x2);
end
}
end
</pre>}}
</pre>}}


<div class="thumb tright"><div><flash>file=Macierz.swf</flash><div.thumbcaption>Eliminacja Gaussa z wyborem elementu głównego</div></div></div>
Mamy bowiem


Poniżej zapisujemy algorytm rozwiązania układu równań, gdy znany jest
<center><math>\begin{align} fl_\nu(\Delta(p,q)) &= \Big(p^2(1+\alpha)^2(1+\epsilon_1)-q(1+\beta)\Big)
już rozkład LU wykonany z wyborem w kolumnie.
                          (1+\epsilon_2) \\
    &= \left( p^2-q\frac{(1+\beta)}{(1+\alpha)^2(1+\epsilon_1)}\right)
          (1+\epsilon_2)(1+\alpha)^2(1+\epsilon_1) \\
    &= \Big(p^2-q(1+\delta)\,\Big)(1+\gamma) \,=\,
          \Delta(p,q(1+\delta))(1+\gamma),
\end{align}</math></center>


{{algorytm|Rozwiązywanie układu równań metodą eliminacji Gaussa z wyborem elementu głównego w kolumnie||
gdzie <math>|\delta|,|\gamma|\leq 4\nu</math>. Wyróżnik obliczony w <math>fl_\nu</math>
<pre>
jest więc nieco zaburzonym wyróżnikiem dokładnym dla danych
<math>p</math> i <math>q_\nu=q(1+\delta)</math>. W szczególności


znajdź rozkład <math>\displaystyle PA = LU</math>;
<center><math>\mbox{sgn} (fl_\nu(\Delta(p,q)))= \mbox{sgn} (\Delta(p,q_\nu))</math></center>
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>;
Jeśli <math>p\ge 0</math> to
</pre>}}
 
<center><math>\begin{align} fl_\nu(x1(p,q)) &= \Big(p(1+\alpha)+
        \sqrt{fl_\nu(\Delta(p,q))}(1+\epsilon_3)\Big)(1+\epsilon_4) \\
  &= \Big(p(1+\alpha)+\sqrt{\Delta(p,q_\nu)
      (1+\gamma)}(1+\epsilon_3)\Big)(1+\epsilon_4) \\
  &= \left(p+\sqrt{\Delta(p,q_\nu)}\frac{\sqrt{1+\gamma}(1+\epsilon_3)}
        {1+\alpha}\right)(1+\epsilon_4)(1+\alpha) \\
  &= \left(p+\sqrt{\Delta(p,q_\nu)}\right)(1+e_1),
\end{align}</math></center>


Dla niektórych ważnych klas macierzy wiadomo, że rozkład LU jest wykonalny <strong>bez wyznaczania elementu głównego</strong>, co istotnie może zmniejszyć całkowity czas
gdzie <math>|e_1|\leq 6\nu</math>. Zauważmy, że ostatnia równość
działania algorytmu. Jest tak m.in. dla macierzy
zachodzi dlatego, że dodajemy liczby tego samego znaku. (Inaczej
* symetrycznych, dodatnio określonych: <math>\displaystyle A=A^T</math> oraz <math>\displaystyle x^TAx > 0</math>, <math>\displaystyle \forall x
<math>|e_1|</math> mogłaby być dowolnie duża i tak byłoby w algorytmie
\neq 0</math>,
szkolnym.) Dla drugiego pierwiastka mamy
* silnie diagonalnie dominujących: macierz <math>\displaystyle A</math> (lub <math>\displaystyle A^T</math>) spełnia


<center><math>\displaystyle |a_{ii}| > \sum_{j\neq i} |a_{ij}|, \qquad \forall i.
<center><math>fl_\nu(x2(p,q))\,=\,\frac {q(1+\beta)}{fl_\nu(x1(p,q))}(1+\epsilon_5)
  \,=\,\frac{q_\nu}{fl_\nu(x1(p,q))}(1+e_2),  
</math></center>
</math></center>
gdzie <math>|e_2|\le 8\nu</math>.
Podobny wynik otrzymalibyśmy dla <math>p<0</math>. Algorytm zmodyfikowany
jest więc numerycznie poprawny, gdyż otrzymane w <math>fl_\nu</math> pierwiastki
są nieco zaburzonymi dokładnymi pierwiatkami dla danych
<math>p_\nu=p</math> i <math>q_\nu=q(1+\delta)</math>.
Aby oszacować błąd algorytmu, wystarczy zbadać uwarunkowanie
zadania ze względu na zaburzenie danej <math>q</math>, ponieważ pokazaliśmy,
że zaburzenia <math>p</math> można przenieść na zaburzenia <math>q</math> i wyniku.
Niestety, choć algorytm jest numerycznie poprawny, zaburzenia
<math>q</math> mogą sprawić, że nawet znak wyróżnika <math>\Delta</math> może być
obliczony nieprawidłowo. Na przykład dla <math>p=1</math> i <math>q=1\pm 10^{t+1}</math>
mamy <math>\Delta(p,q)=\mp 10^{t+1}</math>, ale
<math>\Delta(rd_\nu(p),rd_\nu(q))=\Delta(1,1)=0</math>. Ogólnie
<center><math>|fl_\nu(\Delta(p,q))-\Delta(p,q)|\,\leq\,4\nu(p^2+2|q|)</math>,</center>
a więc tylko dla <math>|\Delta(p,q)|=|p^2-q|>4\nu (p^2+2|q|)</math>
możemy być pewni obliczenia właściwego znaku <math>\Delta</math>. Przy
tym warunku oraz <math>\Delta>0</math> błąd danych przenosi się w
normie euklidesowej na błąd wyniku następująco:
<center><math>\begin{align} \lefteqn{ \Big( (x1(p,q) - x1(p,q_\nu))^2
                +(x2(p,q) - x2(p,q_\nu))^2 \Big)^{1/2} } \\
  &= \frac{\sqrt 2 |\delta q|} {\sqrt{p^2-q}+\sqrt{p^2-q_\nu}}
  \,\leq\, 2\sqrt 2 \nu\frac{|q|}{\sqrt{p^2-q}} \\
  &= 2\sqrt 2 \nu\frac{|q|/p^2}{\sqrt{1-q/p^2}
        \max(\eta/|p|,\sqrt{2(1+(1-q/p^2))}) } \\
  & & \qquad\qquad\qquad\cdot\max(\eta,(x1(p,q)^2+x2(p,q)^2)^{1/2}).
\end{align}</math></center>
Stąd widać, że zadanie jest dobrze uwarunkowane dla <math>q/p^2\ll 1</math>
i może być źle uwarunkowane dla <math>q/p^2\approx 1</math>. W ostatnim
przypadku nie możemy być pewni otrzymania dobrego wyniku w <math>fl_\nu</math>.
</div></div>
-->
==Literatura==
W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 2.3</b> w
* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.
Warto także przejrzeć rozdział 2 w
* <span style="font-variant:small-caps">P. Deulfhard, A. Hohmann</span>, <cite>Numerical Analysis in Modern Scientific Computing</cite>, Springer, 2003,
omawiający zagadnienia uwarunkowania i numerycznej poprawności algorytmów.
Nieocenioną monografią na ten temat jest
* <span style="font-variant:small-caps">N. Higham</span>, <cite>Accuracy and Stability of Numerical Algorithms</cite>, SIAM, 2002.

Aktualna wersja na dzień 21:46, 11 wrz 2023


Własności zadania obliczeniowego i algorytmu numerycznego

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

Uwarunkowanie zadania obliczeniowego

Jak zobaczyliśmy w poprzednich przykładach, dane, jakimi dysponujemy wykonując zadanie obliczeniowe, są z natury rzeczy wartościami zaburzonymi. Źródłem tych zaburzeń są:

  • błąd reprezentacji danych w arytmetyce zmiennoprzecinkowej (np. 0.1 nie jest równe dokładnie 1/10)
  • błędy w parametrach będących rezultatem poprzednich obliczeń (np. chcemy rozwiązać równanie f(x)=a, ale a jest rezultatem innej symulacji), a także
  • błędy pomiarowe w zadaniach praktycznych (np. chcemy policzyć numeryczną prognozę pogody, ale temperaturę wyjściową znamy tylko z dokładnością do 0.1 stopnia, co gorsza --- z niektórych stacji w ogóle nie mamy danych)

Okazuje się, że powszechna intuicja, że małe zaburzenia danych powinny dawać małe zaburzenia wyniku, nie znajduje potwierdzenia nawet w bardzo prostych przypadkach. Z drugiej strony, umiejętność oceny jakościowego wpływu zaburzenia danych na wynik jest kapitalna w świecie obliczeń numerycznych w ogólności, a w szególności --- inżynierskich.

Wprowadza się pojęcie uwarunkowania zadania, to znaczy jego podatności na zaburzenia danych. Dla przejrzystości przypuśćmy, że nasze zadanie obliczeniowe polega na wyznaczeniu f(x) dla danego x.


Jak bardzo będzie odległe f(x~), gdy x~x? Rozważa się dwa przypadki:

  • uwarunkowanie względne: jak względne zaburzenie danych wpływa na błąd względny wyniku:
    ||f(x)f(x~)||||f(x)||condrel(f,x)||xx~||||x||

Najmniejszy mnożnik condrel(f,x) spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (względnego) zadania obliczenia f(x) dla danego x.

  • uwarunkowanie bezwzględne: jak bezwzględne zaburzenie danych wpływa na błąd bezwzględny wyniku:
    ||f(x)f(x~)||condabs(f,x)||xx~||

Najmniejszy mnożnik condabs(f,x) spełniający powyższą nierówność nazywamy współczynnikiem uwarunkowania (bezwzględnego) zadania obliczenia f(x) dla danego x.

Powiemy, że zadanie f(x) jest

  • dobrze uwarunkowane w punkcie x, gdy cond(f,x)1,
  • źle uwarunkowane w punkcie x, gdy cond(f,x)1,
  • źle postawione w punkcie x, gdy cond(f,x)=+.

Teraz już rozumiemy, że redukcja cyfr przy odejmowaniu jest tylko odzwierciedleniem faktu, że zadanie odejmowania dwóch bliskich liczb jest po prostu zadaniem źle uwarunkowanym!

Przykład: Uwarunkowanie zadania obliczenia sumy

Właściwie już wcześniej sprawdziliśmy, że zadanie obliczenia s(x,y)=x+y ma

condabs(s,(a,b))=1,condrel(s,(a,b))=|a|+|b||a+b|

Tak więc, gdy ab, to condrel(s,(a,b))+ i zadanie jest bardzo źle uwarunkowane. Nawet małe zaburzenie względne danych może skutkować bardzo dużym zaburzeniem wyniku. Zgodnie z prawem Murphy'ego, najczęściej rzeczywiście tak będzie...

Przykład

Dla różniczkowalnej funkcji skalarnej f:RR mamy

|f(x)f(x~)||f(x)||xx~|

i w konsekwencji dla zadania obliczenia f(x) dla danego x mamy, przy założeniu małych zaburzeń,

condabs(f,x)=|f(x)|,condrel(f,x)=|f(x)||x||f(x)|

Jest zrozumiałe, że złe uwarunkowanie jest szkodliwe w praktyce numerycznej: zaburzenia danych są nieuniknione, a źle uwarunkowane zadanie tylko je wzmocni na wyjściu. Jednak za jakiś czas zobaczysz wyjątkową sytuację, gdy złe uwarunkowanie pewnego podzadania nie tylko nie przeszkadza, ale wręcz pomaga szybciej rozwiązać zadanie główne!

Rozkład algorytmu względem informacji

Algorytm to dokładnie określona i dozwolona w danym modelu obliczeniowym sekwencja akcji, pozwalająca na rozwiązanie danego zadania (w sposób dokładny lub przybliżony).

Z każdym algorytmem związany jest operator

𝐀𝐋𝐆:FG,

taki że 𝐀𝐋𝐆(f) jest wynikiem działania algorytmu w arytmetyce idealnej dla danej f.

Zauważmy, że wynik 𝐀𝐋𝐆(f) działania algorytmu nie zależy bezpośrednio od f, ale raczej od informacji o f (uzyskanej dzięki poleceniu 𝒩). Informacja ta może być pełna albo tylko częściowa. Informacja jest pełna gdy, np. f=(f1,,fn)Rn i wczytamy wszystkie współrzędne fi. Informacja może być częściowa, gdy f jest funkcją. Wtedy wiele danych może posiadać tę samą informację, co łatwo zaobserwować na przykładzie zadania całkowania.

Niech N:Fn=0Rn będzie operatorem informacji, tzn.

N(f)=(y1,y2,,yn)

jest informacją o f zebraną przy idealnej realizacji algorytmu. Zauważmy, że informacja jest pełna, gdy N jest przekształceniem różnowartościowym, tzn. jeśli f1f2 implikuje N(f1)N(f2). W przeciwnym przypadku mamy do czynienia z informacją częściową.

Każdy algorytm 𝐀𝐋𝐆 może być przedstawiony jako złożenie operatora informacji i pewnego operatora φ:N(F)G zdefiniowanego równością

φ(N(f))=𝐀𝐋𝐆(f).

Zauważmy, że w przypadku informacji częściowej zwykle nie istnieje algorytm dający dokładne rozwiązanie zadania dla każdej danej fF, ponieważ dla danych o tej samej informacji mogą istnieć różne rozwiązania.

Problem wyboru algorytmu

Wybór algorytmu jest najistotniejszą częścią całego procesu numerycznego rozwiązywania zadania. Kierujemy się przy tym przede wszystkim następującymi kryteriami:

  • dokładnością algorytmu,
  • złożonością algorytmu,
  • własnościami numerycznymi algorytmu.

Przez dokładność algorytmu rozumiemy różnicę między rozwiązaniem dokładnym S(f) a rozwiązaniem 𝐀𝐋𝐆(f) dawanym przez algorytm w arytmetyce idealnej. Jeśli 𝐀𝐋𝐆(f)=S(f), fF, to algorytm nazywamy dokładnym.

Mówiąc o złożoności, mamy na myśli złożoność pamięciową (zwykle jest to liczba stałych i zmiennych używanych przez algorytm), jak również złożoność obliczeniową. Na złożoność obliczeniową algorytmu dla danej f składa się koszt uzyskania infomacji y=N(f) (zwykle jest on proporcjonalny do liczby wywołań polecenia 𝒩), oraz koszt kombinatoryczny przetworzenia tej informacji, aż do uzyskania wyniku φ(y). Koszt kombinatoryczny zwykle mierzymy liczbą operacji arytmetycznych wykonywanych przez algorytm.

Przez własności numeryczne algorytmu rozumiemy jego własności przy realizacji w arytmetyce flν. Temu ważnemu tematowi poświęcimy teraz osobny paragraf.

Numeryczna poprawność algorytmu

Pożądane jest, aby algorytm dawał "dobry" wynik zarówno w arytmetyce idealnej, jak i w arytmetyce flν. Niestety, jak zobaczymy, nie zawsze jest to możliwe. Nawet jeśli algorytm jest dokładny, to w wyniku jego realizacji w flν możemy otrzymać wynik flν(𝐀𝐋𝐆(f)) daleko odbiegający od S(f). W szczególności, prawie zawsze mamy

S(f)flν(𝐀𝐋𝐆(f)).

Zauważmy również, że o ile z reguły znamy dokładne zachowanie się algorytmu w arytmetyce idealnej dla danej informacji, to nie można tego samego powiedzieć o jego zachowaniu się w arytmetyce flν. W związku z tym powstaje pytanie, jak kontrolować błąd algorytmów wynikający z błędów zaokrągleń i jakie algorytmy uznamy za te o najwyższej jakości numerycznej.

Istnienie błędów reprezentacji liczb rzeczywistych powoduje, że informacja y=N(f) o danej f nie jest w ogólności reprezentowana dokładnie. Znaczy to, że zamiast na informacji dokładnej, dowolny algorytm będzie operować na informacji nieco zaburzonej yν, tzn. zaburzonej na poziomie błędu reprezentacji. Tak samo wynik dawany przez algorytm będzie w ogólności zaburzony na poziomie błędu reprezentacji. W najlepszym więc wypadku wynikiem działania algorytmu w flν będzie (φ(yν))ν zamiast φ(y). Algorytmy dające tego rodzaju wyniki uznamy za posiadające najlepsze własności numeryczne w arytmetyce flν i nazwiemy numerycznie poprawnymi.

Powiemy, że ciąg rzeczywisty aν=(aν,1,,aν,n) (a właściwie rodzina ciągów {aν}ν) jest nieco zaburzonym ciągiem a=(a1,,an), jeśli istnieje stała K taka, że dla wszystkich dostatecznie małych ν zachodzi

|aν,jaj|Kν|aj|,1jn,

albo ogólniej

aνaKνa,

gdzie jest pewną normą w Rn. W pierwszym przypadku mówimy o zaburzeniu "po współrzędnych", a w drugim o zaburzeniu w normie .

Zauważmy, że niewielkie zaburzenia po współrzędnych pociągają za sobą niewielkie zaburzenia w normie. Rzeczywiście, zachodzi wówczas

aνa=max1jn|aν,jaj|Kνmax1jn|aj|=Kνa,

i korzystając z faktu, że w przestrzeni skończenie wymiarowej wszystkie normy są równoważne, otrzymujemy dla pewnych stałych K1 i K2

aνaK1aνaK1KνaK2K1Kνa,

czyli nierówność dla zaburzenia w normie, ze stałą K=K2K1K.

Definicja Algorytm numerycznie poprawny

Algorytm 𝐀𝐋𝐆 rozwiązywania zadania nazywamy numerycznie poprawnym w zbiorze danych F0F, jeśli dla każdej danej fF0 wynik flν(𝐀𝐋𝐆(f)) działania algorytmu w arytmetyce flν można zinterpretować jako nieco zaburzony wynik algorytmu w arytmetyce idealnej dla nieco zaburzonej informacji yν=(N(f))νN(F) o f, przy czym poziom zaburzeń nie zależy od f.

Formalnie znaczy to, że istnieją stałe K1, K2 oraz ν0>0 takie, że spełniony jest następujący warunek. Dla dowolnej νν0 oraz informacji yN(F0) można dobrać yνN(F) oraz (φ(yν))ν takie, że

yνyK1νy,
(φ(yν))νφ(yν)K2νφ(yν),

oraz

flν(𝐀𝐋𝐆(f))=flν(φ(N(f)))=(φ(yν))ν.
Numerycznie poprawny algorytm daje w arytmetyce flν wynik ALG(N(x)), który daje się zinterpretować jako mało zaburzony wynik f(y) zadania na mało zaburzonych danych x.

Zauważmy,że jeśli fRn, N(f)=(f1,,fn), oraz algorytm jest dokładny, 𝐀𝐋𝐆φS, to numeryczną poprawność algorytmu można równoważnie zapisać jako

flν(𝐀𝐋𝐆(f))=(S(fν))ν.

Numeryczna poprawność algorytmu jest wielce pożądaną cechą.

Algorytm numerycznie poprawny działa w arytmetyce zmiennoprzecinkowej (niemal) tak, jakby wszystkie obliczenia były wykonywane w arytmetyce dokładnej, a tylko dane i wyniki podlegały reprezentacji w skończonej precyzji.

Rola uwarunkowania zadania

Niech 𝐀𝐋𝐆()=φ(N()) będzie algorytmem numerycznie poprawnym dla danych F0F. Wtedy jego błąd w flν można oszacować następująco:

S(f)flν(𝐀𝐋𝐆(f))=S(f)(φ(yν))νS(f)φ(y)+φ(y)φ(yν)+φ(yν)(φ(yν))νS(f)𝐀𝐋𝐆(f)+φ(y)φ(yν)+K2νφ(yν)S(f)𝐀𝐋𝐆(f)+(1+K2ν)φ(y)φ(yν)+K2νφ(y),

przy czym yνyK1νy. Stąd w szczególności wynika, że jeśli algorytm jest numerycznie poprawny i ciągły ze względu na informację y, to

limν0S(f)flν(𝐀𝐋𝐆(f))=S(f)𝐀𝐋𝐆(f).

To znaczy, że dla dostatecznie silnej arytmetyki algorytm będzie się zachowywał w flν prawie tak, jak w arytmetyce idealnej.

Z powyższych wzorów wynika, że błąd w flν algorytmu numerycznie poprawnego zależy w dużym stopniu od:

  • dokładności algorytmu w arytmetyce idealnej,
  • dokładności ν arytmetyki flν,
  • wrażliwości algorytmu na małe względne zaburzenia informacji y.

Ponieważ dwa pierwsze punkty są raczej oczywiste, poświęcimy trochę więcej uwagi jedynie trzeciemu.

Jeśli φ spełnia warunek Lipschitza ze stałą L, a dokładniej

φ(yν)φ(y)Lyνy,

to

S(f)flν(𝐀𝐋𝐆(f))S(f)𝐀𝐋𝐆(f)+(1+K2ν)Lyνy+K2νφ(y)S(f)𝐀𝐋𝐆(f)+(1+K2ν)LK1νy+K2νφ(y).

W tym przypadku błędy zaokrągleń zwiększają błąd bezwzględny algorytmu proporcjonalnie do ν.

Bardziej jednak interesuje nas błąd względny. Wybierzmy "małe" η0 i przypuśćmy, że

φ(yν)φ(y)MK1νmax(η,φ(y)),

dla pewnej M niezależnej od y, tzn. błąd względny informacji, yνyK1νy, przenosi się na błąd względny wyniku (w arytmetyce idealnej) ze "współczynnikiem wzmocnienia" M, albo na błąd bezwzględny ze współczynnikiem Mη. (Zauważmy, że gdybyśmy wzięli η=0, to dla y takiej, że φ(y)=0, musiałoby być φ(yν)=0 --- co zwykle, choć nie zawsze, nie jest prawdą.) Wtedy

S(f)flν(𝐀𝐋𝐆(f))S(f)𝐀𝐋𝐆(f)+(1+K2ν)MK1νmax(η,φ(y))+K2νφ(y)=S(f)𝐀𝐋𝐆(f)+ν(MK1(1+K2ν)+K2)max(η,φ(y)).

W szczególności, gdy algorytm jest dokładny i korzysta z pełnej informacji o f, tzn. S𝐀𝐋𝐆φ, to błąd

S(f)flν(𝐀𝐋𝐆(f))max(η,S(f))(MK1(1+K2ν)+K2)ν(MK1+K2)ν.

Stąd wynika, że jeśli (MK1+K2)ν1, to błąd względny algorytmu w flν jest mały, o ile S(f)η. Błąd względny jest proporcjonalny do dokładności ν, arytmetyki flν, współczynników proporcjonalności Ki algorytmu numerycznie poprawnego, oraz do wrażliwości M zadania S na małe względne zaburzenia danych.

Zwróćmy uwagę na istotny fakt, że interesują nas właściwie tylko te zaburzenia danych (informacji), które powstają przy analizie algorytmu numerycznie poprawnego. I tak, jeśli algorytm jest numerycznie poprawny z pozornymi zaburzeniami danych w normie, to trzeba zbadać wrażliwość zadania ze względu na zaburzenia danych w normie. Jeśli zaś mamy pozorne zaburzenia "po współrzędnych" (co oczywiście implikuje pozorne zaburzenia w normie) to wystarczy zbadać uwarunkowanie zadania ze względu na zaburzenia "po współrzędnych", itd.

Przykład: Iloczyn skalarny

Załóżmy, że dla danych ciągów rzeczywistych o ustalonej długości n, aj, bj, 1jn, chcemy obliczyć

S(a,b)=j=1najbj

Rozpatrzymy algorytm dokładny zdefiniowany powyższym wzorem i korzystający z pełnej informacji o kolejnych współrzednych.

Oznaczmy przez a~j i b~j reprezentacje liczb aj i bj w flν, a~j=aj(1+αj), b~j=bj(1+βj), oraz przez γj i δj błędy względne powstałe przy kolejnych mnożeniach i dodawaniach. Oczywiście |αj|,|βj|,|γj|,|δj|ν. Otrzymujemy

flν(j=1najbj)=((flν(j=1n1ajbj)+a~nb~n(1+γn))(1+δn)==((a~1b~1(1+γ1)+a~2b~2(1+γ2))(1+δ2)++a~nb~n(1+γn))(1+δn)=a~1b~1(1+γ1)(1+δ2)(1+δn)++a~jb~j(1+γj)(1+δj)(1+δn)=j=1najbj(1+ej),

gdzie w przybliżeniu (tzn. gdy ν0) mamy |e1|(n+2)ν i |ej|(nj+4)ν, 2jn. Algorytm naturalny jest więc numerycznie poprawny w całym zbiorze danych, gdyż wynik otrzymany w flν można zinterpretować jako dokładny wynik dla danych aν,j=aj i bν,j=bj(1+ej), przy czym bνbp(n+2)νbp.

Zobaczmy teraz, jak błąd we współrzędnych bj wpływa na błąd wyniku. Mamy

|j=1najbjflν(j=1najbj)|=|j=1najbjj=1najbj(1+ej)|=|j=1nejajbj|j=1n|ej||ajbj|(n+2)νj=1n|ajbj|.

Stąd dla η0

|j=1najbjflν(j=1najbj)|max(η,|j=1najbj|)Kη(n+2)ν,

gdzie

Kη=Kη(a,b)=j=1n|ajbj|max(η,|j=1najbj|).

Zauważmy, że jeśli iloczyny ajbj są wszystkie dodatnie albo wszystkie ujemne, to Kη=1, tzn. zadanie jest dobrze uwarunkowane, a błąd względny jest zawsze na poziomie co najwyżej nν. W tym przypadku algorytm zachowuje się bardzo dobrze, o ile liczba n składników nie jest horrendalnie duża. W ogólności jednak Kη może być liczbą dowolnie dużą i wtedy nie możemy być pewni uzyskania dobrego wyniku w flν.


Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 2.3 w

  • D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.

Warto także przejrzeć rozdział 2 w

  • P. Deulfhard, A. Hohmann, Numerical Analysis in Modern Scientific Computing, Springer, 2003,

omawiający zagadnienia uwarunkowania i numerycznej poprawności algorytmów. Nieocenioną monografią na ten temat jest

  • N. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 2002.