\[ \def\_#1{\mathbf{#1}} \def\x{\times} \def\R{\mathbb{R}} \def\mat#1{\begin{bmatrix}#1\end{bmatrix}} \def\matr#1{\begin{bmatrix*}[r]#1\end{bmatrix*}} \]

NLLS: Kružnice

Je dáno $m$ bodů v rovině, $\_a_1,\ldots,\_a_m\in\R^2$. Hledáme kružnici se středem $\_c\in\R^2$ a poloměrem $r\ge 0$ takovou, že součet čtverců vzdáleností bodů od kružnice je nejmenší.

Označme jako $\_x=(\_c,r)\in\R^3$ vektor parametrů kružnice. Nechť ${\rm dist}(\_x, \_a)$ je orientovaná vzdálenost bodu $\_a$ od kružnice s parametry $\_x$. Tedy $\lvert{\rm dist}(\_x, \_a)\rvert$ je vzdálenost bodu $\_a$ od kružnice, přičemž pro $\_a$ vně kružnice je ${\rm dist}(\_x,\_a)>0$ a pro $\_a$ uvnitř kružnice je ${\rm dist}(\_x,\_a)<0$. Chceme minimalizovat funkci \begin{equation} f(\_x) = \sum_{i=1}^m {\rm dist}(\_x, \_a_i)^2 \label{eq:obj} \end{equation}

Implementační úkoly

  1. Napište funkci d = dist(x,A), kde x je vektor $3\times 1$, A je matice $2\times m$ obsahující body $\_a_1, \_a_2, \ldots, \_a_m$ a d je vektor $m\times 1$ takový, že d(i) je orientovaná vzdálenost bodu $\_a_i$ od kružnice s parametry $\_x$.
  2. Implementujte funkci x_new = make_GN_iter(x,A), která provede jednu iteraci čisté (tedy s jednotkovou délkou kroku) Gauss-Newtonovy metody.
  3. Implementujte funkci [x_new,success] = make_LM_iter(x,A,mu), která provede jednu iteraci Levenberg-Marquardtovy metody (bližší popis funkcí je v jejich šablonách.)

Poznámky:

  • Pro práci na úloze máte k dispozici skript main.m a funkci fit_circle.m.
  • Doporučujeme si napsat jednoduchý skript, který vám dovolí naklikat body $\_a_i$ a uložit si je do souboru .mat (bude se hodit funkce ginput). Takto si připravte několik vhodných množin bodů, na kterých můžete zkoušet svůj kód.

Úkoly do zprávy

  1. Pro $m\ge3$ a body $\_a_1, \ldots, \_a_m$ v obecné poloze, zkuste nalézt co nejobecnější podmínku charakterizující, ve kterých bodech $\_x$ je funkce $f$ diferencovatelná. Odpověď důvodněte.
  2. Může se zdát, že iterační metody na nelineární nejmenší čtverce bez omezení (např. Gauss-Newtonovu a Levenberg-Marquardtovu metodu) nejde použít, protože máme omezení $r\ge0$. Vadí to? Co se stane, budeme-li toto omezení ignorovat? Můžou metody konvergovat k řešení se záporným $r$? Odpovědi odůvodněte.
  3. Může mít pro nějaké body $\_a_1,\ldots,\_a_m$ funkce $f$ více lokálních minim s různou funkční hodnotou? Pokud odpovíte záporně, odůvodněte. Pokud odpovíte kladně, najděte příklad (body $\_a_1,\ldots,\_a_m$ a dvě lokální minima $\_x^1,\_x^2$ funkce $f$ s různou funkční hodnotou) a odůvodněte, proč $\_x^1,\_x^2$ jsou lokální minima (můžete použít argument, že když např. Gauss-Newtonova metoda zkonverguje do stacionárního bodu, tak tento bod je téměř jistě lokální minimum a není tedy třeba ověřovat definitnost Hessiánu). Do zprávy pak
    • exportujte (pomocí matlabského příkazu print -depsc obr.eps nebo print -dsvg obr.svg) obrázek, ve kterém budou body $\_a_1,\ldots,\_a_m$ a dvě kružnice s parametry $\_x^1$ a $\_x^2$,
    • napište funkční hodnoty $f(\_x^1)$ a $f(\_x^2)$,
    • napište gradienty $\nabla f(\_x^1), \nabla f(\_x^2)$ (pokud by funkce v bodech $\_x^1$ příp. $\_x^2$ nebyla diferencovatelná, diskutujte a zkuste najít jiný argument, proč jsou body lokální minima).
  4. (nepovinný úkol pro nadšence) Zatím jsme prokládání kružnice formulovali jako minimalizaci funkce \eqref{eq:obj} přes tři proměnné $\_x=(\_c,r)$, tedy $$\min_{\_c\in\R^2,\,r\ge0} f(\_c,r) = \min_{\_c\in\R^2} \overbrace{\min_{r\ge0} f(\_c,r)}^{g(\_c)}.$$ Ale pokud je $\_c$ pevné, vnitřní minimalizaci přes $r$ lze vyřešit explicitně a napsat tak vzorec pro funkci $g(\_c)$. Tedy úlohu můžeme vlastně formulovat jako minimalizaci funkce $g$ přes proměnné $\_c\in\R^2$. Rozpracujte tuto myšlenku: odvoďte vzorec pro funkci $g$, implementujte iterační metodu (Gauss-Newtona a/nebo Levenberg-Marquardta) pro minimalizaci funkce $g$ a diskutujte, která formulace (dvě vs. tři proměnné) je výhodnější (např. porovnáním rychlosti/přesnosti iteračních metod příp. pracnosti implementace).

Postup prací:

  • Z gitlabu si stáhněte šablony pro funkce k implementaci a pomocné funkce a skripty.
  • Implementujte požadované matlabské funkce.
  • Napište PDF zprávu a pojmenujte ji report.pdf
  • Zabalte všechny soubory .m a soubor report.pdf do ZIP souboru a nahrajte je do Brute. Udělejte ZIP soubor tak, aby se vaše soubory rozbalily rovnou do aktuálního adresáře, ne do nějakého podadresáře (jinak to nebude fungovat.)
courses/b0b33opt/cviceni/hw/kruznice/start.txt · Last modified: 2020/04/20 00:48 by wernetom