Algorytm Levenberga-Marquardta

Algorytm Levenberga-Marquardta – algorytm optymalizacji nieliniowej. Jest to algorytm iteracyjny, łączący w sobie cechy metody największego spadku i metody Gaussa-Newtona.

Sformułowanie problemu

Mając daną serię danych $(t_{i},y_{i})\in \mathbf {R} ^{2},$ gdzie $i=1,2,\dots ,N,$ szukamy dopasowania ${\bar {y}}=f(t|\mathbf {p} ),$ gdzie $\mathbf {p} \in \mathbf {R} ^{n}$ – wektor parametrów. Zakładamy, że najlepszym dopasowaniem jest to minimalizujące funkcjonał:

$\chi ^{2}(f)=\chi ^{2}(\mathbf {p} )=\sum _{i=1}^{N}[y_{i}-f(t_{i}|\mathbf {p} )]^{2}.$

Algorytm Levenberga-Marquardta w ogólności znajduje rozwiązanie zadania optymalizacji nieliniowej funkcji dającej się zapisać w postaci:

$\Phi (\mathbf {x} )={\frac {1}{2}}\sum _{i=1}^{N}r_{i}^{2}(\mathbf {x} ),$

gdzie $\mathbf {x} \in \mathbf {R} ^{n}$ i zakładamy, że $N\geqslant n.$ Jak łatwo zauważyć, funkcjonał $\chi ^{2}$ daje się zapisać w taki sposób. Dla uproszczenia, przedstawmy funkcje $r_{i}$ jako wektor $\mathbf {r} (\mathbf {x} )=(r_{1}(\mathbf {x} ),\dots ,r_{N}(\mathbf {x} ))$ (zwany wektorem rezydualnym). Wtedy $\Phi (\mathbf {x} )=\|\mathbf {r} (\mathbf {x} )\|^{2}.$ Pochodne funkcji $\Phi$ można zapisać przy użyciu Macierzy Jacobiego funkcji $\mathbf {r} ,$ zdefiniowanego jako $\{{\mathsf {J}}(\mathbf {x} )\}_{ij}={\frac {\partial r_{i}}{\partial x_{j}}}(\mathbf {x} ).$ W ogólnym przypadku gradient funkcji $\Phi$ można zapisać:

$\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} ),$

a jej Macierz Hessego:

$\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} ).$

W przypadku, gdy funkcje $r_{j}$ można aproksymować funkcjami liniowymi w otoczeniu interesującego nas punktu (wtedy $\nabla ^{2}r_{j}(\mathbf {x} )$ jest bliskie zeru), lub gdy $r_{j}(\mathbf {x} )$ jest małe, hesjan funkcji $\Phi$ przyjmuje prostszą postać:

$\nabla ^{2}\Phi (\mathbf {x} )={\mathsf {J}}(\mathbf {x} )^{\mathsf {T}}{\mathsf {J}}(\mathbf {x} ),$

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

Opis metody

Najprostszym podejściem do problemu minimalizacji funkcji $\Phi$ jest metoda największego spadku, opisana schematem:

$\mathbf {x} _{i+1}=\mathbf {x} _{i}-\lambda \nabla \Phi (\mathbf {x} _{i}),$

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 szereg Taylora:

$\nabla \Phi (\mathbf {x} )=\nabla \Phi (\mathbf {x} _{0})+(\mathbf {x} -\mathbf {x} _{0})^{\mathsf {T}}\nabla ^{2}\Phi (\mathbf {x} _{0})+\ldots$

i przyjęcie przybliżenia kwadratowego funkcji $\Phi$ w otoczeniu $\mathbf {x} _{0}$ do rozwiązania równania $\nabla \Phi ({\bar {\mathbf {x} }})=0.$ W ten sposób otrzymujemy metodę Gaussa-Newtona opisaną schematem:

$\mathbf {x} _{i+1}=\mathbf {x} _{i}-(\nabla ^{2}\Phi (\mathbf {x} _{i}))^{-1}\nabla \Phi (\mathbf {x} _{i}),$

gdzie hesjan funkcji $\Phi$ 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:

$\mathbf {x} _{i+1}=\mathbf {x} _{i}-({\mathsf {H}}(\mathbf {x} _{i})+\lambda {\mathsf {I}})^{-1}\nabla \Phi (\mathbf {x} _{i}),$(*)

wraz z następującym algorytmem:

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

W typowych zastosowaniach $k=10.$ W przypadku, gdy $\lambda$ 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:

$\mathbf {x} _{i+1}=\mathbf {x} _{i}-({\mathsf {H}}(\mathbf {x} _{i})+\lambda {\textrm {diag}}[{\mathsf {H}}])^{-1}\nabla \Phi (\mathbf {x} _{i}),$

gdzie:

${\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].$

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 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

• Numerical Recipes in C, Chapter 15.5: Nonlinear models (en). fizyka.umk.pl. [zarchiwizowane z tego adresu (2009-04-19)].