При этом для любого из рассмотренных выше способов выбора

При К = 1, т.е. при иобновлении“ матрицы Гессе на каждой итерации, эта оценка равносильна (5.20), а при довольно редком ?5обновлении“, т.е. при больших значения АГ, она приближается к (5.22).

5.5. Квазиньютоновские методы

Среди алгоритмов многомерной минимизации следует вы­ делить группу алгоритмов, которые объединяют достоинства

метода наискорейш его спуска и мет ода Ньютона. Такие ал­

горитмы принято относить к так называемым к в а з и н ъ ю т о - н о в ск и м м ет о д а м . Особенность этих алгоритмов состоит в том, что при их применении нет необходимости вычислять и обращать матрицу Гессе целевой функции f{x ) и в то же время удается сохранить высокую скорость сходимости алгоритмов, присущую методу Ньютона и его модификациям (см. 5.3 и 5.4).

Элементы релаксационной последоват ельност и {х к} в ал­ горитмах квазиньютоновских методов минимизации непрерыв­ но дифференцируемой в Шп целевой функции f(x ) строят в соответствии с рекуррентным соотношением х к= х к~ 1 + щ рк, но направление спуска на каждой /с-й итерации задают в виде

Здесь w k = — grad f{x k~l) — антиградиент целевой функции в точке х к~1, &Ак — положительно определенная матрица поряд­ ка п, обновляемая на к - й итерации. Отметим, что направление, задаваемое на каждой /с-й итерации вектором р к (5.23), явля­ ется направлением спуска, так как с учетом (5.23)

(grad/( * fc_1),pfc) = - ( w k,A kwk) < 0.

Матрицы {Ak} определяют таким образом, чтобы последо­ вательность {Ak} при к —Уоо имела предел, равный i f - 1 (ж*), где Н(х*) — матрица Гессе целевой функции, вычисленная в точке минимума х* этой функции.

На к-й итерации алгоритма поиска точки минимума про­ исходит спуск из точки х к~х с шагом спуска щ\рк\, причем значение щ выбирают путем исчерпывающего спуска в напра­ влении вектора рк. На первой итерации (к = 1) спуск начинают из выбранной начальной точки х° и при этом обычно в ка­ честве А\ берут единичную матрицу 1п порядка п. Отметим, что в этом случае р1 = ги1 = — grad/(x°), т.е. первая итерация алгоритма квазиньютоновского метода полностью совпадает с итерацией метода наискорейшего спуска (см. 5.1). Если при­ нять А^ = H~l (xk~l) в (5.23) и щ = 1 в (5.14), то придем к „классическому41 методу Ньютона (см. 5.3).

Как было отмечено выше (см. 5.3), если целевая функция является сильно выпуклой, то алгоритмы метода Ньютона обла­ дают квадратичной скоростью сходимости. Поэтому можно ожидать, что алгоритмы квазиньютоновских методов будут иметь достаточно высокую скорость сходимости, если на ка­ ждой /e-й итерации матрица Ak выбрана близкой к матрице Н~1 (хк~г) в точке х к~ 1 Е W1. Используя при конструировании матрицы Ak аппроксимацию матрицы Дг_1 (жА:“ 1) с учетом ин­ формации о градиенте целевой функции в той же точке х к~! , можно существенно упростить процедуру нахождения напра­ вления спуска на fc-й итерации. Именно эти соображения лежат в основе построения алгоритмов квазиньютоновских методов.

Остановимся на вопросе построения последовательности матриц Ak. Эти матрицы строят в соответствии с рекуррент­

где A Ak — положительно определенная матрица порядка п, на­ зываемая поправочной. Способ ее выбора определяет вариант

алгоритма квазиньютоновского метода. На первой итерации, как уже отмечалось, обычно принимают А\ = 1п.

Предположим, что f(x ) = | (Qx, х) + (с, х) — квадратичная функция с положительно определенной матрицей Q. В этом слу­ чае функция f(x ) сильно выпукла, grad/(cc) = Qx + c (см. 3.6),

Дwk = grad f (x k~l) —grad/(icA:) = Q(xk~l —x k),

имы можем записать

где матрица Q совпадает с матрицей Гессе Н(х) функции f(x). В общем случае, проводя аналогию с (5.25), потребуем,

а поправочную матрицу Дл4* в (5.24) выберем, исходя из условия

т.е. так, чтобы последовательность {Ак} имела при к —>оо предел, равный Н~1 (х*).

Подставляя (5.24) в (5.26), получаем &AkAwk = —А хк —

—AkAwk. Непосредственной подстановкой можно убедиться, что для любого к£ N одним из решений этого уравнения будет

где у, z G К" — произвольные элементы из R” . Полагая у = А хк и z = Akwk, получаем с учетом (5.24) формулу

для выбора элементов последовательности {Ак}. Такой способ выбора матриц Акприводит к одному из алгоритмов квазиньютоновского типа, известного в литературе как метод Давидона — Флетчера — Пауэлла, или ДФП-метод.

Рассмотренный способ побновления41 матрицы Аь+i сохра­ няет ее свойство положительной определенности на каждой к-й итерации. Для того чтобы убедиться в этом, воспользуемся ме­ тодом математической индукции. Очевидно, что матрица А\ положительно определенная, если принять А\ = 1 п. Предполо­ жим, что при к > 1 матрица Ак также положительно опреде­ ленная, и докажем, что матрица Ak+1 в (5.27) положительно определенная.

Выберем ненулевой вектор у Е W1. Тогда с учетом (5.27) можно записать

Из курса линейной алгебры известно [IV], что любая положи­ тельно определенная матрица Ак порядка п может быть пред­ ставлена в виде Ак = СкС^, где Ск — невырожденная нижняя треугольная матрица порядка п. Обозначив и = С ^ у и v — = C k A w k, перепишем (5.28) в виде

Покажем, что знаменатель второго слагаемого в правой части (5.29) отрицателен. В самом деле, учитывая, что А хк = = щ рк, а направление спуска рк на к-й итерации определяется

согласно (5.23), получаем

(Awk, A x k) = (grad/(®fc_1) - grad f ( x k), x k - x k_1) =

= -X k(gradf(xk~1),A kgradf(xk~1)) - x k(gradf(xk),p k) < < - x fc(grad/(a?fc),p fe),

поскольку по предположению индукции матрица Ak положи­ тельно определенная. Так как точку х к находим при помощи исчерпывающего спуска из точки х к~ 1 в направлении векто­ ра рк, то в соответствии со свойством (4.57) исчерпывающего спуска правая часть последнего соотношения равна нулю, т.е. (ДкД А хк) < 0.

Всилу неравенства Коши — Буняковского имеем |ii|2|t;|2 ^

^(u, v)2, причем равенство возможно лишь при и = v, что с учетом невырожденности матрицы Ск эквивалентно равенству

у= Awk. Но тогда (А х к, у ) = (Аж*, Дги^) < 0 и в итоге при любом ненулевом векторе у £ Шп из (5.29) следует, что

т.е. матрица Ак+ 1 положительно определенная.

К соотношению (5.23) можно прийти и из других соображе­ ний. Пусть А —: положительно определенная матрица порядка п. Введем в Шп наряду со стандартным скалярным произведе­

= (A~1 gradf{x),h)A + o{\h\)

заключаем, что в евклидовом пространстве Rn со скалярным произведением (5.30) градиент целевой функции принимает вид gradAf(x ) — A~l grad/(x). Отсюда следует, что алгоритм квазиньютоновского метода, в котором используется (5.23), можно

рассматривать как алгоритм наискорейшего спуска в меняю­ щейся на каждой А;-й итерации метрике р(х,у) = (А^1 х, у) при условии, что Ак — положительно определенная матрица. По­ этому квазиньютоновские методы иногда называют методами переменной метрики*.

Можно показать**, что если

(wk,p k) f _—

**(<3р*,р*)’

то при минимизации квадратичной функции с положительно определенной матрицей Q алгоритм ДФП-метода позволяет по­ лучить точку минимума не более чем за п итераций. Векторы рк, к = 1 ,п, определяющие направления спуска, являются сопря­ женными относительно матрицы Q и связаны сооношениями

матрица AkQ имеет п собственных значений, равных единице, а такая матрица является единичной, т.е. AnQ = 1 п. Отсюда Ап = Q-1 , и, следовательно, матрица Ап является обратной к матрице Гессе Н(х*) = Q квадратичной функции, вычисленной в точке минимума х*

Отметим, что если в алгоритме ДФП-метода А\ = / п, то для квадратичной функции траектория поиска совпадает с тра­ екторией, полученной по методу сопряженных направлений.

Например, при минимизации квадратичной функции, рассмо­ тренной в примерах 5.1 и 5.3, алгоритм ДФП-метода приводит к точке х* = (—л/5) —2\/5) минимума этой функции из началь­ ной точки ( - 2, 1 ) за две итерации, причем траектория поиска при выборе А\ = I2 полностью совпадает с полученной в при­ мере 5.3 методом сопряженных направлений (см. рис. 5.4).

При минимизации целевой функции, не являющейся квадра­ тичной, использование ДФП-метода в общем случае не позволя­

*См.: Поляк Б . Т.

**См.: А оки М., а также: Базара М Ш е т т и К.

ет найти точку минимума за конечное число итераций. Как и в методе сопряженных направлений, чтобы ослабить влияние накапливаемых погрешностей на сходимость ДФП-метода и уменьшить вероятность появления после очередных п итераций линейно зависимых направлений спуска, применяют процеду­

ничную матрицу, т.е. полагают Amn+\ = А\ = In, m € N. Опишем алгоритм ДФП-метода в случае, когда (не ква­

дратичная) целевая функция f(x ) дифференцируема в Rn. На предварительном этапе задаем значение £3 параметра точно­ сти поиска в условии прекращения итераций |го*| < £3 и вы­ бираем начальную точку х°. Принимаем А\ = / п, формируем множество 1 о = {п, 2п, ...} моментов обновления алгоритма,

полагаем k = 1 и переходим к основной части алгоритма.

1.На к-й итерации в точке ж*- 1 вычисляем антиградиент wk = — grad/(a?fc-1) целевой функции и проверяем выполнение неравенства |iufc|< £3. Если оно выполнено, то итерации пре­ кращаем и полагаем х* « ж*-1, /(ж*) и /(ж *-1). В противном случае переходим к п. 2.

2.Используя (5.23), вычисляем вектор рк = AkWk, опреде­ ляющий направление спуска из точки х к~1, и, минимизируя

функцию фк{х) = /(ж *- 1 + x Pk)i находим значение х* и точку хк = х к~ 1 + ХкРк. Если к € Jo, то принимаем Afc+i = 1п, по­ лагаем к := к + 1 и возвращаемся к п. 1 . В противном случае переходим к п. 3.

3. Полагаем А хк = х к —х к~к, Awk = wk+1 —wk, по формуле (5.27) вычисляем матрицу Ak+i,. полагаем к := к + 1 и возвра­ щаемся к п. 1 .

Как отмечено выше (см. 5.2), применение процедуры „обно­ вления" алгоритма увеличивает общее число итераций, необхо­ димых для достижения заданной точности поиска. Поэтому на практике используют варианты ДФП-метода, в которых (как и в случае метода сопряженных направлений) не используется

Таблица 5.8

Из результатов вычислений видно, что заданная точность поиска минимума достигнута за 7 итераций, матрицы Ai и Н~1 (х7) на конечном, седьмом шаге поиска близки друг к другу. Матрица Н~1 (х*)) вычисленная в точке минимума ж*, имеет вид

0,5 1,0

Н -\х*)

1,0 2,5

и также мало отличается от А7 и Я - 1 (ж7). Отметим, что для достижения той же точности в методе наискорейшего спуска потребовалось 96 итераций. #

Возможны и другие способы построения последовательно­ сти {Ak} матриц с применением (5.24). Так, например, в алго­ ритме квазиньютоновского метода, называемого в литературе методом Бройдена — Флетчера — Шенно, или БФШ-методом,

используют соотношение

в алгоритме метода Пауэлла —

где Ахк = А хк + АкAwk, к G N, а в алгоритме метода Мак-Кор- мика —

причем, как и в случае ДФП-метода, обычно полагают А\ = /„. Можно показать, что любой из рассмотренных способов „ обновления “ матрицы Ак+\ сохраняет ее свойство положитель­

ной определенности, а последовательность {-A*,} при к оо сходится к Н~1(х*).

Пример 5.9. Сравним различные алгоритмы на примере минимизации функции f(x 1,2:2) = ixi — х2) 2 + ix i — I)2, рассмо­ тренной в ряде предыдущих примеров. Выберем начальную точку а:0= (—1 , -2 ), в которой f(x°) = 13 и параметр точности поиска ез = 10~3. Как и в примере 5.8, в применяемых алгорит­ мах не будем использовать „обновление" через фиксированное число шагов.

В табл. 5.9 приведены координаты точек х к, найденных при помощи трех алгоритмов кваэиньютоновских методов. Траек­ тории поиска точки минимума для этих алгоритмов показаны на рис. 5.11 (о — БФШ-метод; б — метод Пауэлла; в — метод Мак-Кормика).

Таблица 5.9

Сравнительный анализ результатов показывает, что для рассматриваемой функции наилучший результат по количеству итераций, требуемых для достижения заданной точности, дают ДФП-метод и БФШ-метод. Метод Мак-Кормика по этому кри­ терию эффективности им заметно уступает. #

Отметим, что сравнение эффективности алгоритмов мини, мизации принято проводить на специально сконструированных функциях. Графики этих функций имеют четко выраженную овражную структуру. На рис. 5.12 представлены линии уровня унимодальной функции Розенброка

а на рис. 5.13 — функции Химмельблау