Zagadnienia

10.1 Twierdzenie o istnieniu rozwiązań
10.2 Obniżanie kosztu iteracji

10. Wariacje na temat metody Newtona

Podstawowy element metody Newtona — rozwiązywanie równania zlinearyzowanego — w przypadku, gdy $N$ jest duże, może stanowić wąskie gardło całego procesu iteracyjnego. Dlatego w tym i następnym rozdziale poszukamy skutecznych metod ominięcia tego ograniczenia; jednak na początek przytoczymy inną wersję twierdzenia o zbieżności, która (przy silniejszych założeniach) zagwarantuje nam także istnienie rozwiązań.

10.1. Twierdzenie o istnieniu rozwiązań

W dotychczasowej analizie zakładaliśmy, że rozwiązanie $x^{*}$ równania $F(x)=0$ istnieje. Okazuje się, że nieco modyfikując założenia standardowe można udowodnić i istnienie rozwiązania, i lokalną zbieżność metody Newtona do tego rozwiązania.

Twierdzenie 10.1 (Kantorowicza, o istnieniu i zbieżności)

Niech $F:\mathbb{R}^{N}\supset D\rightarrow\mathbb{R}^{N}$ , przy czym $D$ jest niepusty, otwarty i wypukły, i niech $F$ będzie różniczkowalna w $D$ oraz istnieje $L>0$ takie, że

$\| F^{{\prime}}(x)-F^{{\prime}}(y)\|\leq L\| x-y\|\qquad\forall x,y\in D.$

Niech istnieje $x_{0}\in D$ takie, że $F^{{\prime}}(x_{0})$ jest nieosobliwa oraz dla pewnych dodatnich $\beta,\eta$ zachodzi

$\| F^{{\prime}}(x_{0})^{{-1}}F(x_{0})\|\leq\beta,\qquad\| F^{{\prime}}(x_{0})^{{-1}}\|\leq\eta,$

przy czym $\alpha=L\beta\eta\leq 1/2$ . Określmy $\delta=(1-\sqrt{1-2\alpha})/(L\beta)$ i załóżmy, że $\bar{K}(x_{0},\delta)\subset D$ .

Wówczas ciąg $x_{k}$ zdefiniowany metodą Newtona jest dobrze określony, co więcej $x_{k}\in\bar{K}(x_{0},\delta)$ oraz ma granicę $x^{*}$ , która jest jedynym miejscem zerowym $F$ w $\bar{K}(x_{0},\delta)$ . Ponadto

$\| x_{k}-x^{*}\|\leq(L\beta 2^{k})(2\alpha)^{{2^{k}}}$

dla $k\geq 0$ .

Dowód

Dowód tego twierdzenia można znaleźć np. w [12, rozdział 12.6.1].

∎

10.2. Obniżanie kosztu iteracji

Z punktu widzenia praktycznej realizacji, każdy krok metody Newtona ma parę potencjalnie kosztownych momentów, które teraz przedyskutujemy.

Wyznaczenie wzoru na $F^{{\prime}}(x)$ może być trudne: na przykład, gdy $F$ jest zadana bardzo skomplikowaną procedurą (a nie jednolinijkowym wzorem), wtedy ani ręczne, ani automatyczne różniczkowanie $F$ może nie prowadzić do poprawnych czy zadowalających rezultatów. Może też być bardzo żmudne.
Wyznaczenie wartości $F(x)$ , a zwłaszcza: $F^{{\prime}}(x)$ , może być kosztowne. Rzeczywiście, $F^{{\prime}}(x)$ to $N^{2}$ elementów macierzy pochodnej, a dodatkowo — poza najprostszymi przypadkami — pochodna jest zwykle dana bardziej skomplikowanym wyrażeniem niż sama funkcja. Rozsądne jest więc przyjęcie¹³O ile $F^{{\prime}}(x)$ nie jest macierzą rozrzedzoną., że w praktyce koszt wyznaczenia $F^{{\prime}}(x)$ będzie rzędu $c\,\cdot\, N\,\cdot\,\text{(koszt obliczenia }F(x))$ , z raczej dużą stałą $c$ . Nawet wyznaczenie pojedynczej wartości $F(x)$ może być w niektórych zadaniach bardzo kosztowne. W [3] wspomina się o $F:\mathbb{R}^{{50}}\rightarrow\mathbb{R}^{{50}}$ , której jedna wartość kosztowała 100 godzin pracy ówczesnego szybkiego komputera.
Rozwiązanie układu z $F^{{\prime}}(x)$ jest potencjalnie najbardziej kosztownym fragmentem iteracji. Rzeczywiście, gdy $F^{{\prime}}(x)$ jest macierzą gęstą bez żadnych szczególnych własności, to koszt rozwiązania układu

$F^{{\prime}}(x)s=F(x)$

jest rzędu $O(N^{3})$ .

Te obserwacje prowadzą do kilku heurystycznych sposobów modyfikacji metody Newtona tak, by obniżyć koszt jednej jej iteracji. Czasem — choć nie zawsze, i dlatego te sposoby są ważne! — obniżenie kosztu iteracji będzie skutkowało pogorszeniem szybkości zbieżności metody, ale nawet i wówczas może okazać się, że ostatecznie metoda wolniej zbieżna, ale tańsza, będzie bardziej efektywna od metody Newtona.

10.2.1. Uproszczona metoda Newtona

W uproszczonej metodzie Newtona, najdroższy fragment iteracji Newtona — wyznaczenie rozkładu trójkątnego macierzy pochodnej — usuwamy poza pętlę:

$x_{{k+1}}=x_{k}-F^{{\prime}}(x_{0})^{{-1}}F(x_{k}).$

Rzeczywiście, jeśli bowiem na samym początku iteracji dokonamy kosztem $O(N^{3})$ flopów rozkładu LU (lub QR) macierzy $F^{{\prime}}(x_{0})$ i zapamiętamy czynniki rozkładu, to potem każdy z układów równań:

$F^{{\prime}}(x_{0})s=F(x_{k})$

będziemy mogli rozwiązać już tylko kosztem nie wyższym niż $O(N^{2})$ ! Oznacza to, że koszt jednej iteracji spada nam do $O(N^{2})$ , jest więc $N$ -krotnie mniejszy niż w przypadku oryginalnej metody Newtona. Ceną za przyspieszenie jest wyraźne zmniejszenie szybkości zbieżności metody:

Twierdzenie 10.2 (o zbieżności uproszczonej metody Newtona)

Przy standardowych założeniach, istnieją stałe dodatnie $C$ i $\delta$ takie, że jeśli $x_{0}\in K(x^{*},\delta)$ , to uproszczona metoda Newtona jest zbieżna przynajmniej liniowo do $x^{*}$ oraz

$\displaystyle\| x_{{k+1}}-x^{*}\|\leq C\| x_{{0}}-x^{*}\|\| x_{{k}}-x^{*}\|\qquad\forall k\in\mathbb{N}.$

(10.1)

Dowód

Zostanie podany później, gdy będziemy dysponowali wygodnym narzędziem do badania tego typu metod.

∎

Ćwiczenie 10.1

Zaimplementuj uproszczoną metodę Newtona.

Rozwiązanie:

Wystarczy zmodyfikować algorytm realizujący metodę Newtona:

Uproszczona metoda Newtona

function simplerNewton(x, F, stop)

oblicz macierz pochodnej $F^{{\prime}}(x)$ ;

wyznacz rozkład $F^{{\prime}}(x)$ , np. $QR=F^{{\prime}}(x)$ ;

while not stop do begin

rozwiąż układ z macierzą trójkątną, $Rs=Q^{T}F(x)$ ;

x = x - s;

end

return(x);

Ćwiczenie 10.2

Metoda Szamańskiego to metoda, w której po $m$ krokach uproszczonej metody Newtona dokonujemy wyznaczenia i rozkładu macierzy pochodnej (czyli restartujemy metodę uproszczoną, biorąc za przybliżenie początkowe ostatnio wyznaczone przybliżenie rozwiązania). Jest to więc coś pośredniego pomiędzy prawdziwą metodą Newtona (w której w każdej iteracji wyznacza się i rozkłada macierz pochodnej) a uproszczoną metodą Newtona (gdzie macierz pochodnej wyznacza się tylko jeden raz).

Formalnie możemy więc przyjąć, że $m$ kroków uproszczonej metody Newtona to jest ,,jeden duży krok” i określić iterację Szamańskiego następująco:

	$\displaystyle x_{{k+\frac{1}{m}}}$	$\displaystyle=x_{k}-F^{{\prime}}(x_{k})^{{-1}}F(x_{k})$
	$\displaystyle x_{{k+\frac{j+1}{m}}}$	$\displaystyle=x_{{k+\frac{j}{m}}}-F^{{\prime}}(x_{k})^{{-1}}F(x_{{k+\frac{j}{m}}}),\qquad j=1,\ldots,m-1.$

Znajdź wykładnik $p$ taki, że

$\| x_{{k+1}}-x^{*}\|\leq C\| x_{{k}}-x^{*}\|^{p}.$

Rozwiązanie:

Z definicji metody, $x_{{k+1}}$ wyznaczamy jako $m$ -tą iterację uproszczonej metody Newtona startującej z $x_{k}$ . Na mocy oszacowania błędu uproszczonej metody Newtona mamy więc

$\| x_{{k+\frac{j+1}{m}}}-x^{*}\|\leq C\| x_{{k}}-x^{*}\|\| x_{{k+\frac{j}{m}}}-x^{*}\|,$

skąd

$\| x_{{k+1}}-x^{*}\|\leq C\| x_{{k}}-x^{*}\|\| x_{{k+\frac{m-1}{m}}}-x^{*}\|\leq(C\| x_{{k}}-x^{*}\|)^{2}\| x_{{k+\frac{m-2}{m}}}-x^{*}\|\leq\cdots\leq C^{{m}}\| x_{{k}}-x^{*}\|^{{m+1}}.$

A więc metoda ma wykładnik $m+1$ .

Ćwiczenie 10.3

Porównaj efektywność metody Szamańskiego i Newtona w zależności od rozmiaru zadania $N$ , przy następujących modelowych założeniach:

koszt wyznaczenia $F(x)$ jest równy $c_{0}N$
koszt wyznaczenia $F^{{\prime}}(x)$ jest równy $c_{1}N^{2}$
koszt wyznaczenia rozkładu $F^{{\prime}}(x)$ jest równy $c_{r}N^{3}$
koszt wyznaczenia rozwiązania układu z macierzą $F^{{\prime}}(x)$ o danym rozkładzie jest równy $c_{s}N^{2}$

Rozwiązanie:

Oczywiście musimy założyć $m>1$ , w przeciwnym razie obie metody są tożsame. Dla metody Newtona mamy

$E_{{\text{Newton}}}=2^{{1/(c_{r}N^{3}+(c_{1}+c_{s})N^{2}+c_{0}N)}},$

a dla metody Szamańskiego

$E_{{\text{Szamański}}}=(m+1)^{{1/(c_{r}N^{3}+(c_{1}+mc_{s})N^{2}+mc_{0}N)}}.$

Stąd $E_{{\text{Szamański}}}\geq E_{{\text{Newton}}}$ wtedy i tylko wtedy, gdy $\log _{2}E_{{\text{Szamański}}}\geq\log _{2}E_{{\text{Newton}}}$ . Dalej przekształcając, dostajemy warunek wiążący ze sobą $m$ i $N$ :

$\log _{2}(m+1)\geq\dfrac{1+m\, O(\frac{1}{N})}{1+O(\frac{1}{N})},$

co zachodzi dla $N$ dostatecznie dużych względem $m$ .

Ćwiczenie 10.4 (liniowa niezmienniczość iteracji Newtona)

Rozważmy równanie $F(x)=0$ i jego liniową transformację

$G(y)=AF(B^{{-1}}y)=0,$

gdzie $A,B$ są nieosobliwymi macierzami rozmiaru $N\times N$ . Wykaż, że jeśli tylko $y_{0}=Bx_{0}$ , to iteracja Newtona dla $G$ ,

$y_{{k+1}}=y_{k}-G^{{\prime}}(y_{k})^{{-1}}G(y_{k})$

jest równoważna metodzie Newtona dla $F$ w tym sensie, że zawsze zachodzi $y_{k}=Bx_{k}$ .

Rozwiązanie:

Iteracja Newtona dla $F$ to

$x_{{k+1}}=x_{k}-F^{{\prime}}(x_{k})^{{-1}}F(x_{k}).$

Mnożąc lewostronnie przez $B$ i podstawiając gdzie trzeba $x_{k}=B^{{-1}}y_{k}$ mamy

$y_{{k+1}}=y_{k}-BF^{{\prime}}(B^{{-1}}y_{k})^{{-1}}F(B^{{-1}}y_{k})=y_{k}-BF^{{\prime}}(B^{{-1}}y_{k})^{{-1}}A^{{-1}}\, AF(B^{{-1}}y_{k}).$

Na zakończenie wystarczy zauważyć, że $BF^{{\prime}}(B^{{-1}}y_{k})^{{-1}}A^{{-1}}=(AF^{{\prime}}(B^{{-1}}y_{k})B^{{-1}})^{{-1}}=G^{{\prime}}(y_{k})^{{-1}}$ .

Ćwiczenie 10.5

Wykaż, że jeśli spełnione są założenia standardowe oraz metoda Newtona jest zbieżna, to dla dostatecznie dużych $k$ zachodzi

$\| F^{{\prime}}(x_{k})^{{-1}}F(x_{{k+1}})\|\leq\dfrac{1}{4}\| F^{{\prime}}(x_{k})^{{-1}}F(x_{{k}})\|.$

Może to więc być dość prosty test (niezmienniczy ze względu na liniowe transformacje zmiennej niezależnej i zależnej!), pozwalający przypuszczać, że nie zachodzi zbieżność: na przykład, gdy na pewnym kroku metody stwierdzimy $\| F^{{\prime}}(x_{k})^{{-1}}F(x_{{k+1}})\|\geq\dfrac{3}{4}\| F^{{\prime}}(x_{k})^{{-1}}F(x_{{k}})\|$ , możemy wtedy uznać prawdopodobny brak zbieżności i zakończyć iterację.

10.2.2. Wpływ niedokładności wyznaczenia pochodnej lub funkcji na zbieżność metody Newtona

Twierdzenie 10.3 (o lokalnej stabilności metody Newtona)

Przy standardowych założeniach, istnieją stałe dodatnie $C$ , $\delta$ oraz $\tilde{\delta}$ takie, że jeśli $x_{k}\in K(x^{*},\delta)$ oraz $\|\Delta _{k}\|<\tilde{\delta}$ , to

$\displaystyle x_{{k+1}}=x_{k}-\left(F^{{\prime}}(x_{k})+\Delta _{k}\right)^{{-1}}\left(F(x_{k})+\epsilon _{k}\right)$

(10.2)

jest dobrze określony oraz

$\displaystyle\| e_{{k+1}}\|\leq C\left(\| e_{k}\|^{2}+\|\Delta _{k}\|\| e_{k}\|+\|\epsilon _{k}\|\right),$

(10.3)

gdzie $e_{k}=x_{k}-x^{*}$ .

Dowód

Dowód znajdziemy np. w [9].

∎

Z twierdzenia o lokalnej stabilności metody Newtona można wywnioskować wiele twierdzeń o zbieżności tych modyfikacji metody Newtona, które sprowadzają się do zaburzenia oryginalnej metody Newtona, na przykład twierdzenia o zbieżności uproszczonej metody Newtona.

Dowód

(Twierdzenia 10.2, o zbieżności uproszczonej metody Newtona.) Jeden krok uproszczonej metody Newtona

$x_{{k+1}}=x_{k}-F^{{\prime}}(x_{0})^{{-1}}F(x_{k})$

możemy zapisać w języku twierdzenia o lokalnej stabilności metody Newtona, gdzie

$\Delta _{k}=F^{{\prime}}(x_{0})-F^{{\prime}}(x_{k}),\qquad\epsilon _{k}=0.$

Niech więc $\delta$ będzie dostatecznie małe tak, by dla $x_{k}$ zachodziło twierdzenie o lokalnej stabilności metody Newtona w kuli $K(x^{*},\delta)$ . Dodatkowo załóżmy, że $x_{0}\in K(x^{*},\delta)$ . Wtedy

$\|\Delta _{k}\|=\| F^{{\prime}}(x_{0})-F^{{\prime}}(x_{k})\|\leq L\| x_{0}-x_{k}\|\leq L\left(\| x_{0}-x^{*}\|+\| x_{k}-x^{*}\|\right)$

i na mocy właśnie tw. o lokalnej stabilności,

$\| e_{{k+1}}\|\leq\tilde{C}\left(\| e_{k}\|^{2}+L\left(\| e_{0}\|+\| e_{k}\|\right)\| e_{k}\|\right)\leq C\| e_{0}\|\| e_{k}\|.$

Ponieważ z założenia $\| e_{0}\|<\delta$ , to z powyższego wynika, że

$\| e_{{k+1}}\|\leq C\delta\| e_{k}\|$

jeśli więc $\delta$ jest na tyle małe, by dodatkowo $C\delta<1$ , to ciąg $x_{k}$ zawiera się w $K(x^{*},\delta)$ oraz $e_{k}$ musi być zbieżny do zera przynajmniej liniowo.

∎

10.2.3. Metoda z przybliżoną pochodną

Jak już wspominaliśmy, innym kłopotliwym momentem w realizacji metody Newtona jest wyznaczanie macierzy pochodnej. Opracowanie dokładnej formuły obliczania $F^{{\prime}}(x)$ może być żmudne, podatne na ludzkie pomyłki. Uzyskany wzór może być ciężki w implementacji i kosztowny w użyciu.

Ale, ponieważ metoda Newtona jest jedynie metodą przybliżoną, można zaryzykować użycie przybliżonej pochodnej — najlepiej: przybliżonej niskim kosztem. Jednym z pomysłów mogłoby być wyznaczenie przybliżenia pochodnej na podstawie ilorazów różnicowych, zgodnie z poniższym twierdzeniem.

Twierdzenie 10.4 (o różnicowej aproksymacji pochodnej kierunkowej)

Przyjmijmy, że są spełnione założenia standardowe. Niech $x$ będzie ustalonym punktem w $D$ i niech $w\in\mathbb{R}^{N}$ . Niech $h\in\mathbb{R}$ będzie na tyle małe, by $x+hw\in D$ . Wtedy

$\|\dfrac{F(x+hw)-F(x)}{h}-F^{{\prime}}(x)w\|\leq\dfrac{L}{2}\| w\|^{2}\cdot|h|.$

Twierdzenie to gwarantuje więc, że pochodną w kierunku wektora $w$ można aproksymować ilorazem różnicowym

$F^{{\prime}}(x)w\approx\dfrac{F(x+hw)-F(x)}{h}$

z błędem rzędu $O(h)$ .

Dowód

Twierdzenie jest prostą konsekwencją twierdzenia o wartości średniej i założenia lipschitzowskości pochodnej:

	$\displaystyle\dfrac{1}{h}\left(F(x+hw)-F(x)\right)$	$\displaystyle=\dfrac{1}{h}\int _{0}^{1}F^{{\prime}}(x+thw)\cdot hw\, dt=\int _{0}^{1}F^{{\prime}}(x)w+(F^{{\prime}}(x+thw)-F^{{\prime}}(x))w\, dt$
		$\displaystyle=F^{{\prime}}(x)w+\int _{0}^{1}(F^{{\prime}}(x+thw)-F^{{\prime}}(x))w\, dt.$

Stąd

$\|\dfrac{1}{h}\left(F(x+hw)-F(x)\right)-F^{{\prime}}(x)w\|\leq\int _{0}^{1}\| F^{{\prime}}(x+thw)-F^{{\prime}}(x)\|\| w\|\, dt\leq\int _{0}^{1}\| thw\|\| w\|\, dt=\dfrac{L}{2}\| w\|^{2}\cdot|h|.$

∎

Z powyższego wynika pomysł na przybliżenie całej macierzy pochodnej ilorazami różnicowymi: wystarczy wziąć $F^{{\prime}}(x_{k})\approx A_{k}$ , gdzie $j$ -ta kolumna $A_{k}$ , $(A_{k})_{j}$ , jest wyznaczona jako iloraz różnicowy

$(A)_{j}=\dfrac{1}{h}\left(F(x+he_{j})-F(x)\right)\approx F^{{\prime}}(x)e_{j}.$

Twierdzenie 10.5 (o zbieżności metody z pochodną przybliżoną ilorazami różnicowymi)

Przy standardowych założeniach, istnieją dodatnie $\delta$ i $h$ (dostatecznie małe) takie, że jeśli $(h_{k})$ jest ciągiem takim, że $0<|h_{k}|\leq h$ oraz $x_{0}\in K(x^{*},\delta)$ , to ciąg zdefiniowany wzorem

$x_{{k+1}}=x_{k}-A_{k}^{{-1}}F(x_{k}),$

w którym kolumny macierzy $A_{k}$ zadane są ilorazami różnicowymi:

$(A_{k})_{j}=\dfrac{1}{h}\left(F(x_{k}+h_{k}e_{j})-F(x_{k})\right),\qquad j=1,\ldots,N,$

jest dobrze określony i zbieżny przynajmniej liniowo do $x^{*}$ . Co więcej,

jeśli $h_{k}\rightarrow 0$ , to $x_{k}\rightarrow x^{*}$ superliniowo;
jeśli $|h_{k}|\leq M\| x_{k}-x^{*}\|$ dla pewnej stałej $M$ , to $x_{k}\rightarrow x^{*}$ przynajmniej kwadratowo;
jeśli $|h_{k}|\leq M\| F(x_{k}\|$ dla pewnej stałej $M$ , to $x_{k}\rightarrow x^{*}$ przynajmniej kwadratowo.

Dowód

Dla wygody, dowód poprowadzimy w normie $\|\cdot\| _{1}$ — tzw. kolumnowej normie macierzowej (w przestrzeni skończenie wymiarowej i tak wszystkie normy są równoważne). Mamy więc

$\| A_{k}\| _{1}=\max _{j}\|(A_{k})_{j}\| _{1},$

oraz

$\| A_{k}-F^{{\prime}}(x_{k})\| _{1}=\max _{j}\|(A_{k})_{j}-F^{{\prime}}(x_{k})e_{j}\| _{1}.$

Tymczasem, na mocy lematu o różnicowej aproksymacji pochodnej,

$\|(A_{k})_{j}-F^{{\prime}}(x_{k})e_{j}\| _{1}\leq\dfrac{L}{2}\| e_{j}\| _{1}^{2}|h_{k}|=\dfrac{L}{2}|h_{k}|.$

Oznaczmy $\Delta _{k}=A_{k}-F^{{\prime}}(x_{k})$ , wtedy $\|\Delta _{k}\| _{1}\leq\dfrac{L}{2}|h_{k}|$ . Na mocy twierdzenia o lokalnej stabilności metody Newtona z (dla $\delta,h$ dostatecznie małych) mamy, że $A_{k}$ jest nieosobliwa oraz zachodzi

$\| e_{{k+1}}\| _{1}\leq C\left(\| e_{k}\| _{1}^{2}+|h_{k}|\| e_{k}\| _{1}\right).$

Stąd już wynika teza twierdzenia. Ostatni warunek kwadratowej zbieżności wynika z oszacowania residuum przez błąd,

$\| F(x_{k})\| _{1}\leq 2\| F^{{\prime}}(x^{*})\| _{1}\| e_{k}\| _{1},$

gdy tylko $x_{k}$ jest dostatecznie blisko $x^{*}$ .

∎

10.2.4. Niedokładna metoda Newtona

Jak pamiętamy, najkosztowniejszą częścią jednej iteracji Newtona jest rozwiązywanie układu równań z macierzą pochodnej. Alternatywą dla metody bezpośredniej rozwiązywania równania poprawki

$F^{{\prime}}(x_{k})s=F(x_{k})$

mogłaby być, zwłaszcza w przypadku gdy $N$ jest duże, metoda iteracyjna (mielibyśmy zatem do czynienia z iteracją wewnątrz iteracji). Na $k$ -tym kroku metody Newtona moglibyśmy zatrzymywać wewnętrzną iterację stosując np. residualne kryterium stopu z parametrem wymuszającym $\eta _{k}$ ,

$\| F(x_{k})-F^{{\prime}}(x_{k})s\|\leq\eta _{k}\| F(x_{k})\|.$

Taką modyfikację metody Newtona nazywa się niedokładną metodą Newtona.

Twierdzenie 10.6 (o błędzie niedokładnej metody Newtona)

Przy standardowych założeniach, istnieją dodatnie stałe $C$ i $\delta$ takie, że jeśli $x_{k}\in K(x^{*},\delta)$ oraz $s$ spełnia warunek

$\| F^{{\prime}}(x_{k})s-F(x_{k})\|\leq\eta _{k}\| F(x_{k})\|,$

to następna iteracja niedokładnej metody Newtona dana wzorem

$x_{{k+1}}=x_{k}-s$

spełnia

$\| e_{{k+1}}\|\leq C\left(\| e_{k}\|+\eta _{k}\right)\| e_{k}\|.$

Dowód

Niech $\delta$ będzie na tyle małe, by zachodził lemat 9.2 o oszacowaniu funkcji i pochodnej, a także by zachodziło twierdzenie 9.3 o zbieżności metody Newtona z punktem startowym $x_{k}$ . Aby udowodnić twierdzenie, najpierw postaramy się wyłuskać związek pomiędzy naszą metodą, a prawdziwą metodą Newtona (o której już sporo wiemy). Mamy bowiem

$x_{{k+1}}=x_{k}-s=x_{k}-F^{{\prime}}(x_{k})^{{-1}}F(x_{k})+F^{{\prime}}(x_{k})^{{-1}}F(x_{k})-s=x_{{k+1}}^{{\text{Newton}}}+F^{{\prime}}(x_{k})^{{-1}}r,$

gdzie $r$ jest residuum rozwiązania równania na poprawkę,

$r=F(x_{k})-F^{{\prime}}(x_{k})s.$

Stąd wynika, że

$x_{{k+1}}-x^{*}=x_{{k+1}}^{{\text{Newton}}}-x^{*}+F^{{\prime}}(x_{k})^{{-1}}r$

i konsekwentnie

$\| x_{{k+1}}-x^{*}\|\leq\| x_{{k+1}}^{{\text{Newton}}}-x^{*}\|+\| F^{{\prime}}(x_{k})^{{-1}}r\|.$

Ponieważ na mocy twierdzenia o zbieżności metody Newtona mamy $\| x_{{k+1}}^{{\text{Newton}}}-x^{*}\|\leq\tilde{C}\| x_{{k}}-x^{*}\|^{2}$ , to wystarczy oszacować ostatni człon nierówności. Z definicji $r$ mamy

$\| F^{{\prime}}(x_{k})^{{-1}}r\|\leq\| F^{{\prime}}(x_{k})^{{-1}}\|\| F(x_{k})-F^{{\prime}}(x_{k})s\|\leq 2\| F^{{\prime}}(x^{*})^{{-1}}\|\cdot\eta _{k}\| F(x_{k})\|\leq 4\operatorname{cond}(F^{{\prime}}(x^{*}))\eta _{k}\| x_{k}-x^{*}\|.$

(W ostatniej nierówności skorzystaliśmy z lematu o oszacowaniu funkcji i pochodnej.)

∎

W powyższym dowodzie dostaliśmy rezultat, w którym stała $C$ w oszacowaniu błędu zależała od uwarunkowania $F^{{\prime}}(x^{*})$ , skąd moglibyśmy wysnuć wniosek, że jeśli $F^{{\prime}}(x^{*})$ jest źle uwarunkowana, to $\eta _{k}$ należy brać patologicznie małe. W rzeczywistości nie jest aż tak źle i można osłabić założenia na ciąg $(\eta _{k})$ , ale pod warunkiem zmiany normy, w której badamy błąd, tak, by uwzględniać zachowanie się pochodnej: $\|\cdot\|=\| F^{{\prime}}(x^{*})\cdot\|$ (zob. [9]).

Wniosek 10.1 (twierdzenie o zbieżności niedokładnej metody Newtona)

Przy standardowych założeniach, istnieją $\delta$ i $\eta$ (dostatecznie małe) takie, że jeśli $x_{0}\in K(x^{*},\delta)$ oraz $0\leq\eta _{k}\leq\eta$ to niedokładna metoda Newtona z parametrem wymuszającym $\eta _{k}$ jest zbieżna przynajmniej liniowo do $x^{*}$ . Ponadto,

jeśli $\eta _{k}\rightarrow 0$ , to zbieżność jest superliniowa;
jeśli $\eta _{k}\leq M_{\eta}\| F(x_{k})\|$ dla pewnego ustalonego $M_{\eta}$ , to zbieżność jest przynajmniej kwadratowa.

Ćwiczenie 10.6

Udowodnij powyższy wniosek.

Niedokładna metoda Newtona jest bardzo popularną metodą w przypadku, gdy $N$ jest bardzo duże: wówczas trudno myśleć o sposobach rozwiązywania równania poprawki innych niż jakaś metoda iteracyjna. Dodatkowym plusem tej metody jest to, że znakomicie można ją łączyć z przybliżaniem pochodnej kierunkowej ilorazami różnicowymi; rzeczywiście, głównym składnikiem metody iteracyjnej jest mnożenie $F^{{\prime}}(x_{k})w$ , co można tanio przybliżać w sposób opisany w rozdziale 10.2.3. W ten sposób oszczędzamy na kilku polach na raz:

nie musimy wykonywać, zazwyczaj bardzo kosztownego, rozkładu macierzy pochodnej;
nie musimy nawet wyznaczać macierzy pochodnej, lecz w zamian wystarczy, że w każdej wewnętrznej iteracji, jedną dodatkową wartość $F(x_{k}+hw)$ ;
ponieważ po kilku iteracjach zbliżamy się do rozwiązania, że $x_{{k+1}}\approx x_{{k}}$ , w ten sposób od razu dysponujemy dobrym przybliżeniem startowym dla iteracji wewnętrznej: w sprzyjających okolicznościach metody Kryłowa mogą z tego zrobić dodatkowy pożytek;
możemy ograniczyć koszt rozwiązywania układu z macierzą pochodnej (liczbę iteracji wewnętrznych), stosując łagodne kryterium stopu, gdy jesteśmy daleko od rozwiązania.

Przykład 10.1 (Numeryczne eksperymenty z wariantami metody Newtona dla równania Allena–Cahna)

Porównamy jakość pracy kilku metod:

klasycznej metody Newtona
uproszczonej metody Newtona
metody Szamańskiego o $m=4$ krokach
metody Newtona z pochodną przybliżoną ilorazami różnicowymi (ze stałym, ale bardzo małym krokiem, $h\approx\sqrt{\epsilon}$ , gdzie $\epsilon$ to precyzja arytmetyki.

Do porównania wybierzemy:

wykresy uzyskanego rozwiązania, $u(x)$ , dla każdej z metod
wykresy residuum, $|F(u(x))|$ , dla każdej z metod
historię zbieżności metody, $\| F(u_{j})\| _{2}$
liczbę iteracji, czas pracy, końcową wartość residuum.

%
% rozwiazania zadania z macierza rozrzedzona
% 
disp('Matematyka obliczeniowa II');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x, resid, info, output] = newton(nazwa_f, nazwa_df, x0, atol, rtol, step, maxit)
%
%	[x, resid, info, output] = newton(nazwa_f, x0, atol, rtol, step, maxit)
%
% Rozwiazuje uklad rownan nieliniowych F(x) = 0, gdzie F jest zadana M-funkcja o 
% nazwie ”nazwa_f”.
% Parametry wejscia:
% 	nazwa_f - nazwa M-funkcji zadajacej funkcje F(). 
% 		Ta funkcja *musi* przyjmowac takze argumenty macierzowe!
% 	x0 - przyblizenie poczatkowe rozwiazania
%	atol, rtol - tolerancja bledu residualnego
%	step - co ile iteracji uaktualniac macierz pochodnej; step = 1 daje standardowa metode Newtona, step > maxit daje metode cieciw
%	maxit - maksymalna dopuszczalna liczba iteracji
% Wyjscie:
%	x - koncowe przyblizenie rozwiazania
%	resid - wektor zawierajacy historie zbieznosci, resid(i) - wartosc
% 		normy residuum w i-tej iteracji
%	info - kod zakonczenia: = 1 - sukces, = 0 - przekroczona maxit
%	output - struktura zawierajaca statystyki

% kilka pomocniczych definicji
epsq = sqrt(eps);	% pierwiastek z podwojonej precyzji maszyny
unit = ones(size(x0));	% wektor zlozony z samych jedynek

% poczatkowe rozwiazanie przyblizone
x = x0;
	
% wyznaczenie residuum poczatkowego
fx = feval(nazwa_f, x);
iter = 1;
resid(iter) = norm(fx);

while( (resid(iter) >= (atol + rtol*resid(1))) && (iter < maxit) )

if (mod(iter-1, step) == 0)
		% aproksymacja pochodnej - wszystkie kolumny jednoczesnie
		h = epsq*max(abs(x),unit);
		xh = x*unit' + diag(h);
		if exist(nazwa_df) > 0 
			Df = feval(nazwa_df, x);
		else
			Df = (feval(nazwa_f, xh) - fx*unit')./(unit*h');
		end
		% rozklad macierzy pochodnej do wykorzystania
		[L, U, p] = lu(Df, 'vector');
	end
	% poprawka przyblizonego rozwiazania, x = x - Df\fx;
	x = x - U\(L\fx(p));
	
	% wyznaczenie residuum
	fx = feval(nazwa_f, x);
	iter = iter + 1;
	resid(iter) = norm(fx);
end
info = ~(iter == maxit);
output.iterations = iter - 1;
end

function y = acfun(u) % lewa strona rownania Allena-Cahna
global T; % macierz laplasjanu
global delta;
global x; % wektor wezlow zmiennej przestrzennej
y = T*u + delta*u.*(1-u.^2) - repmat(delta*sqrt(x),[1, size(u,2)]);
end

function y = acdfun(u) % pochodna acfun
global T; % macierz laplasjanu
global delta;
y = T + delta*spdiag(1-3*u.^2);
end

global T;
global delta;
global x;

atol = 1e-8; rtol = 0; maxit = 300; delta = 10;
Nrange = [250]; % zestaw wartosci testowych rozdzielczosci dyskretyzacji
for k = 1:length(Nrange)
	N = Nrange(k); % rozmiar macierzy
	x = linspace(0,1,N+2)'; x = x(2:end-1);
	e = ((N+1)^2)*ones(N,1); T = spdiags([-e, 2*e, -e], [-1,0,1], N,N);
	u0 = ones(N,1);

% metoda Newtona
	tic; [un, residn] = newton('acfun', 'acdfun', u0, atol, rtol, 1, maxit);
	timen = toc;
	% historia zbieznosci metody 
	semilogy(residn,'b-o;|F(u_{Newton})|;'); pause(1);

% uproszczona metoda Newtona
	tic; [uun, residun] = newton('acfun', 'acdfun', u0, atol, rtol, maxit+1, maxit);
	timeun = toc;
	% historia zbieznosci metody 
	semilogy(residun,'r-o;|F(u_{uproszcz Newton})|;'); pause(1);
	
	% metoda Szamanskiego z odswiezeniem co 4 iteracje
	tic; [usz, residsz] = newton('acfun', 'acdfun', u0, atol, rtol, 4, maxit);
	timesz = toc;
	% historia zbieznosci metody 
	semilogy(residsz,'g-o;|F(u_{Szamanski})|;'); pause(1);
	
	% metoda Newtona z przyblizona pochodna
	tic; [unfd, residnfd] = newton('acfun', ”, u0, atol, rtol, 1, maxit);
	timenfd = toc;
	% historia zbieznosci metody 
	semilogy(residnfd,'m-o;|F(u_{Newton przybl poch})|;'); pause(1);

% rozwiazanie
	plot(x,[un,uun,usz,unfd], ['b-;u_{Newton};';'r-;u_{Newton upr};';'g-;u_{Szamanski};';'m-;u_{Newton przybl poch};']); pause(1);
	
	% residuum
	semilogy(x,abs(acfun([un,uun,usz,unfd])), ['b-;|F(u_{Newton})|;';'r-;|F(u_{Newton upr})|;';'g-;|F(u_{Szamanski})|;';'m-;|F(u_{Newton przybl poch})|;']); pause(1);
	
	
	[length(residn), length(residun), length(residsz), length(residnfd)]
	[timen,timeun,timesz,timenfd]
	[residn(end), residun(end), residsz(end), residnfd(end)]
end

Zwróć uwagę na to, że odpowiednio dobierając częstość uaktualniania macierzy pochodnej, można znacząco poprawić efektywność metody (w porównaniu do metody Newtona). Z drugiej strony, zbyt rzadka aktualizacja może przeszkodzić w uzyskaniu zbieżności, jak to dzieje się w naszym przykładzie w przypadku uproszczonej metody Newtona.

Ćwiczenie 10.7

Sprawdź, jak zmienią się wyniki (uzyskane rozwiązanie, jego dokładność, a także — koszt metody i szybkość jej zbieżności), gdy zmienisz jeden z poniższych parametrów

osłabisz nieliniowość, biorąc $\delta=1$ ,
nieco wzmocnisz nieliniowość, biorąc $\delta=15$ ,
zwiększysz rozmiar zadania $N$ do 400,
zmniejszysz $N$ do 20,
zmniejszysz tolerancję błędu do $10^{{-12}}$ .

Wskazówka:

Zdaje się, że dla $N=20$ mamy niejednoznaczność rozwiązania? Jak to zweryfikować?

Ćwiczenie 10.8

Zaimplementuj niedokładną metodę Newtona i wykorzystaj ją w skrypcie rozwiązującym (najlepiej: dwuwymiarowe) równanie Allena–Cahna. Jak wpłynęło to na efektywność metody, w porównaniu do standardowej metody Newtona?

Przykład 10.2

Monografii Kelley'a towarzyszy zestaw skryptów, pozwalających przetestować działanie różnych metod iteracyjnych w kilku praktycznych zadaniach [9].

Poniższy skrypt, którego podstawowe kody źródłowe pochodzą ze strony http://www4.ncsu.edu/~ctk/newton, umożliwi Ci zapoznanie się z zadaniem rozwiązania równania Chandrasekhara [9, rozdział 5.6].

Aby w pełni wykorzystać możliwości skryptu, zachęcamy do uruchomienia go na własnym (podłączonym do internetu) komputerze i następnie do wnikliwej obserwacji zmiany wyników przy zmianie, czy to parametrów zadania, czy to parametrów pracy solvera.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.