Algorytm Levenberga-Marquardta

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacji, wyszukiwania

Algorytm Levenberga-Marquardta jest obecnie jednym z najczęściej wykorzystywanych algorytmów optymalizacji nieliniowej, w szczególności w nieliniowym zadaniu najmniejszych kwadratów. Jest to algorytm iteracyjny, łączący w sobie cechy metody największego spadku i metody Gaussa-Newtona.

Sformułowanie problemu[edytuj | edytuj kod]

Mając daną serię danych (t_i, y_i) \in \mathbf{R}^{2}, gdzie i = 1, 2, \ldots, 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^2_i(\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}), \ldots, 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[edytuj | edytuj kod]

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ż[edytuj | edytuj kod]

Linki zewnętrzne[edytuj | edytuj kod]