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

В математике и вычислительной технике алгоритм Левенберга-Марквардта ( LMA или просто LM ), также известный как метод демпфированных наименьших квадратов ( 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 раза и уменьшением в 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$.

См. также

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

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

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

Внешние ссылки

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