Загребаев Методы матпрограммирования 2007

λ* −

k 1

где	xk = xk −1 + λ*k −1Sk −1 ,	(3.12)
λ*k −1	= arg(min f (xk −1 + λk −1Sk −1 )) .
	λk −1

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

Таким образом, метод сопряженных градиентов сходится для квадратичных функций за число шагов, не превышающих n.

Для того чтобы распространить метод сопряженных градиентов на случай поиска экстремума неквадратичных функций (в том числе заданных неявно), рассмотрим способ выбора направления спуска, не использующий в явном виде матрицу Гессе Q .

Для этого воспользуемся выражением градиента квадратичной функции:

f (xr) = Qx + b .

Тогда, принимая во внимание итерационное соотношение (3.12), можно записать,

f (xk ) − f (xk −1 ) = Q(xk − xk −1 ) = λ*k −1QSk −1 .

(3.13)

Откуда получаем

QSk −1 = f (xk ) − f (xk −1 ) .

Подставим полученное выражение в формулу (3.11) для βk −1 :

βk −1 = < f (xk ), f (xk ) − f (xk −1 ) > ,

<Sk −1 , f (xk ) − f (xk −1 ) >

или, раскрывая скалярные произведения и производя несложные операции, получим

rr 644471 4448

βk −1 = < f (xk ), f (xk ) > −< f (xk ), f (xk −1 ) > .

<Sk −1, f (xk ) > − < Sk −1, f (xk −1 ) >

1442443 144424443

2

3

211

Рассмотрим каждый из занумерованных членов.

2-й член. Так как xk – точка минимума на направлении Sk −1 , то

градиент в этой точке ортогонален к направлению поиска (рис. 3.26). Откуда получаем:

< Sk −1 , f (xk ) >= 0 .

f (xk )

Sk −1

	xk	Рис. 3.26. Ортогональ-
		ность градиента f ( xk )
		и направления S k −1

3-й член. Подставим вместо Sk −1 его выражение, определяемое формулой (3.10). Тогда получим:

< Sk −1 , f (xk −1 ) >=< − f (xk −1 ), f (xk −1 ) > +βk −2 < Sk −2 , f (xk −1 ) > .

Очевидно, второе слагаемое равно нулю, как скалярное произведение двух ортогональных векторов. Таким образом,

< Sk −1 , f (xk −1 ) >= − < f (xk −1 ), f (xk −1 ) > .

1-й член. Выразим f (xk −1 ) в соответствии с формулой (3.10), тогда получим

< f (xk ), f (xk −1 ) >=< f (xk ),−Sk −1 + βk −2 Sk −2 >= = − < f (xk ), Sk −1 > +βk −2 < f (xk ), Sk −2 > .

Разрешим (3.13) относительно f (xk ) и подставим в последнее выражение. Тогда получим:

212

< f (x ), f (x − ) >= β − < f (x − ) + λ* − QS − , S − >=

k k 1 k 2 k 1 k 1 k 1 k 2

= β − < f (x − ), S − > +β − λ* − < QS − , S − >.

k 2 k 1 k 2 k 2 k 1 k 1 k 2

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

< f (xk ), f (xk −1 ) >= 0 .

Это означает, что вектора f (xk ), f (xk −1 ) ортогональны. Вообще можно показать, что все градиенты f (x0 ), f (x1 ), ..., f (xk )

взаимно ортогональны.

Итак, в окончательном виде имеем:

βk −1

=

< f (xk ), f (xk ) >

=

f (xk )

2

.

< f (xk −1 ), f (xk −1 )

>

f (xk −1 )

2

Последнюю формулу можно применять для неквадратичных функций, так как в нее не входит матрица Q .

Метод сопряженных градиентов является одним из наиболее эффективных методов минимизации достаточно гладких функций. При определенных условиях метод обладает квадратичной скоростью сходимости:

xk −1 − x* ≤ c xrk − x* 2 , c > 0 .

3.3.3.Методы второго порядка

3.3.3.1.Метод Ньютона (метод Ньютона – Рафсона)

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

213

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

Метод Ньютона является прямым обобщением метода отыскания корня уравнения ϕ(x) = 0 , где ϕ(x) – функция скалярной пе-

ременной. Разложение в ряд Тейлора до первого члена позволяет переписать уравнение в следующем виде:

0 = ϕ(xk +1 ) ≈ ϕ(xk ) + ϕ′(xk )(xk +1 − xk ) .

Тогда при определенных условиях можно улучшить приближение к значению корня следующим образом:

xk +1 = xk − ϕ′(xk ) .

ϕ(xk )

Нас интересует n-мерная задача оптимизации, сводящаяся фактически к определению корня уравнения f (x) = 0 . Разложение в ряд Тейлора в этом случае дает:

0 = f (xk +1 ) ≈ f (xk ) + Q(xk )(xk +1 − xk ) ,

где Q(xk ) – матрица Гессе.

Отсюда

xk +1 = xk − Q −1 (xk ) f (xk )

при условии, что существует обратная матрица Q −1 (xk ) . Эта фор-

мула и определяет метод Ньютона минимизации функции. Если функция квадратичная

f (x) = a + x тb + 12 x T Qx ,

где Q – положительно определенная матрица, то, исходя из произвольной начальной точки x0 , с помощью полученной выше формулы можно получить следующую точку:

x1 = x0 − Q−1 (b + Qx0 ) = −Q−1b .

Эта точка является точкой минимума квадратичной функции.

214

Таким образом, для квадратичной функции метод Ньютона определяет точку минимума за один шаг.

3.3.3.2. Сходимость метода Ньютона

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

Q −1 (xk ) f (xk ) , а не уменьшаться. Кроме того, матрица Гессе на

протяжении всего процесса должна быть невырожденной, так как необходимо существование обратной матрицы.

Если матрица удовлетворяет этим условиям и, кроме того, является ограниченной и удовлетворяет условию Липшица, то сущест-

вует некоторая окрестность точки минимума x* , такая, что для любой начальной точки x0 из этой окрестности метод Ньютона сходится с квадратичной скоростью:

xk −1 − x* ≤ c xk − x* 2 , c ≥ 0 .

Таким образом, для сходимости метода начальная точка x0 должна выбираться достаточно близко к искомой точке минимума x* .

Подводя итог, можно следующим образом сформулировать достоинства и недостатки метода.

Недостатки:

необходимость вычислять и обращать матрицу вторых производных. В ряде задач трудоемкость итерации метода Ньютона может за счет этого оказаться непомерно большой;

сходимость метода зависит от выбора начальной точки x0 . В

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

Достоинства:

вблизи точки минимума метод обеспечивает более быструю сходимость, чем градиентные методы;

215

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

Ниже приводятся некоторые модификации метода Ньютона, направленные на устранение указанных выше недостатков.

3.3.3.3.Метод Ньютона с регулировкой шага

Внекоторых случаях можно управлять сходимостью метода Ньютона, изменяя размеры шагов, т.е. определяя последовательность точек по формуле

xk +1 = xk − ρk Q −1 (xk ) f (xk ) ,

где 0 < ρk ≤1 (обычному методу Ньютона соответствует ρk =1 ).

В этом случае метод называется методом Ньютона с регулировкой шага, или обобщенным методом Ньютона.

Существуют два способа выбора параметра ρk . 1. Первый способ заключается в следующем:

1) выбирается ρk =1 и вычисляется точка x = xk + pk , где pk = −Q−1 (xk ) f (xk ) ;

2)вычисляется f (x) = f (xk + pk ) ;

3)производится проверка неравенства

f (x) − f (xk ) ≤ (ερk < f (xk ), pk >) , 0 < ε < 12 ,

4) если это неравенство выполняется, то в качестве искомого значения ρk берем ρk =1 . В противном случае производится

дробление ρk до тех пор, пока неравенство не выполнится.

2. В другом варианте метода, значение ρk выбирается из условия минимума функции в направлении движения:

f (xk − ρk Q−1 (xk ) f (xk )) = min f {xk − ρQ−1 (xk ) f (xk )} .

ρ≥0

216

Преимущество данной модификации метода Ньютона состоит в том, что метод сходится при любом выборе начального приближения x0 . Скорость сходимости при этом остается квадратичной,

трудоемкость каждой итерации увеличивается. Однако при первом способе выбора параметра ρk трудоемкость каждой итерации не-

сколько увеличивается.

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

Замечания. Как было отмечено в п. 3.3.3.2, на сходимость метода Ньютона большое влияние оказывает матрица Гессе. Если она не является положительно определенной, то метод расходится. Если она положительно определенная (т.е. выпуклая f (x) ), то метод сходится в окрестности точки минимума.

Определить, является ли матрица положительно определенной можно либо используя критерий Сильвестра, либо собственные числе матрицы Гессе. Матрица будет положительно определенная, если все угловые миноры матрицы положительны (критерий Сильвестра), или все собственные числа матрицы положительны.

В качестве примера рассмотрим следующие задачи. Задача 1. Рассмотрим функцию

f (x) = 2x12 + x22 + sin(x1 + x2 ) .

Проверить является ли для данной функции матрица Гессе положительно определенной?

Решение. Градиент функции имеет вид:

f (x) = 4x1 + cos(x1 + x2 ) . 2x2 + cos(x1 + x2 )

Матрица Гессе:

Q(x) =		4 − sin(x1 + x2 )				− sin(x1 + x2 )				.
		− sin(x	+ x	2	)	2 − sin(x	+ x	2	)
		1				1

217

Угловыми минорами будут:

∆1 = 4 −sin(x1 + x2 ) и ∆2 = 8 − 6 sin(x1 + x2 ) > 0 .

Так как sin(x1 + x2 ) ≤1 , то ∆1 > 0 и ∆2 > 0 .

Следовательно, согласно критерию Сильвестра Q(x) – положи-

тельно определенная матрица.

Задача 2. Проверить положительную определенность матрицы Гессе функции:

f (x) = 4x12 − x22 − 2x1x2 + 6x1 − x2 − 2 .

Решение. Легко установить, что матрица Гессе данной функции будет:

Q =		8	− 2		.
		− 2	− 2

По критерию Сильвестра матрица Гессе не является положи-

тельно определенной, так как её главный определитель

∆ = −16 − 4 = −20 < 0 .

Проверим полученный результат, используя собственные числа матрицы Гессе:

		8	− λ		− 2			= λ2 − 6λ − 20 = 0 ,
		− 2		− 2 − λ
			6 +			> 0 , λ2 = 6 −
λ1 =				2	116				2	116	< 0 .

Как видим, не все собственные числа положительны. Следовательно, матрица Гессе – не положительно определенная. Метод Ньютона расходится.

3.3.4. Метод переменной метрики (метод Девидона)

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

218

случаях вычисление матрицы Гессе может быть связано с большими трудностями (например, она может быть получена только численными методами). Кроме того, известно, что число умножений при обращении матрицы размера [n ×n], приблизительно, пропор-

ционально n3, что представляет собой большую величину уже при относительно небольших значениях n.

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

Основная итерационная формула метода имеет вид:

xk +1 = xk − λk ηk f (xk ),

(3.14)

где ηk – матрица, аппроксимирующая обратную матрицу Гессе, а λk – шаг, который выбирается путем минимизации функции в на-

правлении (− ηk f (xk )).

Кстати, формула (3.14) является общей для градиентных методов и метода Ньютона. В методе наискорейшего спуска роль ηk

играет единичная матрица, а в методе Ньютона – обратная матрица Гессе.

В рассматриваемом методе матрица ηk вычисляется по рекур-

рентной формуле:

ηk +1 = ηk + Ak

− Bk ,

(3.15)

где

∆

(∆

)

т

η ∆

(∆

)тη

т

x

x

g

k

g

k

A =

k

k

,

B =

k

k

,

k

(∆

)т

∆

k

(∆

)тη ∆

x

g

k

g

k

g

k

k

k

∆

k =

k +1 −

k ,

x

x

x

∆

k = f (

k +1) − f (

k ).

g

x

x

В качестве начального значения η0

для

организации рекур-

рентного процесса (3.15) принимается η0

= E , где E – единичная

матрица. Таким образом, начальное направление поиска экстремума совпадает с направлением поиска в методе наискорейшего спус-

219

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

f (x) = a + xrтb + 12 x тQx

с положительно определенной матрицей Гессе Q . Оказывается, что роль матрицы Ak состоит в том, чтобы обеспечить сходимость

матрицы ηk +1 к обратной матрице Гессе Q−1 , а роль матрицы Bk в том, чтобы обеспечить положительную определенность матрицы ηk +1 и, в пределе, исключить влияние произвольного задания мат-

рицы η0 . Действительно,

η0 = E,

η1 = E + A0 − B0 ,

η2 = η1 + A1 − B1 = E + ( A0 + A1 ) − (B0 + B1 ),

M

kk

ηk +1 = E + ∑ Ai − ∑Bi .

i=0 i=0

В случае квадратичной функции, сумма матриц Ai ( i =1, k ), при

k = n −1 , равна обратной матрице Q−1 , а сумма матриц Bi равна

матрице E . Это легко доказать, если в начале показать с помощью математической индукции, что вырабатываемые методом Девидона векторы направлений ∆x0 , ∆x1 , …, ∆xn−1 образуют множество Q -

сопряженных векторов.

Данный метод можно рассматривать как один из вариантов метода сопряженных направлений. Он решает задачу минимизации квадратичной функции за конечное число шагов, не превосходящее n .

n−1

Покажем, что ∑ Ai = Q−1 .

i=0

Из соотношений (3.15), учитывая, что f (x) = Qx + b , получим:

∆gk = f (xk +1 ) − f (xk ) = Q(xk +1 − xk ) = Q∆xk .

220