Алгоритм Левенберга – Марквардта - Levenberg–Marquardt algorithm

В математика и вычисления, Алгоритм Левенберга – Марквардта (LMA или просто LM), также известный как метод наименьших квадратов с затуханием (DLS), используется для решения нелинейный метод наименьших квадратов проблемы. Эти проблемы минимизации возникают особенно в наименьших квадратов подгонка кривой.

LMA используется во многих программных приложениях для решения общих задач аппроксимации кривой. Однако, как и во многих алгоритмах подгонки, LMA находит только местный минимум, что не обязательно глобальный минимум. LMA интерполирует между Алгоритм Гаусса – Ньютона (GNA) и метод градиентный спуск. LMA больше крепкий чем GNA, что означает, что во многих случаях он находит решение, даже если оно начинается очень далеко от конечного минимума. Для корректных функций и разумных начальных параметров LMA, как правило, немного медленнее, чем GNA. LMA также можно рассматривать как Гаусс – Ньютон используя регион доверия подход.

Алгоритм был впервые опубликован в 1944 г. Кеннет Левенберг,^[1] во время работы в Франкфордский армейский арсенал. Он был открыт заново в 1963 году. Дональд Марквардт,^[2] кто работал статистик в DuPont, и независимо от Жирара,^[3] Wynne^[4] и Моррисон.^[5]

Проблема

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

${\displaystyle {\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},}$ который предполагается непустым.

Решение

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

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

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

куда

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

это градиент (вектор-строка в данном случае) ${\displaystyle f}$ относительно ${\displaystyle {\boldsymbol {\beta }}}$.

Сумма ${\displaystyle S\left({\boldsymbol {\beta }}\right)}$ квадратичных отклонений имеет минимум в нуле градиент относительно ${\displaystyle {\boldsymbol {\beta }}}$. Приведенная выше аппроксимация первого порядка ${\displaystyle f\left(x_{i},{\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)}$ дает

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

или в векторной записи,

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

Взяв производную от ${\displaystyle S\left({\boldsymbol {\beta }}+{\boldsymbol {\delta }}\right)}$ относительно ${\displaystyle {\boldsymbol {\delta }}}$ и установка результата на ноль дает

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

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

Вклад Левенберга состоит в замене этого уравнения «версией с затуханием»:

${\displaystyle \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],}$

куда ${\displaystyle \mathbf {I} }$ - единичная матрица, дающая в качестве приращения ${\displaystyle {\boldsymbol {\delta }}}$ к оцениваемому вектору параметров ${\displaystyle {\boldsymbol {\beta }}}$.

(Неотрицательный) коэффициент демпфирования ${\displaystyle \lambda }$ настраивается на каждой итерации. Если сокращение ${\displaystyle S}$ быстро, можно использовать меньшее значение, приближая алгоритм к Алгоритм Гаусса – Ньютона, тогда как если итерация дает недостаточное уменьшение остатка, ${\displaystyle \lambda }$ можно увеличить, сделав шаг ближе к направлению градиентного спуска. Обратите внимание, что градиент из ${\displaystyle S}$ относительно ${\displaystyle {\boldsymbol {\beta }}}$ равно ${\displaystyle -2\left(\mathbf {J} ^{\mathrm {T} }\left[\mathbf {y} -\mathbf {f} \left({\boldsymbol {\beta }}\right)\right]\right)^{\mathrm {T} }}$. Следовательно, при больших значениях ${\displaystyle \lambda }$, шаг будет сделан примерно в направлении, противоположном градиенту. Если либо длина расчетного шага ${\displaystyle {\boldsymbol {\delta }}}$ или уменьшение суммы квадратов от последнего вектора параметров ${\displaystyle {\boldsymbol {\beta }}+{\boldsymbol {\delta }}}$ упадет ниже предопределенных пределов, итерация остановится, а вектор последнего параметра ${\displaystyle {\boldsymbol {\beta }}}$ считается решением.

Недостатком алгоритма Левенберга является то, что если значение коэффициента демпфирования ${\displaystyle \lambda }$ большой, инвертирующий ${\displaystyle \mathbf {J} ^{\text{T}}\mathbf {J} +\lambda \mathbf {I} }$ вообще не используется. Флетчер представил, что мы можем масштабировать каждый компонент градиента в соответствии с кривизной, чтобы было большее движение вдоль направлений, где градиент меньше. Это позволяет избежать медленного схождения в направлении небольшого градиента. Поэтому Флетчер в своей статье 1971 г. Модифицированная подпрограмма Марквардта для нелинейных наименьших квадратов заменил единичную матрицу ${\displaystyle \mathbf {I} }$ с диагональной матрицей, состоящей из диагональных элементов ${\displaystyle \mathbf {J} ^{\text{T}}\mathbf {J} }$, что делает масштаб решения инвариантным:

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

Аналогичный коэффициент демпфирования появляется в Тихоновская регуляризация, который используется для решения линейных некорректно поставленные проблемы, а также в регресс гребня, оценка техника в статистика.

Выбор параметра демпфирования

Были выдвинуты различные более или менее эвристические аргументы в пользу лучшего выбора параметра демпфирования. ${\displaystyle \lambda }$. Существуют теоретические аргументы, показывающие, почему некоторые из этих вариантов гарантируют локальную сходимость алгоритма; тем не менее, этот выбор может ухудшить глобальную сходимость алгоритма из-за нежелательных свойств крутой спуск, в частности, очень медленная сходимость, близкая к оптимальной.

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

Если использовать коэффициент демпфирования ${\displaystyle \lambda /\nu }$ приводит к уменьшению квадрата остатка, тогда это принимается как новое значение ${\displaystyle \lambda }$ (и новое оптимальное местоположение принимается как полученное с этим коэффициентом демпфирования), и процесс продолжается; при использовании ${\displaystyle \lambda /\nu }$ привел к худшему остатку, но с использованием ${\displaystyle \lambda }$ привело к лучшему остатку, тогда ${\displaystyle \lambda }$ остается неизменным, а новый оптимум принимается как значение, полученное с ${\displaystyle \lambda }$ как коэффициент демпфирования.

Эффективная стратегия управления параметром демпфирования, называемая запоздалое признание, состоит из небольшого увеличения параметра для каждого шага вверх и большого уменьшения для каждого шага вниз. Идея, лежащая в основе этой стратегии, заключается в том, чтобы избежать слишком быстрого спуска в начале оптимизации, тем самым ограничивая шаги, доступные в будущих итерациях, и, следовательно, замедляя сходимость.^[6] Было показано, что увеличение в 2 раза и уменьшение в 3 раза является эффективным в большинстве случаев, в то время как для больших задач могут работать более экстремальные значения с увеличением в 1,5 раза и уменьшением в раз. из 5.^[7]

Геодезическое ускорение

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

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

куда ${\displaystyle {\boldsymbol {a}}_{k}}$ это решение

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

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

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

куда ${\displaystyle f({\boldsymbol {x}})}$ и ${\displaystyle {\boldsymbol {J}}}$ уже вычислены алгоритмом, поэтому требуется только одна дополнительная оценка функции для вычисления ${\displaystyle f({\boldsymbol {x}}+h{\boldsymbol {\delta }})}$. Выбор конечно-разностного шага ${\displaystyle h}$ может повлиять на стабильность алгоритма, и обычно значение около 0,1 является разумным.^[7]

Поскольку ускорение может указывать в направлении, противоположном скорости, чтобы предотвратить остановку метода в случае слишком малого демпфирования, добавляется дополнительный критерий ускорения, чтобы принять шаг, требующий, чтобы

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

куда ${\displaystyle \alpha }$ обычно устанавливается на значение меньше 1, с меньшими значениями для более сложных задач.^[7]

Добавление члена геодезического ускорения может обеспечить значительное увеличение скорости сходимости, и это особенно полезно, когда алгоритм движется через узкие каньоны в ландшафте целевой функции, где разрешенные шаги меньше и выше точность из-за второго порядка термин дает значительные улучшения.^[7]

Пример

Плохой подбор

Лучше подходит

Наиболее подходящий

В этом примере мы пытаемся подобрать функцию ${\displaystyle y=a\cos \left(bX\right)+b\sin \left(aX\right)}$ с использованием алгоритма Левенберга – Марквардта, реализованного в GNU Octave как leasqr функция. Графики показывают все более точное соответствие параметрам. ${\displaystyle a=100}$, ${\displaystyle b=102}$ в исходной кривой. Только тогда, когда параметры на последнем графике выбраны наиболее близко к исходному, кривые точно соответствуют. Это уравнение является примером очень чувствительных начальных условий для алгоритма Левенберга – Марквардта. Одной из причин такой чувствительности является наличие нескольких минимумов - функция ${\displaystyle \cos \left(\beta x\right)}$ имеет минимумы при значении параметра ${\displaystyle {\hat {\beta }}}$ и ${\displaystyle {\hat {\beta }}+2n\pi }$.

Смотрите также

• Регион доверия
• Метод Нелдера – Мида
• Варианты алгоритма Левенберга – Марквардта также использовались для решения нелинейных систем уравнений.^[9]

Рекомендации

1. Левенберг, Кеннет (1944). «Метод решения некоторых нелинейных задач наименьших квадратов». Квартал прикладной математики. 2 (2): 164–168. Дои:10.1090 / qam / 10666.
2. Марквардт, Дональд (1963). "Алгоритм оценки нелинейных параметров методом наименьших квадратов". Журнал SIAM по прикладной математике. 11 (2): 431–441. Дои:10.1137/0111030. HDL:10338.dmlcz / 104299.
3. Жирар, Андре (1958). "Отрывок из Revue d'optique théorique et instrumentale". Rev. Opt. 37: 225–241, 397–424.
4. Винн, К. Г. (1959). «Проектирование линз на электронно-цифровом компьютере: I». Proc. Phys. Soc. Лондон. 73 (5): 777–787. Bibcode:1959PPS .... 73..777Вт. Дои:10.1088/0370-1328/73/5/310.
5. Моррисон, Дэвид Д. (1960). «Методы нелинейных задач наименьших квадратов и доказательства сходимости». Материалы семинара Лаборатории реактивного движения по программам слежения и определению орбиты: 1–9.
6. Транструм, Марка К; Махта, Бенджамин Б; Сетна, Джеймс П. (2011). «Геометрия нелинейных наименьших квадратов с приложениями к неаккуратным моделям и оптимизации». Физический обзор E. APS. 83 (3): 036701.
7. Транструм, Марка К; Сетна, Джеймс П. (2012). «Улучшения алгоритма Левенберга-Марквардта для нелинейной минимизации методом наименьших квадратов». arXiv:1201.5885.
8. «Нелинейная аппроксимация методом наименьших квадратов». Научная библиотека GNU. Архивировано из оригинал на 2020-04-14.
9. Канцов, Кристиан; Ямасита, Нобуо; Фукусима, Масао (2004). «Методы Левенберга – Марквардта со свойствами сильной локальной сходимости для решения нелинейных уравнений с выпуклыми ограничениями». Журнал вычислительной и прикладной математики. 172 (2): 375–397. Дои:10.1016 / j.cam.2004.02.013.

дальнейшее чтение

• Море, Хорхе Дж .; Соренсен, Дэниел К. (1983). «Вычисление шага доверительной области» (PDF). SIAM J. Sci. Стат. Вычислить. 4 (3): 553–572. Дои:10.1137/0904038.
• Gill, Philip E .; Мюррей, Уолтер (1978). «Алгоритмы решения нелинейной задачи наименьших квадратов». Журнал SIAM по численному анализу. 15 (5): 977–992. Bibcode:1978SJNA ... 15..977G. Дои:10.1137/0715063.
• Пухоль, Хосе (2007). «Решение нелинейных обратных задач и метод Левенберга-Марквардта». Геофизика. SEG. 72 (4): W1 – W16. Bibcode:2007Geop ... 72 Вт ... 1P. Дои:10.1190/1.2732552.^{[постоянная мертвая ссылка ]}
• Нокедаль, Хорхе; Райт, Стивен Дж. (2006). Численная оптимизация (2-е изд.). Springer. ISBN  978-0-387-30303-1.

внешняя ссылка

• Подробное описание алгоритма можно найти в Числовые рецепты на языке C, Глава 15.5: Нелинейные модели
• К. Т. Келли, Итерационные методы оптимизации, SIAM Frontiers in Applied Mathematics, № 18, 1999 г., ISBN  0-89871-433-8. Интернет-копия
• История алгоритма в новостях SIAM
• Учебник Ананта Ранганатана
• К. Мадсен, Х. Б. Нильсен, О. Тинглефф, Методы решения нелинейных задач наименьших квадратов (Учебное пособие по нелинейному методу наименьших квадратов; код L-M: аналитический якобиан секущий )
• Т. Струц: Подгонка данных и неопределенность (практическое введение в взвешенный метод наименьших квадратов и другие аспекты). 2-е издание, Springer Vieweg, 2016 г., ISBN  978-3-658-11455-8.
• Х. П. Гэвин, Метод Левенберга-Марквардта для нелинейных задач аппроксимации кривой методом наименьших квадратов (MATLAB реализация включена)