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

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

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

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

Проблема

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

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

Рішення

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

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

$f(x_i,\mathit{\beta }+\mathit{\delta })\approx f(x_i,\mathit{\beta })+{𝐉}_i\mathit{\delta },$

де

${𝐉}_i=\frac{\partial f(x_i,\mathit{\beta })}{\partial \mathit{\beta }}$

є градієнтом (в даному випадку вектором рядка) Шаблон:Tmath з відношенням до Шаблон:Tmath.

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

$S(\mathit{\beta }+\mathit{\delta })\approx \sum _{i=1}^m{[y_i-f(x_i,\mathit{\beta })-{𝐉}_i\mathit{\delta }]}^2,$

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

$\begin{array}{cc}S(\mathit{\beta }+\mathit{\delta })& \approx {\Vert 𝐲-𝐟\left(\mathit{\beta }\right)-𝐉\mathit{\delta }\Vert }^2\\ & ={[𝐲-𝐟\left(\mathit{\beta }\right)-𝐉\mathit{\delta }]}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)-𝐉\mathit{\delta }]\\ & ={[𝐲-𝐟\left(\mathit{\beta }\right)]}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)]-{[𝐲-𝐟\left(\mathit{\beta }\right)]}^{\mathrm{T}}𝐉\mathit{\delta }-{\left(𝐉\mathit{\delta }\right)}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)]+{\mathit{\delta }}^{\mathrm{T}}{𝐉}^{\mathrm{T}}𝐉\mathit{\delta }\\ & ={[𝐲-𝐟\left(\mathit{\beta }\right)]}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)]-2{[𝐲-𝐟\left(\mathit{\beta }\right)]}^{\mathrm{T}}𝐉\mathit{\delta }+{\mathit{\delta }}^{\mathrm{T}}{𝐉}^{\mathrm{T}}𝐉\mathit{\delta }.\end{array}$

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

$\left({𝐉}^{\mathrm{T}}𝐉\right)\mathit{\delta }={𝐉}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)],$

де $𝐉$ є матриця Якобі, у якій Шаблон:Tmath-й ряд дорівнює ${𝐉}_i$, і де $𝐟\left(\mathit{\beta }\right)$ і $𝐲$ є векторами з Шаблон:Tmath-ї компоненти $f(x_i,\mathit{\beta })$ і $y_i$ відповідно. Наведений вище вираз отримано для Шаблон:Tmath підпадає під метод Гаусса–Ньютона. Матриця Якобі, як визначено вище, є (загалом) не квадратною матрицею, а прямокутною матрицею розміру $m\times n$, де $n$ — кількість параметрів (розмір вектора $\mathit{\beta }$). Матричне множення $\left({𝐉}^{\mathrm{T}}𝐉\right)$ дає необхідну $n\times n$ квадратну матрицю і добуток матриці-вектору з правого боку дають вектор розміру $n$. В результаті виходить набір $n$ лінійних рівнянь, які можна розв'язувати для Шаблон:Tmath.А

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

$({𝐉}^{\mathrm{T}}𝐉+\lambda 𝐈)\mathit{\delta }={𝐉}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)],$

де Шаблон:Tmath є одиничною матрицею, що дає приріст Шаблон:Tmath до розрахункового вектора параметрів Шаблон:Tmath.

Коефіцієнт (невід'ємного) демпфування Шаблон:Tmath коригується на кожній ітерації. Якщо зменшення Шаблон:Tmath є швидким, можна використовувати менше значення, наближаючи алгоритм до алгоритму Гаусса–Ньютона, тоді як, якщо ітерація дає недостатнє зменшення залишку, Шаблон:Tmath можна збільшити, наблизивши на крок до напрямку градієнтного спуску. Зверніть увагу, що градієнт Шаблон:Tmath з відношенням до Шаблон:Tmath дорівнює $-2{\left({𝐉}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)]\right)}^{\mathrm{T}}$. Тому для великих значень Шаблон:Tmath, крок буде зроблено приблизно в напрямку, протилежному градієнту. Якщо будь-яка довжина обчисленого кроку Шаблон:Tmath або зменшення суми квадратів з останнього вектора параметрів Шаблон:Tmath падають нижче попередньо визначених меж, ітерація зупиняється та останній вектор параметра Шаблон:Tmath вважається рішенням.

Коли коефіцієнт демпфування Шаблон:Tmath є великим відносно $\Vert {𝐉}^{\mathrm{T}}𝐉\Vert$, інвертувати ${𝐉}^{\mathrm{T}}𝐉+\lambda 𝐈$ не потрібно, оскільки оновлення добре апроксимується невеликим кроком градієнта ${\lambda }^{-1}{𝐉}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)]$.

Щоб зробити масштаб рішення інваріантним, алгоритм Марквардта розв'язав модифіковану задачу з кожним компонентом градієнта, масштабованим відповідно до кривизни. Це забезпечує більший рух уздовж напрямків, де градієнт менший, що дозволяє уникнути повільної збіжності в напрямку малого градієнта. Флетчер у своїй статті 1971 року Модифікована підпрограма Марквардта для нелінійних найменших квадратів спростив форму, замінивши ідентичну матрицю Шаблон:Tmath з діагональною матрицею, що складається з діагональних елементів Шаблон:Tmath :

$[{𝐉}^{\mathrm{T}}𝐉+\lambda \mathrm{diag}\left({𝐉}^{\mathrm{T}}𝐉\right)]\mathit{\delta }={𝐉}^{\mathrm{T}}[𝐲-𝐟\left(\mathit{\beta }\right)].$

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

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

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

Абсолютні значення будь-якого вибору залежать від того, наскільки добре масштабована початкова проблема. Марквардт рекомендував починати зі значення Шаблон:Tmath і фактор Шаблон:Tmath. Спочатку налаштування $\lambda ={\lambda }_0$ і обчислення залишкової суми квадратів $S\left(\mathit{\beta }\right)$ через один крок від початкової точки з коефіцієнтом демпфування $\lambda ={\lambda }_0$ а по-друге з Шаблон:Tmath. Якщо обидва ці показники гірші за початкову точку, то згасання збільшується шляхом послідовного множення на Шаблон:Tmath поки не буде знайдено кращу точку з новим коефіцієнтом демпфування Шаблон:Tmath для деяких Шаблон:Tmath.

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

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

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

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

${𝒗}_k+\frac{1}{2}{𝒂}_k$

де ${𝒂}_k$ є рішенням

${𝑱}_k{𝒂}_k=-f_{vv}.$

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

$\begin{array}{cc}{f}_{vv}^i& \approx \frac{f_i(𝒙+h\mathit{\delta })-2f_i(𝒙)+f_i(𝒙-h\mathit{\delta })}{h^2}\\ & =\frac{2}{h}(\frac{f_i(𝒙+h\mathit{\delta })-f_i(𝒙)}{h}-{𝑱}_i\mathit{\delta })\end{array}$

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

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

$\frac{2\Vert {𝒂}_k\Vert }{\Vert {𝒗}_k\Vert }\le \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$.

Див. також

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

Література

• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journalШаблон:Dead link
• Шаблон:Cite book

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

• 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. 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.
• H. P. Gavin, The Levenberg-Marquardt method for nonlinear least-squares curve-fitting problems (MATLAB implementation included)