Tekst źródłowy strony Algorytm Levenberga-Marquardta

'''Algorytm Levenberga-Marquardta''' – [[algorytm]] [[optymalizacja (matematyka)|optymalizacji nieliniowej]]. Jest to algorytm iteracyjny, łączący w sobie cechy [[metoda najszybszego spadku|metody największego spadku]] i [[metoda Newtona (optymalizacja)|metody Gaussa-Newtona]].

== Sformułowanie problemu ==
Mając daną serię danych <math>(t_i, y_i) \in \mathbf{R}^2,</math> gdzie <math>i = 1, 2, \dots, N,</math> szukamy dopasowania <math>\bar{y} = f(t|\mathbf{p}),</math> gdzie <math>\mathbf{p} \in \mathbf{R}^n</math> – wektor parametrów. Zakładamy, że ''najlepszym dopasowaniem'' jest to minimalizujące [[funkcjonał]]:
: <math>\chi^2(f) = \chi^2(\mathbf{p}) = \sum_{i=1}^N [y_i - f(t_i|\mathbf{p})]^2.</math>

Algorytm Levenberga-Marquardta w ogólności znajduje rozwiązanie zadania optymalizacji nieliniowej funkcji dającej się zapisać w postaci:
: <math>\Phi(\mathbf{x}) = \frac{1}{2}\sum_{i=1}^N r^2_i(\mathbf{x}),</math>

gdzie <math>\mathbf{x} \in \mathbf{R}^n</math> i zakładamy, że <math>N \geqslant n.</math> Jak łatwo zauważyć, funkcjonał <math>\chi^2</math> daje się zapisać w taki sposób. Dla uproszczenia, przedstawmy funkcje <math>r_i</math> jako wektor <math>\mathbf{r}(\mathbf{x}) = (r_1(\mathbf{x}), \dots, r_N(\mathbf{x}))</math> (zwany ''wektorem rezydualnym''). Wtedy <math>\Phi(\mathbf{x}) = \|\mathbf{r}(\mathbf{x})\|^2.</math> Pochodne funkcji <math>\Phi</math> można zapisać przy użyciu [[Macierz Jacobiego|Macierzy Jacobiego]] funkcji <math>\mathbf{r},</math> zdefiniowanego jako <math>\{\mathsf{J}(\mathbf{x})\}_{ij} = \frac{\partial r_i}{\partial x_j}(\mathbf{x}).</math> W ogólnym przypadku [[Gradient (matematyka)|gradient]] funkcji <math>\Phi</math> można zapisać:
: <math>\nabla\Phi(\mathbf{x}) = \sum_{i=1}^N r_i(\mathbf{x})\nabla r_i(\mathbf{x}) = \mathsf{J}(\mathbf{x})^\mathsf{T} \mathbf{r}(\mathbf{x}),</math>

a jej ''[[Macierz Hessego]]'':
: <math>\nabla^2\Phi(\mathbf{x}) = \mathsf{J}(\mathbf{x})^\mathsf{T}\mathsf{J}(\mathbf{x})+\sum_{i=1}^N r_j(\mathbf{x})\nabla^2 r_j(\mathbf{x}).</math>

W przypadku, gdy funkcje <math>r_j</math> można aproksymować funkcjami liniowymi w otoczeniu interesującego nas punktu (wtedy <math>\nabla^2 r_j(\mathbf{x})</math> jest bliskie zeru), lub gdy <math>r_j(\mathbf{x})</math> jest małe, hesjan funkcji <math>\Phi</math> przyjmuje prostszą postać:
: <math>\nabla^2\Phi(\mathbf{x}) = \mathsf{J}(\mathbf{x})^\mathsf{T}\mathsf{J}(\mathbf{x}),</math>

a więc hesjan można otrzymać wprost mając dany jakobian wektora rezydualnego <math>\mathbf{r}(\mathbf{x}),</math> co jest charakterystyczne dla zadania najmniejszych kwadratów.

== Opis metody ==
Najprostszym podejściem do problemu minimalizacji funkcji <math>\Phi</math> jest [[Metoda najszybszego spadku|metoda największego spadku]], opisana schematem:
: <math>\mathbf{x}_{i+1} = \mathbf{x}_i - \lambda \nabla\Phi(\mathbf{x}_i),</math>

która jest, w ogólnym przypadku, wolno zbieżna. Aby poprawić jej zbieżność, można skorzystać z wiedzy o drugiej pochodnej minimalizowanej funkcji w badanym punkcie. Jednym z możliwych podejść jest rozwinięcie gradientu minimalizowanej funkcji w [[Wzór Taylora|szereg Taylora]]:
: <math>\nabla\Phi(\mathbf{x}) = \nabla\Phi(\mathbf{x}_0) + (\mathbf{x}-\mathbf{x}_0)^\mathsf{T}\nabla^2 \Phi(\mathbf{x}_0) + \ldots</math>

i przyjęcie przybliżenia kwadratowego funkcji <math>\Phi</math> w otoczeniu <math>\mathbf{x}_0</math> do rozwiązania równania <math>\nabla\Phi(\bar{\mathbf{x}}) = 0.</math> W ten sposób otrzymujemy [[metoda Newtona (optymalizacja)|metodę Gaussa-Newtona]] opisaną schematem:
: <math>\mathbf{x}_{i+1}=\mathbf{x}_i - (\nabla^2\Phi(\mathbf{x}_i))^{-1}\nabla\Phi(\mathbf{x}_i),</math>

gdzie hesjan funkcji <math>\Phi</math> nie musi być znany dokładnie i często wystarczy podane wcześniej przybliżenie. Niestety, szybkość zbieżności tej metody zależy od wyboru punktu początkowego, a konkretnie od liniowości minimalizowanej funkcji w otoczeniu punktu startowego. [[Kenneth Levenberg]] zauważył, że opisane metody (największego spadku i Gaussa-Newtona) nawzajem się uzupełniają i zaproponował następującą modyfikację kroku metody:
: <math>\mathbf{x}_{i+1} = \mathbf{x}_i - (\mathsf{H}(\mathbf{x}_i)+\lambda\mathsf{I})^{-1}\nabla\Phi(\mathbf{x}_i),</math>(*)

wraz z następującym algorytmem:
# oblicz wartość <math>\mathbf{x}_{i+1}</math> na podstawie <math>\mathbf{x}_i</math> i równania (*),
# oblicz wartość błędu w punkcie <math>\mathbf{x}_{i+1},</math>
# jeśli błąd wzrósł, wróć do wartości <math>\mathbf{x}_i,</math> zwiększ wartość <math>\lambda</math> <math>k</math>-krotnie i wróć do kroku 1 (przybliżenie liniowe minimalizowanej funkcji w otoczeniu <math>\mathbf{x}_i</math> okazało się nie dość ścisłe, więc zwiększamy „wpływ” metody największego spadku),
# jeśli błąd zmalał, zaakceptuj ten krok i zmniejsz wartość <math>\lambda</math> <math>k</math>-krotnie (założenie o liniowości minimalizowanej funkcji w otoczeniu <math>\mathbf{x}_i</math> okazało się wystarczająco ścisłe, więc zwiększamy „wpływ” metody Gaussa-Newtona).

W typowych zastosowaniach <math>k = 10.</math> W przypadku, gdy <math>\lambda</math> jest duże, hesjan w zasadzie nie jest wykorzystywany. [[Donald Marquardt]] zauważył, że nawet w takiej sytuacji można wykorzystać informację zawartą w drugiej pochodnej minimalizowanej funkcji, poprzez skalowanie każdego komponentu wektora gradientu w zależności od krzywizny w danym kierunku (co pomaga w źle uwarunkowanych zadaniach minimalizacji typu ''error valley''). Po uwzględnieniu poprawki Marquardta otrzymujemy następującą postać kroku metody:
: <math>\mathbf{x}_{i+1} = \mathbf{x}_i - (\mathsf{H}(\mathbf{x}_i)+\lambda\textrm{ diag}[\mathsf{H}])^{-1}\nabla\Phi(\mathbf{x}_i),</math>

gdzie:
: <math>\textrm{diag}[\mathsf{H}] = \left[\begin{array}{cccc}h_{11} & 0 & \ldots & 0 \\
0 & h_{22} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & h_{nn}\end{array}\right].</math>

Największą zaletą algorytmu Levenberga-Marquardta jest jego szybka zbieżność, w porównaniu z konkurencyjnymi metodami. Najkosztowniejszą operacją jest natomiast wyznaczenie macierzy odwrotnej, które w praktyce jest przeprowadzane w sposób przybliżony, na przykład przy użyciu metody [[Rozkład według wartości osobliwych|SVD]]. Tym niemniej, nawet w najoszczędniejszych przypadkach koszt jednego kroku rośnie niedopuszczalnie wraz ze wzrostem rozmiaru zadania powyżej tysiąca parametrów. Z drugiej dla zadań o umiarkowanej ilości parametrów (rzędu kilkuset), metoda Levenberga-Marquardta jest dużo szybsza od metody największego spadku.

== Zobacz też ==
* [[optymalizacja (matematyka)]]

== Linki zewnętrzne ==
* {{Cytuj stronę | url = http://www.fizyka.umk.pl/nrbook/c15-5.pdf | tytuł = Numerical Recipes in C, Chapter 15.5: Nonlinear models (en) | opublikowany = fizyka.umk.pl | archiwum = https://web.archive.org/web/20090419053311/http://www.fizyka.umk.pl/nrbook/c15-5.pdf | zarchiwizowano = 2009-04-19}}

{{SORTUJ:Levenberg-Marquardt}}
[[Kategoria:Algorytmy numeryczne]]
[[Kategoria:Metody optymalizacji]]