Levenberg-Marquardt method for unconstrained optimization

Full Text

Введение В этой работе исследуются численные методы решения задачи безусловной оптими- зации f(x) ! min; x 2 Rn; (1) с по крайней мере дважды дифференцируемой целевой функцией f : Rn ! R, и при условиях допускающих, что задача (1) может иметь неизолированные решения. 62 А. Ф. Измаилов, А. С. Куренной, П. И. Стецюк Один успешный подход к решению нелинейных задач с возможно неизолированны- ми решениями основан на методе Левенберга-Марквардта [1], [2], который, однако, в принципе предназначен для решения систем нелинейных уравнений, а не задач опти- мизации. Для уравнения f0(x) = 0; (2) описывающего стационарные точки задачи (1), и для текущего приближения xk 2 Rn , этот метод определяет следующее приближение как xk+1 = xk + pk , где pk является единственным решением линейного уравнения (H2 k + kI)p = Hkf0(xk); (3) в котором базовым выбором симметричной n n матрицы Hk является f00(xk) , а k > 0 играет роль параметра регуляризации. Из результатов в [3], [4] следует, что при правильном управлении этим параметром метод Левенберга-Марквардта облада- ет локальной квадратичной сходимостью к стационарной точке задачи (1) при очень слабых предположениях, допускающих неизолированность стационарных точек. Более того, сходимость метода может быть глобализована одномерным поиском для квадрата невязки уравнения (2) в роли функции качества. Однако, в данном (оптимизационном) контексте указанный способ глобализации сходимости метода Левенберга-Марквардта не вполне удовлетворителен, поскольку по- лучаемый глобализованный алгоритм направлен на поиск стационарных точек задачи (1), а не ее решений. В частности, такой алгоритм ѕне различаетї максимумы и мини- мумы функции f , как и любые стационарные точки задачи (1). В связи с этим обсто- ятельством в настоящей работе разрабатывается ѕболее оптимизационныйї алгоритм, использующий направления pk метода Левенберга-Марквардта, но с одномерным по- иском для целевой функции f задачи (1) вместо невязки уравнения (2). 1. Глобализованный алгоритм Ключевое наблюдение заключается в следующем. Направление pk метода Левен- берга-Марквардта можно явно выразить из (3) как pk = Qkf0(xk); (4) где Qk = (H2 k + kI) 1Hk: Предположим, что матрица Hk положительно определена; тогда матрица Qk невы- рождена (как произведение двух невырожденных матриц), и Q 1 k = H 1 k (H2 k + kI) = Hk + kH 1 k : (5) Но полученная таким образом матрица является симметричной и положительно опре- деленной, откуда следует, что таковой является и матрица Qk , а значит, pk является направлением убывания функции f в точке xk , если f0(xk) 6= 0 . МЕТОД ЛЕВЕНБЕРГА-МАРКВАРДТА ДЛЯ ЗАДАЧ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ 63 Разумеется, матрица Hk = f00(xk) может не быть положительно определенной, даже при xk сколь угодно близком к множеству решений. Поэтому глобализованный алго- ритм должен включать в себя механизм выбора подходящего Hk , причем допускающий базовый выбор вблизи решений. А л г о р и т м. Выбираем параметры 1; 2 > 0 , 1 > 0 , 2 > 1 , > 0 , q > 0 , "; 2 (0; 1) . Выбираем x0 2 Rn и полагаем k = 0 . 1. Если f0(xk) = 0 , стоп. 2. Полагаем Hk = f00(xk) . Если kHkf0(xk)k 1kf0(xk)k1 ; (6) вычисляем pk как решение уравнения (3) при k = minf; kf0(xk)kqg . Иначе переходим к шагу 4. 3. Если hf0(xk); pki 2kpkk2 ; (7) переходим к шагу 5. 4. Вычисляем симметричную n n матрицу Hk и соответствующее решение pk задачи (3), удовлетворяющие (6) и (7). 5. Вычисляем k = j , где j наименьшее неотрицательное целое число, при ко- тором выполняется неравенство Армихо f(xk + jpk) f(xk) + "jhf0(xk); pki: (8) 6. Полагаем xk+1 = xk + kpk , увеличиваем k на 1 , и переходим к шагу 1. Ключевой вопрос состоит в том как реализовать шаг 4. Очевидно, что (6) может выполняться для любого 1 > 0 , если Hk достаточно положительно определена, а именно, ее минимальное собственное значение не меньше 1 . Более того, из (4) и (5) сле- дует, что (7) также достигается достаточной положительной определенностью Hk . По- следнего же можно добиться посредством модифицированного разложения Холецкого [5, разд. 4.4.2.2], [6], или посредством модифицированного симметричного знаконеопре- деленного разложения [7] матрицы Гессе f00(xk) ; см. также [8, разд. 3.4]. При этом на каждой итерации алгоритма потребуется решить не более двух систем линейных урав- нений, как и в случае, когда для второй системы берется просто Hk = I с достаточно большим > 0 . Альтернативным образом можно выбирать Hk последовательно, заменяя Hk на Hk + !I с некоторым ! > 0 , до тех пор, пока (6) и (7) не окажутся выполнены. Разумеется, при таком подходе итерация алгоритма может потребовать решения более двух систем линейных уравнений, но получаемое в результате направление может быть лучше, так как выбранная в итоге матрица Hk может быть ближе к истинной матрице Гессе функции f . 64 А. Ф. Измаилов, А. С. Куренной, П. И. Стецюк 2. Глобальная сходимость Следующая теорема описывает свойства глобальной сходимости предложенного ал- горитма. Теорема 1. Пусть f дважды дифференцируема на Rn . Тогда алгоритм корректно определен, и либо останавливается после конечного числа итераций на шаге 1 в некоторой стационарной точке xk задачи (1), либо ге- нерирует такую бесконечную последовательность fxkg , что, в случае ограниченно- сти fHkg , любая предельная точка этой последовательности является стационар- ной точкой задачи (1). Д о к а з а т е л ь с т в о. Корректная определенность алгоритма вытекает из его конструкции. Допустим, что алгоритм не останавливается на шаге 1, а значит, генерирует беско- нечную последовательность fxkg . Предположим далее, что последовательность fHkg ограничена. Покажем, что последовательность fpkg направлений поиска является рав- номерно градиентной в терминологии [9, c. 24], т. е. если некоторая подпоследователь- ность fxkjg сходится к x 2 Rn , причем f0(x) 6= 0 , то последовательность fpkjg огра- ничена и lim sup j!1 hf0(xkj ); pkj i < 0: (9) Из (7) получаем, что для всех k kf0(xk)kkpkk 2kpkk2: Отсюда и из непрерывности f0 , которая автоматически следует из двукратной диффе- ренцируемости f , вытекает ограниченность последовательности fpkjg . Кроме того, в силу (3) и (6), для всех k справедливо соотношение k(H2 k + kI)pkk = kHkf0(xk)k 1kf0(xk)k: Из ограниченности fHkg и fkjg (где последнее следует из определения k и, опять же, непрерывности f0 ) вытекает, что если fpkjg имеет 0 в качестве предельной точки, то f0(x) = 0 , что противоречит сделанному предположению. Следовательно, норма pkj отделена от нуля, и тогда из (7) вытекает (9). Таким образом, fpkg является равномерно градиентной последовательностью на- правлений поиска. А значит, согласно [9, теорема 1.8], каждая предельная точка fxkg является стационарной точкой задачи (1). 3. Квадратичная скорость сходимости Обозначим через S множество стационарных точек задачи (1): S = fx 2 Rn j f0(x) = 0g: МЕТОД ЛЕВЕНБЕРГА-МАРКВАРДТА ДЛЯ ЗАДАЧ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ 65 Лемма 1. Пусть f дважды непрерывно дифференцируема в окрестности стаци- онарной точки x задачи (1), и пусть выполняется следующая локальная липшицева оценка расстояния до множества стационарных точек: dist(x; S) = O(kf0(x)k) (10) при x ! x . Тогда существует такое > 0 , что kf00(x)f0(x)k kf0(x)k выполняется для любого x 2 Rn; достаточно близкого к x . Д о к а з а т е л ь с т в о. Согласно [10, следствие 2], в сделанных предположениях выполняется dist(x; S) = O(kf00(x)f0(x)k) при x ! x . Нужное утверждение получается комбинированием данной оценки с оцен- кой, следующей из непрерывности f00 : kf0(x)k = O(dist(x; S)) при x ! x . Таким образом, в предположениях этой леммы, если в алгоритме либо 1 > 1 , либо 1 = 1 и 1 > 0 достаточно мало, то при Hk = f00(xk) (6) выполняется для всех xk 2 Rn достаточно близких к x . Во избежание ненужных сложностей, в утверждениях ниже будет удобно формально предполагать, что если xk 2 S (и, соответственно, k = 0 ), то pk = 0 . Заметим, что если xk 2 S , то алгоритм останавливается на шаге 1, и вычисления pk на самом деле больше не выполняются. С другой стороны, pk = 0 всегда является решением (3) при xk 2 S , и для целей анализа удобно предполагать, что в таких случаях выбирается именно это решение. Лемма 2. Пусть в дополнение к предположениям леммы 1 вторая производная f удовлетворяет условию Липшица в окрестности x . Тогда для любого q 2 [1; 2] существует такое > 0 , что для решения pk уравне- ния (3) при Hk = f00(xk) выполнено kf00(x)pkk kpkk для всех xk 2 Rn; достаточно близких к x . Д о к а з а т е л ь с т в о. Согласно [4, леммы 2.1, 2.2], в сделанных предположениях выполняются оценки kpkk = O(dist(xk; S)); (11) 66 А. Ф. Измаилов, А. С. Куренной, П. И. Стецюк dist(xk + pk; S) = O (dist(xk; S))(2+q)=2 (12) при xk ! x . Точнее, в [4, лемма 2.2] дополнительно предполагается, что xk+pk остается достаточно близким к x . Однако, в доказательстве в [4, теорема 2.1] установлено, что любая необходимая близость xk + pk к x гарантируется достаточной близостью xk к x , а значит, это предположение можно опустить. От противного: предположим, что существует последовательность fxkg Rn , схо- дящаяся к x и такая, что f00(x)pk = o(kpkk) при k ! 1. Тогда согласно (11) f00(x)pk = o(dist(xk; S)); (13) в то время как согласно (12) f0(xk + pk) = O(dist(xk + pk; S)) = o(dist(xk; S)): (14) С помощью теоремы о среднем, вновь привлекая (11), получаем f0(xk + pk) f0(xk) f00(x)pk = o(kpkk) = o(dist(xk; S)); а значит, принимая во внимание (13) и (14), f0(xk) = o(dist(xk; S)) при k ! 1, что противоречит (10). Лемма 3. Пусть f дважды дифференцируема в окрестности локального реше- ния x задачи (1), а ее вторая производная удовлетворяет условию Липшица в этой окрестности. Допустим, что выполняется локальная липшицева оценка расстояния (10) при x ! x . Тогда для любого q 2 [1; 2] существует такое > 0 , что для решения pk уравне- ния (3) при Hk = f00(xk) выполнено hf00(x)pk; pki kpkk2 для всех xk 2 Rn , достаточно близких к x . Д о к а з а т е л ь с т в о. В силу необходимого условия оптимальности второго по- рядка f00(x) является неотрицательно определенной матрицей. Но тогда hf00(x); i > 0 для всех 2 Rn , при которых f00(x) 6= 0 . Требуемый результат теперь получается из леммы 2 с привлечением свойств однородности и компактности. Следующий результат вытекает из (12) и леммы 3 посредством рассуждения в [11, следствие 5.1]. МЕТОД ЛЕВЕНБЕРГА-МАРКВАРДТА ДЛЯ ЗАДАЧ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ 67 Лемма 4. В предположениях леммы 3, для каждого q 2 [1; 2] , если либо 2 > 2 , либо 2 = 2 и 2 > 0 достаточно мало, то (7) выполняется для решения pk уравнения (3) с Hk = f00(xk) при любом xk 2 Rn , достаточно близком к x . Более того, если " 2 (0; 1=2) , то (8) выполняется при j = 0 , т. е. на шаге 5 алгоритма принимается k = 1 . Теорема 2. Пусть в предположениях леммы 3 последовательность fxkg полу- чена алгоритмом при q 2 [1; 2] , " 2 (0; 1=2) и либо при 1 > 1 , либо при 1 = 1 и достаточно малом 1 > 0 , и либо при 2 > 2 , либо при 2 = 2 и достаточно малом 2 > 0 . Предположим, что при некотором k приближение xk оказывается достаточно близким к x . Тогда последовательность fxkg сходится к некоторой стационарной точке задачи (1), и для любого достаточно большого k выполняется Hk = f00(xk) , k = 1 , причем скорость сходимости квадратичная. Д о к а з а т е л ь с т в о. Это следует из лемм 1 и 4, а также из [4, теорема 2.2], где нужно принять во внимание, что, согласно рассуждению в [4, теорема 2.1], доста- точная близость xk к x обеспечивает нужную близость всех последующих приближе- ний к x . Покажем теперь, что оценка расстояния (10) обеспечивается выполнением условия квадратичного роста в точке x , что означает существование такого > 0 , что f(x) f(bx) (dist(x; S))2 (15) для любого x 2 Rn , достаточно близкого к x , и для некоторой проекции bx точки x на S . П р е д л о ж е н и е 1. Пусть f дважды дифференцируема в окрестности ста- ционарной точки x задачи (1), а ее вторая производная непрерывна в x , и предполо- жим, что в точке x выполнено условие квадратичного роста. Тогда имеет место локальная липшицева оценка расстояния (10) при x ! x . Д о к а з а т е л ь с т в о. По теореме о среднем, и принимая во внимание, что bx ! x при x ! x , получаем оценки f0(x) f00(bx)(x bx) = f0(x) f0(bx) f00(bx)(x bx) = o(dist(x; S)); f(x) f(bx) 1 2 hf00(bx)(x bx); x bxi = f(x) f(bx) hf0(bx); x bxi 1 2 hf00(bx)(x bx); x bxi = o((dist(x; S))2) 68 А. Ф. Измаилов, А. С. Куренной, П. И. Стецюк при x ! x . Следовательно, в силу (15), (dist(x; S))2 f(x) f(bx) = 1 2 hf0(x); x bxi + o((dist(x; S))2) 1 2 kf0(x)kkx bxk + o((dist(x; S))2) = 1 2 kf0(x)k dist(x; S) + o((dist(x; S))2); откуда очевидным образом следует (10). 4. Численные примеры Помимо глобализации метода Левенберга-Марквардта посредством одномерного по- иска для квадрата евклидовой невязки уравнения (2), естественным конкурентом пред- ложенного в этой работе алгоритма является регуляризованный метод Ньютона, на- правление pk которого в текущем приближении xk вычисляется из линейного уравне- ний (Hk + kI)p = f0(xk); (16) а глобализация сходимости осуществляется посредством одномерного поиска для функ- ции f ; см. [12], [13]. Этот метод также может требовать модификации базового выбора Hk = f00(xk) для обеспечения того, чтобы pk было направлением убывания функции f в точке xk . Более того, с учетом неотрицательной определенности матрицы H2 k , есть основания ожидать, что итерация метода Левенберга-Марквардта будет требовать меньшего количества последовательных модификаций Hk , а значит, решения меньшего количества систем линейных уравнений. Помимо этого, важным теоретическим преиму- ществом метода Левенберга-Марквардта является то, что его итерационные системы (3) всегда разрешимы, поскольку их матрицы положительно определены (а значит, невырождены) при k > 0 . Таким образом, в этом разделе рассматриваются следующие алгоритмы. 1. Алгоритм из разд. 1 с последовательной модификацией матрицы Hk по формуле Hk = Hk + !I на шаге 4, в котором (3) заменено на итерационное уравнение (16) регуляризованного метода Ньютона, а тест (6) опущен, но процедура модификации Hk инициируется и в тех случаях, когда уравнение (16) решить не удается. Такой алгоритм соответствует глобализованным версиям регуляризованного метода Ньютона в [12], [13]. Варианты этого алгоритма с выбором q = 1 и q = 2 будем обозначать RNM1 и RNM2, соответственно. 2. Алгоритм из разд. 1 с последовательной модификацией матрицы Hk по той же формуле. Варианты этого алгоритма с выбором q = 1 и q = 2 будем обозначать L-MM1-obj и L-MM2-obj, соответственно. 3. Алгоритм из разд. 1, но с опущенными тестами (6) и (7), и опущенным шагом 4, a также с функцией f на шаге 5, замененной на квадрат невязки уравнения (2), т. е. МЕТОД ЛЕВЕНБЕРГА-МАРКВАРДТА ДЛЯ ЗАДАЧ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ 69 на функцию ' : R ! R, '(x) = 1 2 kf0(x)k2: Такой алгоритм соответствует глобализованным версиям метода Левенберга-Марквар- дта в [3], [4]. Варианты этого алгоритма с выбором q = 1 и q = 2 будем обозначать L-MM1-res и L-MM2-res, соответственно. В алгоритмах использовались следующие значения параметров: 1 = 2 = 10 9 , 1 = 1:1 , 2 = 2:1 , = 1 , " = 0:01 , = 0:5 , ! = 10 . Запуск считался успешным, если для некоторого k 500 реализовалось неравенство kf0(xk)k < 10 8: Если на некоторой итерации на шаге 5 алгоритма становилось меньше 10 12 , то запуск считался неудачным и прекращался. Для каждой из рассматриваемых ниже трех задач выполнялось 1000 запусков каж- дого алгоритма из случайных начальных точек, равномерно распределенных в области kxk1 100 . Значения столбцов в приводимых ниже таблицах таковы - ѕSї (от ѕSucessesї): процент успешных запусков; - ѕIї (от ѕIterationsї): среднее количество итераций на один успешный запуск; - ѕLSї (от ѕLinear Systemsї): среднее количество решенных линейных систем на один успешный запуск; - ѕOVї (от ѕObjective function Valuesї): средняя величина натуральных логариф- мов от значений целевой функции в точках завершения алгоритма (вне зависимо- сти от успешности запусков). Эта величина характеризует качество получаемых приближений: чем она меньше, тем в среднем лучше получаемые приближения. При этом она более устойчива по отношению к единичным статистическим вы- бросам, чем средняя величина самих значений целевой функции. Таблица 1: Результаты для примера 1 Method S (%) I LS OV RNM1 100 32 32 -61.81 RNM2 95 32 32 -61.57 L-MM1-obj 100 32 32 -61.47 L-MM2-obj 100 32 32 -61.64 L-MM1-res 100 32 32 -61.78 L-MM2-res 96 32 32 -59.92 70 А. Ф. Измаилов, А. С. Куренной, П. И. Стецюк П р и м е р 1. Пусть n = 2 , f(x) = ((x21 + x22 )2 2(x21 x22 ))2 . Тогда множество решений задачи (1) это лемниската Бернулли. Лемниската Бернулли содержится в круге с центром в нуле радиуса p 2 . В опреде- ленном смысле это объясняет то, что при запусках из удаленных от нуля начальных точек все рассматриваемые алгоритмы ведут себя в этом примере схожим образом; см. таблицу 1. Таблица 2: Результаты для примера 2 Method S (%) I LS OV RNM1 100 20 21 -44.60 RNM2 100 26 27 -27.36 L-MM1-obj 100 18 18 -53.29 L-MM2-obj 100 18 18 -51.81 L-MM1-res 100 18 18 -53.61 L-MM2-res 100 18 18 -53.16 П р и м е р 2. Пусть n = 2 , f(x) = x21 x22 . Тогда множество решений задачи (1) задается уравнением x1x2 = 0 . Как видно из таблицы 2, в этом примере алгоритмы, основанные на методе Левенбер- га-Марквардта, ведут себя схожим образом и превосходят RNM как по эффективности, так и по качеству получаемых приближений. При этом все алгоритмы демонстрируют абсолютную робастность. Таблица 3: Результаты для примера 3 Method S (%) I LS OV RNM1 100 28 28 -27.44 RNM2 100 27 27 -27.54 L-MM1-obj 100 17 17 -57.65 L-MM2-obj 100 19 19 -52.57 L-MM1-res 100 17 17 -57.82 L-MM2-res 99 19 19 -53.45 П р и м е р 3. Пусть n = 3 , f(x) = (x21 + x22 x23 )2 . Тогда множество решений задачи (1) задается уравнением x21 + x22 = x23 . Как видно из таблицы 3, в этом примере ситуация та же, что и в примере 2. МЕТОД ЛЕВЕНБЕРГА-МАРКВАРДТА ДЛЯ ЗАДАЧ БЕЗУСЛОВНОЙ ОПТИМИЗАЦИИ 71 Чтобы продемонстрировать различие в поведении двух рассмотренных глобализа- ций методов Левенберга-Марквардта, рассмотрим еще один простой (одномерный) при- мер. Рис. 1: Пример 4 Таблица 4: Результаты для примера 4 Method S (%) I LS OV CS (%) RNM1 77 6 6 -113.97 100 RNM2 77 5 6 -113.39 100 L-MM1-obj 80 5 6 -112.14 100 L-MM2-obj 80 5 6 -112.24 100 L-MM1-res 100 4 5 -45.72 49 L-MM2-res 100 4 5 -44.80 48 П р и м е р 4. Пусть n = 1 , f(x) = x4=2 104x2 . Тогда множество решений задачи (1) состоит из двух точек: 100 и 100 . Кроме того, функция f имеет локальный максимум в точке 0 (см. рис. 1). Помимо данных, сообщавшихся в таблицах для предыдущих примеров, таблица 4 имеет еще столбец ѕCSї (от ѕConvergences to Solutionї), в котором указан процент 72 А. Ф. Измаилов, А. С. Куренной, П. И. Стецюк успешных запусков со сходимостью к одному из решений, а именно, закончившихся в точке, значение функции f в которой отличалось от оптимального значения 108=2 не более чем на 10 5 . Из таблицы 4 видно, что алгоритмы с одномерным поиском для целевой функции f существенно превосходят алгоритмы с одномерным поиском для квадрата невязки уравнения (2) в плане качества получаемых приближений. На рис. 1 также показано распределение начальных точек, из которых имела ме- сто сходимость к одному из решений (черные точки), и тех, из которых имела место сходимость к локальному максимуму (белые точки). Распределение становится весьма сложным вблизи 100= p 3 . Данные приведены для L-MM1-res; для L-MM2-res картина аналогичная. Разумеется, приведенные численные результаты предназначены лишь для иллю- страции теории, разработанной в разд. 1-3, и, в частности, не подразумевают каких- либо далеко идущих выводов. Систематическое сравнительное численное тестирование описанных алгоритмов составит одно из направлений дальнейшей работы авторов.

About the authors

Alexey F. Izmailov

Lomonosov Moscow State University

Email: izmaf@ccas.ru
Doctor of Physics and Mathematics, Professor of the Operations Research Department VMK Faculty, Leninskiye Gory, Moscow 119991, Russian Federation

Alexey S. Kurennoy

Tambov State University named after G.R. Derzhavin

Email: akurennoy@cs.msu.ru
Candidate of Physics and Mathematics, Researcher 33 Internatsionalnaya St., Tambov 392000, Russian Federation

Petr I. Stetsyuk

V. M. Glushkov Institute of Cybernetics of NAS of Ukraine

Email: stetsyukp@gmail.com
Doctor of Physics and Mathematics, Head of the Nonsmooth Optimization Methods Department 40 Akademika Glushkova Ave., Kiev 03187, Ukraine

References

Supplementary files