Algorytm Levenberga-Marquardta

Algorytm Levenberga-Marquardta to technika optymalizacji nieliniowej, która działa w sposób iteracyjny, łącząc aspekty metody największego spadku oraz metody Gaussa-Newtona.

Sformułowanie problemu

Rozważmy zestaw danych opisanych jako:

(t_i, y_i) ∈ R², gdzie i = 1, 2, …, N. Szukamy dopasowania, które można zapisać jako:

ȳ = f(t | p),

gdzie p ∈ Rⁿ oznacza wektor parametrów. Zakładamy, że najlepszym dopasowaniem jest to, które minimalizuje funkcjonał:

χ²(f) = χ²(p) = ∑_i=1^N [y_i − f(t_i | p)]².

Algorytm Levenberga-Marquardta ogólnie znajduje rozwiązanie problemu optymalizacji dla funkcji, którą można zapisać w formie:

Φ(x) = 1/2 ∑_i=1^N r_i²(x),

gdzie x ∈ Rⁿ oraz zakładamy, że N ≥ n. Jak można zauważyć, funkcjonał χ² można zapisać w ten sposób. Dla uproszczenia, przedstawmy funkcje r_i jako wektor r(x) = (r₁(x), …, r_N(x)). Wtedy:

Φ(x) = ‖r(x)‖².

Pochodne funkcji Φ można wyrazić za pomocą macierzy Jacobian funkcji r, zdefiniowanej jako:

{J(x)}_ij = ∂r_i/∂x_j(x).

W ogólnym przypadku gradient funkcji Φ można zapisać:

∇Φ(x) = ∑_i=1^N r_i(x) ∇r_i(x) = J(x)^Tr(x),

a jej macierz Hessego:

∇2Φ(x) = J(x)^TJ(x) + ∑_i=1^N r_j(x) ∇2r_j(x).

W przypadku, gdy funkcje r_j można aproksymować funkcjami liniowymi w otoczeniu interesującego nas punktu, to hesjan funkcji Φ przyjmuje prostszą postać:

∇2Φ(x) = J(x)^TJ(x),

co jest charakterystyczne dla zadań najmniejszych kwadratów.

Opis metody

Najprostszym podejściem do minimalizacji funkcji Φ jest metoda największego spadku, opisana w następujący sposób:

x_i+1 = x_i − λ ∇Φ(x_i),

która w ogólnym przypadku jest wolno zbieżna. Aby poprawić zbieżność, można wykorzystać wiedzę o drugiej pochodnej minimalizowanej funkcji w badanym punkcie. Jednym z podejść jest rozwinięcie gradientu minimalizowanej funkcji w szereg Taylora:

∇Φ(x) = ∇Φ(x₀) + (x − x₀)^T∇2Φ(x₀) + …

I przyjęcie przybliżenia kwadratowego funkcji Φ w otoczeniu x₀ do rozwiązania równania:

∇Φ(x̄) = 0.

W ten sposób uzyskujemy metodę Gaussa-Newtona, opisaną jako:

x_i+1 = x_i − (∇2Φ(x_i))⁻¹∇Φ(x_i),

gdzie hesjan funkcji Φ nie musi być znany dokładnie, a często wystarczy podane wcześniej przybliżenie. Niestety, szybkość zbieżności tej metody zależy od wyboru punktu startowego oraz liniowości minimalizowanej funkcji w tym otoczeniu. Kenneth Levenberg zauważył, że metody największego spadku i Gaussa-Newtona się uzupełniają, i zaproponował modyfikację kroku metody:

x_i+1 = x_i − (H(x_i) + λI)⁻¹∇Φ(x_i),

gdzie λ to parametr regularyzacji, a I to macierz jednostkowa. W ten sposób uzyskujemy algorytm, w którym obliczamy wartość x_{i+1 na podstawie x_i i równania.}

Jeżeli błąd wzrasta, wracamy do wartości x_i, zwiększamy wartość λ i powtarzamy krok. W przeciwnym przypadku, akceptujemy krok i zmniejszamy λ. W typowych zastosowaniach przyjmuje się k = 10. Gdy λ jest duże, hesjan nie jest wykorzystywany, a Donald Marquardt zauważył, że w takiej sytuacji można skalować każdy komponent wektora gradientu w zależności od krzywizny w danym kierunku, co jest korzystne w zadaniach minimizacji.