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 ${\displaystyle (t_{i},y_{i})\in \mathbf {R} ^{2},}$ gdzie ${\displaystyle i=1,2,\dots ,N,}$ szukamy dopasowania ${\displaystyle {\bar {y}}=f(t|\mathbf {p} ),}$ gdzie ${\displaystyle \mathbf {p} \in \mathbf {R} ^{n}}$ – wektor parametrów. Zakładamy, że najlepszym dopasowaniem jest to minimalizujące funkcjonał:

${\displaystyle \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:

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

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

${\displaystyle \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:

${\displaystyle \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 ${\displaystyle r_{j}}$ można aproksymować funkcjami liniowymi w otoczeniu interesującego nas punktu (wtedy ${\displaystyle \nabla ^{2}r_{j}(\mathbf {x} )}$ jest bliskie zeru), lub gdy ${\displaystyle r_{j}(\mathbf {x} )}$ jest małe, hesjan funkcji ${\displaystyle \Phi }$ przyjmuje prostszą postać:

${\displaystyle \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 ${\displaystyle \mathbf {r} (\mathbf {x} ),}$ co jest charakterystyczne dla zadania najmniejszych kwadratów.

Opis metody

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

${\displaystyle \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:

${\displaystyle \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 ${\displaystyle \Phi }$ w otoczeniu ${\displaystyle \mathbf {x} _{0}}$ do rozwiązania równania ${\displaystyle \nabla \Phi ({\bar {\mathbf {x} }})=0.}$ W ten sposób otrzymujemy metodę Gaussa-Newtona opisaną schematem:

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

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

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

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

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

gdzie:

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