\[ \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*}} \def\dist{\mathop{\text{dist}}} \]
Hledáme kružnici, která nejlépe prokládá dané body $(x_1,y_1),\ldots,(x_n,y_n)\in\R^2$ v tom smyslu, že minimalizuje součet čtverců vzdáleností těchto bodů od kružnice. Tedy minimalizujeme funkci \begin{equation}\label{objective} f(x_0,y_0,r) = \sum_{i=1}^n \dist(x_i,y_i,x_0,y_0,r)^2, \end{equation} kde $\dist(x,y,x_0,y_0,r)$ označuje vzdálenost bodu $(x,y)$ od kružnice se středem $(x_0,y_0)$ a poloměrem $r$ (vymyslete, jak se tato vzdálenost spočítá!).
Najít globální minimum funkce \eqref{objective} je obecně velmi obtížné. Budeme proto hledat pouze její lokální minimum, a to pomocí iteračních metod. Použijeme gradientní metodu a protože je $f$ ve tvaru nelineárních nejmenších čtverců, můžeme použít i Gauss-Newtonovu a Levenberg-Marquardtovu metodu.
Napište funkce vykonávajíní jednu iteraci gradientní metody, Gauss-Newtonovy metody a Levenberg-Marquardtovy metody na minimalizaci funkce $f$. U gradientní metody musíme volit délku kroku. K tomu lze použít line-search – to ale dělat nebudete a délku kroku dostane vaše funkce jako argument. Vyzkoušejte si sami, jaký vliv má tato délka kroku na konvergenci metody.
Jako úkoly napište tyto funkce:
d = dist(X, x0, y0, r)
, kde řádky matice $\_X\in\R^{n\times 2}$ jsou body $(x_i, y_i)$ a $\_d\in\R^{1\times n}$ je řádkový vektor orientovaných (uvnitř kružnice záporné, vně kladné) vzdáleností bodů od kružnice se středem $(x_0, y_0)$ a poloměrem $r$. Tuto funkci použijete při implementai následujících funkcí.
[x y r] = grad_iter(X, x0, y0, r0, a)
, kde (x0,y0,r0)
jsou parametry kružnice před iterací, (x,y,r)
jsou parametry po vykonání iterace a a
je délka kroku gradientní metody.
[x y r] = GN_iter(X, x0, y0, r0)
[x y r success] = LM_iter(X, x0, y0, r0, mu)
, kde mu
je parameter LM-metody a success
indikuje jestli došlo ke zlepšení účelové funkce touto iterací.
f = get_objective_function(X)
, kde f
je funkce (případně handle na funkci), která přijímá vektor parametrů c=[x0,y0,r]
a hodnota f( c )
odpovídá hodnotě $f(x_0,y_0,r)$ pro body $X$, jak jsme ji definovali v \eqref{objective}. Můžeme tedy například provést volání f([1,2,3])
a bude platit f([1,2,3]) == sum(dist(X,1,2,3).^2)
, případně f([1,2,3]) == sum(dist(X,1,2,3)**2)
.
Templaty, včetně skriptu pro testování si stáhněte zde: template pro matlab, template pro python. Pro python spouštějte skript nlls.py
, pro matlab nlls.m
. Animace se pauzuje po každém kroku, musíte klikat pro další iterace.
V řídícím skriptu se přijme krok gradientní metody pouze pokud se zlepší účelová funkce a na základě toho se upraví délka kroku, podobně jako u LM. Nenechte se tím zmást, tímto zaručujeme konvergenci metody. Více to diskutujeme v úkolech pro radost.