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

Алгоритм Левенберга-Марквардта прост. Алгоритм Левенберга-Марквардта эффективен.

А еще о нем говорят, что он где-то посередине между градиентным спуском и методом Ньютона, что бы это не значило. Ну, с методом Ньютоном и его связью с градиентным спуском вроде как разобрались. Но что имеют ввиду когда произносят эту глубокомысленную фразу? Попробуем слегка подразобраться.

В своих статьях товарищ Левенберг [Levenberg, K. A Method for the Solution of Certain Problems in Last Squares. Quart. Appl. Math. 1944. Vol. 2. P. 164—168.], а после него гражданин Марквардт [Marquardt, Donald (1963). «An Algorithm for Least-Squares Estimation of Nonlinear Parameters». SIAM Journal on Applied Mathematics. 11 (2): 431–441.] рассмотрели задачу наименьших квадратов, которая выглядит так:

e478ffbd37ce67d12a8efa37326ce496.gif,

которую можно записать проще в векторной форме

43ae3d87f8c0662510f39e6a7e1453a0.gif.

А можно еще проще, полностью забив на наименьшие квадраты. Это никак не повлияет на повествование.

Итак, рассматривается задача

573d8c243893d6b70f894859eebc6e22.gif.

Такая задача возникает настолько часто, что важность нахождения эффективного метода ее решения сложно переоценить. Но мы начнем с другого. В предыдущей статье было показано, что широко известный метод градиентного спуска, и не только он, может быть получен из следующих соображений. Допустим, что мы пришли в некоторую точку 7790dd0efb4a03a4c876741804d9b559.gif, в которой минимизируемая функция имеет значение 90340615fd75f4a3550a82c374838b6b.gif. Определим в этой точке вспомогательную функцию 8bf1d54e1f36dd4c9dfd5720437af51c.gif, а также некоторую ее модель c5b1604002da7b2c951dd57929933d24.gif. Для данной модели мы ставим вспомогательную задачу

5ab24f2cdd939fb1186c5ebc91807432.gif

где 7e08381eec55a31db8263ce4d9b04120.gif — некоторое наперед заданное множество допустимых значений, выбираемое так, чтобы задача имела простое решение и при этом функция 462957dda265f4fb8be04327f1c12b0f.gif достаточно точно аппроксимировала da77c5b4891cf3d059f1b04a28b230ef.gif на 7e08381eec55a31db8263ce4d9b04120.gif. Такую схему называют методом доверительного региона, а множество 7e08381eec55a31db8263ce4d9b04120.gif, на котором минимизируется значение модельной функции — доверительным регионом этой функции. Для градиентного спуска мы брали c7e85c48eb16123c23a9e08714f50a0e.gif, для метода Ньютона 7cfb4b9203c4faccd18c1837b9c0e59f.gif, а в качестве модели для da77c5b4891cf3d059f1b04a28b230ef.gif выступала линейная часть разложения в ряд Тейлора a98a0849fa75ac2a6bbcbbc2bbda3054.gif.

Посмотрим, что будет, если усложнить модель, взяв

b2ad63f0162c3cad2512d187035cbf2d.gif.

Минимизируем эту модельную функцию на эллиптическом доверительном регионе e99c58ecc60ed62c044f7690e444d1d2.gif (множитель добавлен для удобства вычислений). Применив метод множителей Лагранжа, получим задачу

44c4c6050ba6f69f94a8c222d68c95f3.gif,

решение которой удовлетворяет равенству

4a3851c834bc50c9b670673f6665f792.gif

или

2c47f53d4e43da14ea21b9fe9d7bb533.gif

Здесь, в отличие от того, что мы видели раньше при использовании линейной модели, направление p оказывается зависимым не только от метрики dabea901c4b1a4079aa96d47bcee4e75.gif, но и от выбора размера доверительного региона a2b0686adbc103ad9f96be85cca5d418.gif, а значит методика линейного поиска неприменима (по крайней мере обоснованно). Также оказывается проблемным определить в явном виде величину 99d394e7d0b74248114405067e0ffd51.gif, соответствующую величине a2b0686adbc103ad9f96be85cca5d418.gif. Однако вполне очевидно, что при увеличении 99d394e7d0b74248114405067e0ffd51.gif длина 4b28c13d5f5d658adb7478fbc9efc923.gif будет уменьшаться. Если при этом мы еще наложим условие d6414f5fbe0850f1d6cd0710c69a89fe.gif, то длина шага будет не больше, чем та, которую дал бы метод Ньютона (всамделешный, без модификаций и условий).

Таким образом, мы можем вместо того, чтобы для заданного a2b0686adbc103ad9f96be85cca5d418.gif искать подходящее значение 99d394e7d0b74248114405067e0ffd51.gif, поступить с точностью до наоборот: найти такое d6414f5fbe0850f1d6cd0710c69a89fe.gif, при котором выполняется условие 55380bdc5a434366df6d181078d6a8b7.gif. Это своего рода замена почившему в данном случае линейному поиску. Марквардт предложил следующую простую процедуру:

  1. если для некотрого значения 99d394e7d0b74248114405067e0ffd51.gif условие eca77f2af789bf09851ed403e71813c5.gif выполнено, то повторять cc35f335e8bf25ccc19e392f3311591a.gif до тех пор, пока 3aadf929956784b2e081f8a634cff90f.gif
  2. если же 22b75c29456939ef738c81e53f52db6b.gif, то принять 434e0a742c4ad3ae7e4cea8c1289d363.gif и повторить.


Здесь 4f268e81e07eb6a87b5798903d82f2c0.gif и 0ead144f2f1591e5747e51404639631f.gif — константы, являющиеся параметрами метода. Умножение на 389a9983ea24ad0b3af0559c2aca381b.gif соответствует расширению доверительного региона, а умножение на 76d0eb69ba026a58bbe3edd275fee712.gif — его сужению.

Указанная методика может быть применена к любой целевой функции. Заметьте, здесь уже не требуется положительная определенность гессиана в отличие от рассмотренного ранее случая, когда метод Ньютона представлялся частным случаем метода последовательного спуска. Не требуется даже его невырожденность, что в ряде случаев очень важно. Однако в этом случае увеличивается цена поиска направления, поскольку каждое изменение 99d394e7d0b74248114405067e0ffd51.gif приводит к необходимости решать линейную систему для определения 4b28c13d5f5d658adb7478fbc9efc923.gif.

Посмотрим что будет, если применить данный подход к задаче о наименьших квадратах.

Градиент функции f58e0b7a5457a6f69dd3c16e2d5aadc5.gif, ее гессиан 9da93ef0ae8eda326c977d890c97beac.gif, где c7a31446512540e721403b61686e91ba.gif. Подставляем и получаем следующую систему, определяющую направление поиска

148b5606f695c8a52b629f69f8c167c5.gif.

Вполне приемлемо, но вычислять вторые производные вектор-функции может быть довольно накладно. Марквардт для обхода этой проблемы предложил использовать не саму функцию 188ee644e8202aad30eac11166858841.gif, а ее линейную аппросимацию 205c97755529c6a19a04851f7ec3d7f7.gif, при которой матрица 338ec0451e1b4b7e7decd0b4443a8828.gif обращается в ноль. Если теперь в качестве dabea901c4b1a4079aa96d47bcee4e75.gif взять единичную матрицу 80932643e100c59ad091cdf4d90a2bd5.gif, то получим стандартную форму метода Левенберга-Марквардта для решения задачи наименьших квадратов:

d3e6e448ae5465d70cded093a02fdb69.gif.

Для данного способа определения направления спуска Марквардтом же была доказана теорема о том, что при устремлении 99d394e7d0b74248114405067e0ffd51.gif к бесконечности направление 4b28c13d5f5d658adb7478fbc9efc923.gif стремится к антиградиенту. Строгое доказательство заинтересованный читатель сможет найти в базовой статье, но надеюсь, что само это утверждение стало достаточно очевидным из логики построения метода. Оно в определенной мере оправдывает вездесущую отсылку к тому, что при увеличении лямбды (которую по непонятной мне причине часто называют параметром регуляризации) мы получаем градиентный спуск. На самом деле ничего подобного — мы получили бы его только в пределе, в том самом, где длина шага стремится к нулю. Намного важнее то, что при достаточно большом значении лямбды направление, которое мы получаем, будет являться направлением спуска, а значит мы получаем глобальную сходимость метода. А вот вторая часть утверждения, что при устремлении лямбды к нулю мы получаем метод Ньютона, со всей очевидностью верна, но только если принять вместо 188ee644e8202aad30eac11166858841.gif ее линейную аппроксимацию 9dc1faa9cf4e8d569afaf161ac627568.gif.

Казалось бы, всё. Минимизируем норму вектор-функции в эллиптической метрике — используем Левенберга-Марквардта. Имеем дело с функцией общего вида и имеем возможность вычислить матрицу вторых производных — велкам использовать метод доверительного региона общего вида. Но есть же извращенцы…

Иногда методом Левенберга-Марквардта для минимизации функции 188ee644e8202aad30eac11166858841.gif называют выражение вот такого вида:

438a2f22ba381c85087a26d2125d1dac.gif.

Вроде все то же самое, но здесь bc819017bab0b9f9d995f262f3f76a42.gif — матрица вторых! производных функции 188ee644e8202aad30eac11166858841.gif. Формально это имеет право на существование, но является извращением. И вот почему. Тот же Марквардт в своей статье предложил метод решения системы уравнений 0a0ec780406efe57ca6444290ccfde09.gif путем минимизации функции 128aa60b3bc0593a797bd9ebd308b402.gif описанным методом. Если в качестве 59b4640f5bad6b14066d718cf44e9f9c.gif взять градиент целевой функции, то действительно получим приведенное выражение. А извращение это потому, что 

решается задача минимизации, порождаемая системой нелинейных уравнений, порождаемых задачей минимизации.

Даблстрайк. Такое выражение, как минимум, не лучше первого уравнения сферического доверительного региона, а вообще намного хуже как с точки зрения производительности (лишние операции по умножению, а в нормальных реализациях — факторизации), так и с точки зрения устойчивости метода (умножение матрицы на себя ухудшает ее обусловленность). Здесь иногда возражают, что e98c6bbd7a088cd46c2ed68aac4943ba.gif гарантированно положительно определена, но в данном случае это не имеет никакого значения. Давайте посмотрим на метод Левенберга-Марквардта с позиций метода последовательного спуска. В этом случае мы, получается, хотим в качестве метрики использовать матрицу 64bdbedcbd61e0741f92025fb7c4a1f8.gif, и чтобы она могла выступать в этом качестве, значение 99d394e7d0b74248114405067e0ffd51.gif должно обеспечивать ее положительную определенность. Учитывая, что dabea901c4b1a4079aa96d47bcee4e75.gif положительно определена, нужное значение 99d394e7d0b74248114405067e0ffd51.gif всегда может быть найдено —, а значит никакой необходимости требовать от bc819017bab0b9f9d995f262f3f76a42.gif положительной определенности не наблюдается.

В качестве матрицы dabea901c4b1a4079aa96d47bcee4e75.gif не обязательно брать единичную, но для квадратичной модели целевой функции указать адекватный доверительный регион уже не так просто, как для линейной модели. Если брать эллиптический регион, индуцированный гессианом, то метод вырождается в метод Ньютона (ну, почти)

f31f2c82f904fd2c95ab33daa112168b.gif

Если, конечно, матрица Гессе положительно определена. Если нет — то как и раньше можно в качестве метрики использовать исправленный гессиан, либо некоторую матрицу, к нему в каком-либо смысле близкую. Встречается также рекомендация использовать в качестве метрики матрицу e83ae41347e08ce41ff17ee556a7e06b.gif, которая по построению гарантированно является положительно определенной. К сожалению, мне не известно хоть сколь-нибудь строгого обоснования данного выбора, но в качестве эмпирической рекомендации он упоминается довольно часто.

В качестве иллюстрации давайте посмотрим, как ведет себя метод на все той же функции Розенброка, причем мы будем рассматривать ее в двух ипостасях — как простую функцию, записанную в форме

b78272e1a4bce3b65e55a7737fd2a1ee.gif,

и как задачу наименьших квадратов

05cc5a49e94d86e2929c86264cd76f6f.gif

0dbobokk7lwzjkd69j9d-wmotks.gif
Так ведет себя метод со сферическим доверительным регионом.
l957xbscmthhqno0yneburansf4.gif
Так тот же метод ведет себя в том случае, если форма доверительного региона задается матрицей, построенной по правилу Давидона-Флетчера-Пауэла. Влияние на сходимость имеется, но куда скромнее, чем в аналогичном случае при использовании линейной модели целевой функции.
rnn2uzvfkpd7fdrxa26qvdlkcfs.gif
А это уже поведение метода, примененного к задаче наименьших квадратов. Сходится за 5 итераций. Только пожалуйста, не делайте из этого вывода, что вторая формулировка для функций такого вида всегда лучше первой. Это не так, просто в этом конкретном случае случае так вышло.

Заключение

Метод Левенберга-Марквардта является, на сколько мне известно, первым методом, основанным на идее доверительного региона. Он весьма неплохо показал себя на практике при решении задачи наименьших квадратов. Сходится метод в большинстве случаев (виденных мной) довольно быстро (о том, хорошо это или плохо я говорил в предыдущей статье). Однако при минимизации функций общего вида выбирать сферу в качестве доверительного региона — вряд ли наилучший из возможных вариантов. Кроме того, существенным недостатком метода (в его базовой формулировке, которая здесь и была описана) является то, что размер доверительного региона задается неявно. Недостаток проявляется в том, что зная значение 99d394e7d0b74248114405067e0ffd51.gif мы, конечно, можем в текущей точке посчитать a2b0686adbc103ad9f96be85cca5d418.gif просто вычислив длину шага 4b28c13d5f5d658adb7478fbc9efc923.gif. Однако когда мы перейдем в новую точку, этому же значению 99d394e7d0b74248114405067e0ffd51.gif будет уже соответствовать совсем другая величина доверительного региона. Таким образом, мы лишаемся возможности определения «характерной для задачи» величины доверительного региона и вынуждены в каждой новой точке определять его размер по-новому. Это может быть существенным, когда для сходимости требуется достаточно большое количество итераций, а вычисление значения функции обходится недешево. Подобные проблемы решаются в более продвинутых методах, основанных на идее доверительного региона.

Но это уже совсем другая история.

© Habrahabr.ru