Numerické metody v R

$$ \xdef\scal#1#2{\langle #1, #2 \rangle} \xdef\norm#1{\left\lVert #1 \right\rVert} \xdef\dist{\rho} \xdef\and{\&}\xdef\brackets#1{\left\{ #1 \right\}} \xdef\parc#1#2{\frac {\partial #1}{\partial #2}} \xdef\mtr#1{\begin{pmatrix}#1\end{pmatrix}} \xdef\bm#1{\boldsymbol{#1}} \xdef\mcal#1{\mathcal{#1}} \xdef\vv#1{\mathbf{#1}}\xdef\vvp#1{\pmb{#1}} \xdef\ve{\varepsilon} \xdef\l{\lambda} \xdef\th{\vartheta} \xdef\a{\alpha} \xdef\tagged#1{(\text{#1})} \xdef\tagged*#1{\text{#1}} \xdef\tagEqHere#1#2{\href{#2\#eq-#1}{(\text{#1})}} \xdef\tagDeHere#1#2{\href{#2\#de-#1}{\text{#1}}} \xdef\tagEq#1{\href{\#eq-#1}{(\text{#1})}} \xdef\tagDe#1{\href{\#de-#1}{\text{#1}}} \xdef\T#1{\htmlId{eq-#1}{#1}} \xdef\D#1{\htmlId{de-#1}{\vv{#1}}} \xdef\conv#1{\mathrm{conv}\, #1} \xdef\cone#1{\mathrm{cone}\, #1} \xdef\aff#1{\mathrm{aff}\, #1} \xdef\lin#1{\mathrm{Lin}\, #1} \xdef\span#1{\mathrm{span}\, #1} \xdef\O{\mathcal O} \xdef\ri#1{\mathrm{ri}\, #1} \xdef\rd#1{\mathrm{r}\partial\, #1} \xdef\interior#1{\mathrm{int}\, #1} \xdef\proj{\Pi} \xdef\epi#1{\mathrm{epi}\, #1} \xdef\grad#1{\mathrm{grad}\, #1} \xdef\hess#1{\nabla^2\, #1} \xdef\subdif#1{\partial #1} \xdef\co#1{\mathrm{co}\, #1} $$

Rychlost konvergence

Definice $\D{3.1}$

Nechť jsou dány 2 posloupnosti $\brackets{e_ k}_ {k = 0}^\infty$ a $\brackets{h_ k}_ {k = 0}^\infty$ takové, že $$ e_k \in [0, \infty), \quad e_k \to 0 \quad \and \quad h_k \in [0, \infty), \quad h_k \to 0. $$ Řekneme, že posloupnost $\brackets{e_k}$ konverguje rychleji (pomaleji) než $\brackets{h_k}$, pokud existuje index $\tilde k \in \N_0$ takový, že $$ e_k \leq_{(\geq)} h_k \quad \forall k \in [\tilde k, \infty) \cap \N_0 $$

Definice $\D{3.2}$ (Rychlost konvergence)

Nechť je dána posloupnost $\brackets{e_ k}_ {k = 0}^\infty$ splňující $e_k \in [0, \infty)$ a $e_k \to 0$. Řekneme, že posloupnost $\brackets{e_k}$ konverguje

• alespoň lineárně s rychlostí $\beta \in (0,1)$, pokud konverguje rychleji než geometrická posloupnost se členy tvaru $q \bar \beta^k$, kde $q > 0$ a $\bar \beta \in (\beta, 1)$.
  • čím větší $\bar \beta$, tím pomaleji jde tato geometrická posloupnost k nule - tj. konverguje rychleji než geometrická posloupnost s $\bar \beta$ větší než $\beta$
• nejvýše lineárně s rychlostí $\beta \in (0,1)$, pokud konverguje pomaleji než geometrická posloupnost se členy tvaru $q \bar \beta^k$, kde $q > 0$ a $\bar \beta \in (0, \beta)$.
• lineárně s rychlostí $\beta \in (0,1)$, pokud konverguje nejvýše a současně alespoň lineárně s rychlostí $\beta$.
• superlineárně (sublineárně), pokud konverguje rycheji (pomaleji) než libovolná geometrická posloupnost se členy tvaru $q \beta^k$, kde $q > 0$ a $\beta \in (0,1)$.

Definice $\D{3.4}$

Nechť je dána posloupnost $\brackets{e_ k}_ {k = 0}^\infty$ splňující $e_k \in [0, \infty)$ a $e_k \to 0$, přičemž $\brackets{e_k}$ konverguje superlineárně. Řekneme, že posloupnost $\brackets{e_k}$ konverguje

• alespoň superlineárně s řádem $p > 1$, pokud konverguje rychleji než všechny posloupnosti se členy tvaru $q \beta^{\bar p^k}$, kde $q > 0$ a $\beta \in (0,1)$ a $\bar p \in (1, p)$
  • čím větší $\bar p$, tím rychleji posloupnost $q \beta^{\bar p^k}$ konverguje - tj. $\brackets{e_k}$ konverguje rychleji než všechny posloupnosti s menším $p$
• nejvýše superlineárně s řádem $p > 1$, pokud konverguje pomaleji než všechny posloupnosti se členy tvaru $q \beta^{\bar p^k}$, kde $q > 0$ a $\beta \in (0,1)$ a $\bar p \in (p, \infty)$
• superlineárně s řádem $p > 1$, pokud konverguje nejvýše a současně alespoň superlineárně s řádem $p$.
• superlineárně s řádem $p = 1$, pokud konverguje pomaleji než všechny posloupnosti se členy tvaru $q \beta^{\bar p^k}$, kde $q > 0$ a $\beta \in (0,1)$ a $\bar p \in (1, \infty)$

Metody

Definice $\D{3.1.1}$ (Unimodální funkce)

Nechť je dán interval $I \subset \R$ a funkce $f: I \to \R$. Řekneme, že $f$ je unimodální na $I$, jestliže existuje $x^* \in I$ takové, že

• $f(x_1) > f(x_2)$ pro libovolná $x_1, x_2 \in I$ splňující $x^* > x_1 > x_2$
• $f(x_1) < f(x_2)$ pro libovolná $x_1, x_2 \in I$ splňující $x^* < x_1 < x_2$

Jinými slovy, unimodální funkce je klesající na $(-\infty, x^* ) \cap I$ (tj. nalevo od $x^* $) a rostoucí na $I \cap (x^* , \infty)$ (tj. napravo od $x^* $).

Unimodalita neimplikuje konvexnost (ani se spojitostí), pouze kvazikonvexnost

Naopak, konvexní funkce nemusí nutně být unimodální (ale ostrá konvexnost $\implies$ unimodalita)

Konvexní funkce nemusí být např. jen rostoucí, ale i neklesající

V této části budeme řešit úlohu $$ f(x) \to \min, \qquad x \in I := [a,b] \tag{\T{3.1.1}} $$

Lemma $\D{3.1.2}$

Nechť $f: I \to \R$ je unimodální na $I$ a $x_1, x_2 \in I$ jsou takové, že $x_1 < x_2$.

• Je-li $f(x_1) \leq f(x_2)$, pak $x^* \leq x_2$
• Je-li $f(x_1) \geq f(x_2)$, pak $x^* \geq x_1$

Upozornění:
Dále uvažujme POUZE UNIMODÁLNÍ funkce.

Navíc nechť $N$ značí povolený počet vyčíslení a přesnost těchto metod je dáno jako $|\bar x - x^* |$, kde $x^* $ je přesné řešení úlohy $\tagEq{3.1.1}$ a $\bar x$ jeho nalezená aproximace.

Metoda prostého dělení

Tato metoda není efektivní a je to de facto hrubá síla

Podle parity $N$ určíme dělící body intervalu $I$.

kde $\delta$ je vhodné malé číslo.

Poté vyčíslíme $f(x_1), \dots, f(x_N)$ (což v případě $N$ sudého a $\delta \in \set{0, {b-a \over k+ 1}}$ znamená pouze $k$ vyčíslení) a nechť v $x_j$ nastává nejmenší hodnota, tj. $$ f(x_j) = \min_{1 \leq i \leq N} f(x_i) $$ Pak z Lemma $\tagDe{3.1.2}$ plyne, že $x^* \in [x_{j-1}, x_{j+1}]$ a tento interval nazveme interval lokalizace minima (ILM) a za aproximaci $x^* $ vezmeme střed ILM, tj. $\bar x := {x_{j-1} + x_{j+1} \over 2}$

Pro délku $l_N$ intervalu lokalizace minima platí $$ l_N := \max_{1 \leq i \leq N} (x_{i+1} - x_{i-1}) = \begin{cases} 2 {b - a \over N+1}, & N = 2k - 1 \ {b - a \over (N/2) + 1} + \delta, & N = 2k \end{cases} $$

Pro $N$ sudé je poslední interval delší, proto dostáváme takový tvar $l_N$

Přesnost této metody je dána polovinou ILM, tj. ${l_N \over 2}$

Rychlost konvergence této metody je sublineární, navíc je tento algoritmus pasivní, tj. volba $x_{m+1}$ nezáleží na $x_1, \dots, x_m$ (závisí pouze na $N$, či na $N$ a $\delta$).

Při rovnosti funkčních hodnot preferujeme konec (teoreticky by to mělo být jedno)

Metoda půlení intervalu

Nechť nyní $N = 2k$. Položme $a_0 = a$, $b_0 = b$ a $$ x_i^- := {a_{i-1} + b_{i-1} \over 2} - \delta \quad \and \quad x_i^+ := {a_{i-1}+b_{i-1} \over 2} + \delta, $$ kde $\delta > 0$ je dostatečně malé a $i = 1, \dots, k$. Vyčíslíme funkci v $x_i^-, x_i^+$, tj. dostaneme $f(x_i^-), f(x_i^+)$. Pak

• jestliže $f(x_i^-) < f(x_i^+)$, pak podle Lemma $\tagDe{3.1.2}$ je ILM $[a_{i-1}, x_i^+] \implies a_i = a_{i-1}, b_i = x_i^+$
• jestliže $f(x_i^-) > f(x_i^+)$, pak podle Lemma $\tagDe{3.1.2}$ je ILM $[x_i^-, b_{i-1}] \implies a_i = x_i^-, b_i = b_{i-1}$

Takto můžeme tento proces opakovat ($k$-krát, jelikož máme $N = 2k$ povolených vyčíslení), kdy za $a,b$ volíme krajní body ILM pro každý krok. Zřejmě, jako aproximaci $x^* $ v $k$-tém kroku bereme střed ILM pro $k$-tý krok.

Délka ILM je v tomto případě $$ l_k = {b - a\over 2^k} + {(2^k - 1)\delta \over 2^{k-1}}, $$ přičemž $\delta \in [0, {b-a\over 2}]$ a navíc pro $k \to \infty$ je $l_k \to 2\delta$.

Z tohoto vyplývá, že čím menší $\delta$, tím je metoda přesnější. Nicméně ve skutečnosti se můžeme dostat k zaokrouhlovacím chybám, které dokonce mohou způsobit, že špatně určíme velikosti $f(x_i^-), f(x_i^+)$ (tím pádem bychom řekli, že $x^* $ je v opačném intervalu, než je ve skutečnosti)

Tato metoda konverguje lineárně s rychlostí $1/2$.

Při rovnosti funkčních hodnot zapomínáme konec (teoreticky by to mělo být jedno)

Metoda zlatého řezu

Myšlenka metody zlatého řezu "vylepšuje" metodu půlení intervalu tak, že každá další iterace umožňuje pouze jedno další vyčíslení.

Zde $\tau$ je řešení rovnice $\tau^2 - \tau - 1 = 0$, tj. a $\frac 1 \tau \approx 0.618$

Mějme funkci $f$, interval $[a,b]$, přesnost $\ve$ nebo počet vyčíslení $N \geq 2$:

1. (Inicializace) Položíme $a_0 := a, b_0 := b$ a $k := 1$. Vypočteme
   $$\l_1 := a_0 + {b_0 - a_0 \over \tau^2} \quad \and \quad \mu_1 := a_0 + {b_0 - a_0 \over \tau}$$
2. Je-li $k = N$, pokračujeme částí 5., jinak následuje krok 3.
3. Vyčíslíme $f(\l_k)$ a $f(\mu_k)$. Jestliže $f(\l_k) \geq f(\mu_k)$:
   1. Položíme $a_k := \l_k, b_k = b_{k-1}, \l_{k+1} := \mu_k$ a
      $$f(\l_{k+1}) := f(\mu_k), \quad \mu_{k+1} := a_k + {b_k - a_k \over \tau}$$ a pokračujeme na krok 4.
   2. Položíme $a_k := a_{k-1}, b_k := \mu_k, \mu_{k+1} := \l_k$ a
      $$f(\mu_{k+1}) := f(\l_k), \quad \l_{k+1} := a_k + {b_k - a_k \over \tau^2}$$ a pokračujeme na krok 4.
4. Položíme $k := k+1$ a pokračujeme krokem 2.
5. Stanovíme poslední ILM jako $[a_{k-1}, b_{k-1}]$ a vypočteme $\bar x := {a_{k-1} + b_{k-1} \over 2}$. KONEC

Tato metoda konverguje lineárně s rychlostí $\frac 1 \tau \approx 0.618$

Toto neznamená, že by na stejný počet vyčíslení byla tato metoda horší než metoda půlení intervalu

Fibonacciho metoda

V této poslední metodě uvažujme, že zkrácení $\delta$ může být jiné v každé kroku metody.

Nechť $F_n$ je $n$-té Fibonacciho číslo

Máme povoleno $N$ vyčíslení, takže $M = N - 1$ a $$ \l_i = a_{i - 1} + {F_{N - i -1} \over F_{N - i + 1}} l_{i-1} = b_{i-1} - {F_{N - 1} \over F_{N - 1 + i}} l_{i-1} $$ $$ \mu_i = a_{i-1} + {F_{N - 1} \over F_{N - 1 + i}} l_{i-1} = b_{i-1} - {F_{N - i -1} \over F_{N - i + 1}} l_{i-1} $$

Tato metoda konverguje lineárně s rychlostí $\frac 1 \tau \approx 0.618$, tj. stejně jako metoda zlatého řezu

Fibonacciho metoda je (mírně) přesnější, než metoda zlatého řezu (která lze vnímat jako limitní varianta Fibonacciho metody). Nicméně u Fibonacciho metody je při změně $N$ potřeba všechny body přepočítat, což u metody zlatého řezu není.