Метод Левенберга-Марквардта

Метод Левенберга-Марквардта – хороший выбор, если вам требуется минимизировать функцию вида F=f₁²(x₁,...,x_n)+ ... +f_m²(x₁,...,x_n). Алгоритм удачно сочетает в себе метод наискорейшего спуска (т.е. минимизации вдоль градиента) и метод Ньютона (т.е. использование квадратичной модели для ускорения поиска минимума функции). От метода наискорейшего спуска алгоритм позаимствовал стабильность работы, от метода Ньютона – ускоренную сходимость в окрестностях минимума. Ниже приведено обсуждение стандартной реализации алгоритма Левенберга-Марквардта, её недостатков, и улучшенной версии алгоритма, входящей в пакет ALGLIB. Перед чтением этого раздела рекомендуется ознакомиться с описанием алгоритма в Википедии или в Numerical Recipes. Далее мы предполагаем, что читающий понимает общие принципы работы алгоритма Левенберга-Марквардта.

Реализация метода Левенберга-Марквардта в ALGLIB

Традиционный вариант алгоритма Левенберга-Марквардта предполагает использование матрицы Якоби для построения квадратичной модели, осуществление одного шага в соответствии с моделью, затем – построение новой модели. У такого подхода есть ряд недостатков – общая неэффективность алгоритма, высокие требования к машинным ресурсам, ухудшение сходимости к минимуму на определенном классе задач. Рассмотрим их подробнее:

Неэффективность алгоритма объясняется влиянием нескольких факторов. Во-первых, в оригинальном алгоритме квадратичная модель функции обновляется после каждого нового шага, однако обновление модели – очень трудоемкая процедура. Можно использовать модель повторно и совершать несколько последовательных шагов. Эффективность таких шагов будет ниже, чем если бы мы обновили модель – но это компенсируется снижением трудоемкости за счет отказа от слишком частого обновления модели. Во-вторых, на ранних этапах оптимизации построение квадратичной модели просто нецелесообразно. Даже если её построить, алгоритм Левенберга-Марквардта всё равно будет совершать шаги, соответствующие методу наискорейшего спуска.
Высокие требования к машинным ресурсам. Очень часто число функций M намного больше размерности задачи N. Хранение матрицы Якоби требует O(M·N) машинной памяти, а получаемая в итоге квадратичная модель обычно имеет размер всего лишь O(N·N). На практике это различие может составлять от нескольких раз до нескольких десятков раз. Если известен гессиан функции F, то построение квадратичной модели на его основе позволяет значительно сократить требования к машинной памяти без существенного изменения быстродействия. А для некоторых задач (например, при построении нейросетевых классификаторов) гессиан может быть вычислен на порядок быстрее якобиана.
Сходимость к минимуму. В Numerical Recipes показано, что простота алгоритма Левенберга-Марквардта объясняется тем, что был отброшен ряд слагаемых второго порядка. В случае, если минимум достигается при малом значении F (good fit), отброшенные слагаемые компенсируют друг друга и сходимость алгоритма близка к квадратичной. Однако если минимум достигается в области больших значений F (bad fit), сходимость становится значительно хуже. Использование точного гессиана (если он доступен) позволяет решить эту проблему.

С учетом вышеизложенных замечаний в пакете ALGLIB реализован улучшенный вариант алгоритма.Новый алгоритм может использовать более полную информацию о задаче: значения функций f_i, матрицу Якоби, градиент функции F, гессиан функции F. В зависимости от того, какая информация доступна, алгоритм может использоваться в следующих вариантах:

FJ (function+Jacobian) – оригинальный алгоритм Левенберга-Марквардта без улучшений.
FGJ (function+gradient+Jacobian) – гибридный алгоритм с повторным использованием квадратичной модели. На первой стадии используется L-BFGS алгоритм для быстрой сходимости к области со сравнительно хорошими значениями функции, после этого используется метод Левенберга-Марквардта. После каждого шага алгоритма построенная квадратичная модель используется в качестве предобуславливателя для нескольких шагов L-BFGS-алгоритма.
FGH (function+gradient+Hessian) – гибридный алгоритм с использованием точного гессиана и с повторным использованием квадратичной модели. В большинстве случаев это – оптимальный выбор (если у вас есть возможность выбирать между использованием Якобиана или Гессиана).

Использование алгоритма и обратная коммуникация

Алгоритм оптимизации в ходе своей работы должен получать значения функции/градиента/... в выбранных им точках. В большинстве программных пакетов эта проблема решается путем передачи указателя на функцию (C++, Delphi) или делегата (C#), который осуществляет эту операцию.

Пакет ALGLIB, в отличие от других библиотек, использует для решения этой задачи обратную коммуникацию. Когда требуется вычислить значение функции (или её производных), состояние алгоритма сохраняется в специальной структуре, после чего управление возвращается в вызвавшую программу, которая осуществляет все вычисления и снова вызывает вычислительную подпрограмму.

Таким образом, работа с алгоритмом оптимизации осуществляется в следующей последовательности:

Подготовка структуры данных LMState при помощи одной из подпрограмм инициализации алгоритма (MinLMFJ, MinLMFGJ, MinLMFGH).
Установка дополнительных параметров алгоритма оптимизации при помощи подпрограмм MinLMSetXXXXX.
Вызов подпрограммы MinLMIteration.
Если подпрограмма вернула False, работа алгоритма завершена и минимум найден (сам минимум может быть получен при помощи подпрограммы MinLMResults).
Если подпрограмма вернула True, подпрограмма требует информацию о функции. В зависимости от того, какие поля структуры LMState установлены в True (ниже этот вопрос рассмотрен более подробно), вычислите функцию/градиент/якобиан/гессиан.
После того, как вся требуемая информация загружена в структуру LMState, требуется повторно вызвать подпрограмму MinLMIteration.

Для обмена информацией с пользователем используются следующие поля структуры LMState:

LMState.X[0..N-1] – массив, хранящий координаты точки, информация о которой запрашивается алгоритмом
LMState.F – в это поле следует поместить значение функции F (если оно было запрошено)
LMState.FI[0..M-1] – в это поле следует поместить вектор функций f_i (если он был запрошен)
LMState.G[0..N-1] – в это поле следует поместить градиент (если он был запрошен)
LMState.J[0..M-1,0..N-1] – в это поле следует поместить матрицу Якоби (если он был запрошен)
LMState.H[0..N-1,0..N-1] – в это поле следует поместить гессиан (если он был запрошен)

В зависимости от того, что именно требуется вычислить, подпрограмма MinLMIteration может устанавливать в True одно и только одно из следующих полей:

NeedF – сигнализирует о том, что требуется вычислить значение функции F (не путайте с вектором f_i)
NeedFG – сигнализирует о том, что требуется вычислить значения функции F и градиента F
NeedFiJ – сигнализирует о том, что требуется вычислить вектор функций f_i и матрицу Якоби
NeedFGH – сигнализирует о том, что требуется вычислить значение функции F, градиент F и гессиан F
XUpdated - установкой этого поля алгоритм сообщает о новой итерации. Новая точка содержится в State.X, а State.F содержит новое значение функции.

Ниже приведена таблица, показывающая, какие из полей могут быть установлены различными вариантами алгоритма Левенберга-Марквардта.

         NeedF   NeedFG  NeedFiJ NeedFGH XUpdated
MinLMFJ     X               X               X
MinLMFGJ    X       X       X               X
MinLMFGH    X       X               X       X