Алгоритм Левенберга — Марквардта

Алгоритм Левенберга–Марквардта (англ. Levenberg–Marquardt algorithm, LMA або просто LM), також відомий як метод сгасних найменших квадратів (англ. damped least-squares, DLS) використовується у математиці та обчислювальній техніці для розв'язування нелінійних задач найменших квадратів. Такі задачі мінімізації особливо актуальні при підборі кривої методом найменших квадратів. LMA інтерполює між алгоритмом Гаусса–Ньютона (GNA) та методом градієнтного спуску. LMA є більш надійним, ніж GNA, що означає, що в багатьох випадках він знаходить рішення, навіть якщо воно починається дуже далеко від кінцевого мінімуму. Для нормальної роботи функцій і розумних стартових параметрів LMA, як правило, повільніше, ніж GNA. LMA також можна розглядати як Гаусса–Ньютона, використовуючи підхід довіри до регіону.

Алгоритм був вперше опублікований у 1944 році Кеннетом Левенбергом, під час роботи у Франкфордському армійському арсеналі. У 1963 році його знову відкрили Дональд Марквардт, який працював статистиком у DuPont, і незалежно Жірард, Вінн і Моррісон.

LMA використовується в багатьох програмних додатках для розв'язання загальних задач допасовування кривої. Використовуючи алгоритм Гаусса–Ньютона, він часто сходиться швидше, ніж методи першого порядку. Однак, як і інші алгоритми ітераційної оптимізації, LMA знаходить лише локальний мінімум, який не обов'язково є глобальним мінімумом.

Проблема

Основне застосування алгоритму Левенберга–Марквардта полягає в задачі допасовування кривої методом найменших квадратів: для заданого набору $m$ емпіричних пар $\left(x_{i},y_{i}\right)$ незалежних і залежних змінних, знайти параметри ${\boldsymbol {\beta }}$ модельної кривої $f\left(x,{\boldsymbol {\beta }}\right)$ так, щоб суму квадратів відхилень $S\left({\boldsymbol {\beta }}\right)$ було зведено до мінімуму:

${\hat {\boldsymbol {\beta }}}\in \operatorname {argmin} \limits _{\boldsymbol {\beta }}S\left({\boldsymbol {\beta }}\right)\equiv \operatorname {argmin} \limits _{\boldsymbol {\beta }}\sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)\right]^{2},$ набір вважається непорожнім.

Рішення

Як і інші алгоритми чисельної мінімізації, алгоритм Левенберга–Марквардта є ітераційною процедурою. Щоб почати мінімізацію, користувач повинен надати початкове припущення для вектора параметрів ${\boldsymbol {\beta }}$. У випадках лише з одним мінімумом, ненавчені стандартні припущення, як ${\boldsymbol {\beta }}^{\text{T}}={\begin{pmatrix}1,\ 1,\ \dots ,\ 1\end{pmatrix}}$ буде працювати нормально; у випадках з кількома мінімумами алгоритм сходить до глобального мінімуму лише в тому випадку, якщо початкове припущення вже дещо близько до остаточного рішення.

На кожному кроці ітерації вектор параметрів ${\boldsymbol {\beta }}$ замінюється на нову оцінку ${\boldsymbol {\beta }}+{\boldsymbol {\delta }}$. Щоб визначити ${\boldsymbol {\delta }}$, функція $f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)$ апроксимується його лінеаризацією:

$f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx f\left(x_{i},{\boldsymbol {\beta }}\right)+\mathbf {J} _{i}{\boldsymbol {\delta }},$

де

$\mathbf {J} _{i}={\frac {\partial f\left(x_{i},{\boldsymbol {\beta }}\right)}{\partial {\boldsymbol {\beta }}}}$

є градієнтом (в даному випадку вектором рядка) $f$ з відношенням до ${\boldsymbol {\beta }}$.

Сума $S\left({\boldsymbol {\beta }}\right)$ квадратних відхилень має свій мінімум при нульовому градієнті по відношенню до ${\boldsymbol {\beta }}$. Наведене вище наближення першого порядку $f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)$ дає

$S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)\approx \sum _{i=1}^{m}\left[y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right)-\mathbf {J} _{i}{\boldsymbol {\delta }}\right]^{2},$

або у векторному позначенні,

${\begin{aligned}S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)&\approx \left\|\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right\|^{2}\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)-\mathbf {J} {\boldsymbol {\delta }}\right]\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]-\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}-\left(\mathbf {J} {\boldsymbol {\delta }}\right)^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]+{\boldsymbol {\delta }}^{\mathrm {T} }\mathbf {J} ^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}\\&=\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]-2\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}+{\boldsymbol {\delta }}^{\mathrm {T} }\mathbf {J} ^{\mathrm {T} }\mathbf {J} {\boldsymbol {\delta }}.\end{aligned}}$

Візьмемо похідну від $S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)$ з відношенням до ${\boldsymbol {\delta }}$ і встановлення результату на нуль, що дає

$\left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right){\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right],$

де $\mathbf {J}$ є матриця Якобі, у якій $i$-й ряд дорівнює $\mathbf {J} _{i}$, і де $\mathbf {f} \left({\boldsymbol {\beta }}\right)$ і $\mathbf {y}$ є векторами з $i$-ї компоненти $f\left(x_{i},{\boldsymbol {\beta }}\right)$ і $y_{i}$ відповідно. Наведений вище вираз отримано для ${\boldsymbol {\beta }}$ підпадає під метод Гаусса–Ньютона. Матриця Якобі, як визначено вище, є (загалом) не квадратною матрицею, а прямокутною матрицею розміру $m\times n$, де $n$ — кількість параметрів (розмір вектора ${\boldsymbol {\beta }}$). Матричне множення $\left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)$ дає необхідну $n\times n$ квадратну матрицю і добуток матриці-вектору з правого боку дають вектор розміру $n$. В результаті виходить набір $n$ лінійних рівнянь, які можна розв'язувати для ${\boldsymbol {\delta }}$.А

Внесок Левенберга полягає в тому, щоб замінити це рівняння на «згасну версію»:

$\left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \mathbf {I} \right){\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right],$

де $\mathbf {I}$ є одиничною матрицею, що дає приріст ${\boldsymbol {\delta }}$ до розрахункового вектора параметрів ${\boldsymbol {\beta }}$.

Коефіцієнт (невід'ємного) демпфування $\lambda$ коригується на кожній ітерації. Якщо зменшення $S$ є швидким, можна використовувати менше значення, наближаючи алгоритм до алгоритму Гаусса–Ньютона, тоді як, якщо ітерація дає недостатнє зменшення залишку, $\lambda$ можна збільшити, наблизивши на крок до напрямку градієнтного спуску. Зверніть увагу, що градієнт $S$ з відношенням до ${\boldsymbol {\beta }}$ дорівнює $-2\left(\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]\right)^{\mathrm {T} }$. Тому для великих значень $\lambda$, крок буде зроблено приблизно в напрямку, протилежному градієнту. Якщо будь-яка довжина обчисленого кроку ${\boldsymbol {\delta }}$ або зменшення суми квадратів з останнього вектора параметрів ${\boldsymbol {\beta }}+{\boldsymbol {\delta }}$ падають нижче попередньо визначених меж, ітерація зупиняється та останній вектор параметра ${\boldsymbol {\beta }}$ вважається рішенням.

Коли коефіцієнт демпфування $\lambda$ є великим відносно $\|\mathbf {J} ^{\mathrm {T} }\mathbf {J} \|$, інвертувати $\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \mathbf {I}$ не потрібно, оскільки оновлення добре апроксимується невеликим кроком градієнта $\lambda ^{-1}\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]$.

Щоб зробити масштаб рішення інваріантним, алгоритм Марквардта розв'язав модифіковану задачу з кожним компонентом градієнта, масштабованим відповідно до кривизни. Це забезпечує більший рух уздовж напрямків, де градієнт менший, що дозволяє уникнути повільної збіжності в напрямку малого градієнта. Флетчер у своїй статті 1971 року Модифікована підпрограма Марквардта для нелінійних найменших квадратів спростив форму, замінивши ідентичну матрицю $\mathbf {I}$ з діагональною матрицею, що складається з діагональних елементів $\mathbf {J} ^{\text{T}}\mathbf {J}$ :

$\left[\mathbf {J} ^{\mathrm {T} }\mathbf {J} +\lambda \operatorname {diag} \left(\mathbf {J} ^{\mathrm {T} }\mathbf {J} \right)\right]{\boldsymbol {\delta }}=\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right].$

Подібний коефіцієнт демпфування з'являється в регуляризації Тихонова, яка використовується для розв'язування лінійних неправильних задач, а також у регресії хребта, методиці оцінки в статистиці.

Вибір параметра демпфування

Для кращого вибору параметра демпфування були висунуті різні більш-менш евристичні аргументи $\lambda$. Існують теоретичні аргументи, які показують, чому деякі з цих варіантів гарантують локальну збіжність алгоритму; однак ці варіанти можуть призвести до того, що глобальна збіжність алгоритму страждає від небажаних властивостей найкрутішого спуску, зокрема, дуже повільної збіжності, близької до оптимальної.

Абсолютні значення будь-якого вибору залежать від того, наскільки добре масштабована початкова проблема. Марквардт рекомендував починати зі значення $\lambda _{0}$ і фактор $\nu >1$. Спочатку налаштування $\lambda =\lambda _{0}$ і обчислення залишкової суми квадратів $S\left({\boldsymbol {\beta }}\right)$ через один крок від початкової точки з коефіцієнтом демпфування $\lambda =\lambda _{0}$ а по-друге з $\lambda _{0}/\nu$. Якщо обидва ці показники гірші за початкову точку, то згасання збільшується шляхом послідовного множення на $\nu$ поки не буде знайдено кращу точку з новим коефіцієнтом демпфування $\lambda _{0}\nu ^{k}$ для деяких $k$.

Якщо використати коефіцієнт демпфування $\lambda /\nu$ це призводить до зменшення квадрата залишку, то це приймається за нове значення $\lambda$ (і нове оптимальне розташування приймається як те, що отримано з цим коефіцієнтом демпфування) і процес продовжується; якщо використовувати $\lambda /\nu$ це призведе до гіршого залишку, але використання $\lambda$ призведе до кращого залишку, тоді $\lambda$ залишається без змін, а новий оптимум приймається за значення, отримане с $\lambda$ як демпфуючий фактор.

Ефективна стратегія для контролю параметра демпфування, що називається відкладеним задоволенням, полягає у збільшенні параметра на невелику величину для кожного кроку на підйом і зменшення на велику величину для кожного кроку вниз. Ідея цієї стратегії полягає в тому, щоб уникнути занадто швидкого спуску на початку оптимізації, отже, обмежуючи кроки, доступні в майбутніх ітераціях, і, отже, сповільнюючи зближення. Збільшення в 2 рази і зменшення в 3 рази виявилося ефективним у більшості випадків, тоді як для великих проблем більш екстремальні значення можуть працювати краще, зі збільшенням у 1,5 раза та зменшенням у 5 раз.

Геодезичне прискорення

При інтерпретації кроку Левенберга–Марквардта як швидкості ${\boldsymbol {v}}_{k}$ вздовж геодезичної траєкторії в просторі параметрів можна покращити метод, додавши член другого порядку, що враховує прискорення ${\boldsymbol {a}}_{k}$ вздовж геодезичної

${\boldsymbol {v}}_{k}+{\frac {1}{2}}{\boldsymbol {a}}_{k}$

де ${\boldsymbol {a}}_{k}$ є рішенням

${\boldsymbol {J}}_{k}{\boldsymbol {a}}_{k}=-f_{vv}.$

Оскільки цей член геодезичного прискорення залежить лише від похідної за напрямком $f_{vv}=\sum _{\mu \nu }v_{\mu }v_{\nu }\partial _{\mu }\partial _{\nu }f({\boldsymbol {x}})$ вздовж напрямку швидкості ${\boldsymbol {v}}$, він не вимагає обчислення повної похідної матриці другого порядку, вимагаючи лише невеликих накладних обчислювальних витрат. Оскільки похідна другого порядку може бути досить складним виразом, може бути зручно замінити її наближенням скінченної різниці

${\begin{aligned}f_{vv}^{i}&\approx {\frac {f_{i}({\boldsymbol {x}}+h{\boldsymbol {\delta }})-2f_{i}({\boldsymbol {x}})+f_{i}({\boldsymbol {x}}-h{\boldsymbol {\delta }})}{h^{2}}}\\&={\frac {2}{h}}\left({\frac {f_{i}({\boldsymbol {x}}+h{\boldsymbol {\delta }})-f_{i}({\boldsymbol {x}})}{h}}-{\boldsymbol {J}}_{i}{\boldsymbol {\delta }}\right)\end{aligned}}$

де $f({\boldsymbol {x}})$ і ${\boldsymbol {J}}$ вже були обчислені за допомогою алгоритму, тому для обчислення потрібна лише одна додаткова оцінка функції $f({\boldsymbol {x}}+h{\boldsymbol {\delta }})$. Вибір скінченного різницевого кроку $h$ може вплинути на стабільність алгоритму, вибір значення близько 0,1 зазвичай є розумним.

Оскільки прискорення може вказувати в напрямку, протилежному швидкості, щоб запобігти зупинці методу, якщо демпфування занадто мале, додається додатковий критерій прискорення, щоб прийняти крок, який вимагає, що

${\frac {2\left\|{\boldsymbol {a}}_{k}\right\|}{\left\|{\boldsymbol {v}}_{k}\right\|}}\leq \alpha$

де $\alpha$ зазвичай фіксується до значення менше ніж 1, з меншими значеннями для складніших задач.

Додавання геодезичного члена прискорення може дозволити суттєво збільшити швидкість збіжності, і це особливо корисно, коли алгоритм рухається через вузькі каньйони в ландшафті цільової функції, де дозволені кроки менші, а точність вища завдяки другому порядку термін дає значні покращення.

Приклад

Погано підходить

Краще підходить

Найкраще підходить

У цьому прикладі ми намагаємося підібрати функції $y=a\cos \left(bX\right)+b\sin \left(aX\right)$ використовуючи алгоритм Левенберга–Марквардта, реалізований в GNU Octave, як функцію leasqr. Графіки показують, що вони все краще відповідають параметрам $a=100$, $b=102$ що використовується на початковій кривій. Лише тоді, коли параметри останнього графіка вибрано найближчими до оригіналу, криві точно підходять. Це рівняння є прикладом дуже чутливих початкових умов для алгоритму Левенберга–Марквардта. Однією з причин такої чутливості є існування кількох мінімумів — функції $\cos \left(\beta x\right)$ що має мінімуми при значенні параметра ${\hat {\beta }}$ і ${\hat {\beta }}+2n\pi$.

Див. також

• Довірча область
• Метод Нелдера – Міда
• Варіанти алгоритму Левенберга–Марквардта також використовувалися для розв'язування нелінійних систем рівнянь.

Література

• Moré, Jorge J.; Sorensen, Daniel C. (1983). Computing a Trust-Region Step (PDF). SIAM J. Sci. Stat. Comput. 4 (3): 553—572. doi:10.1137/0904038.
• Gill, Philip E.; Murray, Walter (1978). Algorithms for the solution of the nonlinear least-squares problem. SIAM Journal on Numerical Analysis. 15 (5): 977—992. Bibcode:1978SJNA...15..977G. doi:10.1137/0715063.
• Pujol, Jose (2007). The solution of nonlinear inverse problems and the Levenberg-Marquardt method. Geophysics. SEG. 72 (4): W1—W16. Bibcode:2007Geop...72W...1P. doi:10.1190/1.2732552.^{[недоступне посилання з 01.02.2020]}
• Nocedal, Jorge; Wright, Stephen J. (2006). Numerical Optimization (вид. 2nd). Springer. ISBN 978-0-387-30303-1.

Зовнішні посилання

• Numerical Recipes in C, Chapter 15.5: Nonlinear models
• C. T. Kelley, Iterative Methods for Optimization, SIAM Frontiers in Applied Mathematics, no 18, 1999, ISBN 0-89871-433-8. Online copy
• History of the algorithm in SIAM news
• A tutorial by Ananth Ranganathan
• K. Madsen, H. B. Nielsen, O. Tingleff, Methods for Non-Linear Least Squares Problems (nonlinear least-squares tutorial; L-M code: analytic Jacobian secant)
• T. Strutz: Data Fitting and Uncertainty (A practical introduction to weighted least squares and beyond). 2nd edition, Springer Vieweg, 2016, ISBN 978-3-658-11455-8.
• H. P. Gavin, The Levenberg-Marquardt method for nonlinear least-squares curve-fitting problems (MATLAB implementation included)