В предыдущих подразделах рассматривались методы Коши и Ньютона. Отмечалось, что метод Коши эффективен при поиске на значительных расстояниях от точки минимума х* и плохо «работает» в окрестности этой точки, тогда как метод Ньютона не отличается высокой надежностью при поиске х* из удаленной точки, однако оказывается весьма эффективным в тех случаях, когда x (k) находится вблизи точки минимума. В этом и последующих подраз­делах рассматриваются методы, которые обладают положительны­ми свойствами методов Коши и Ньютона и основаны на вычислении значений только первых производных. Таким образом, эти методы, с одной стороны, отличаются высокой надежностью при поиске х* из удаленной точки х* и, с другой стороны, быстро сходятся в окрестности точки минимума.

Методы, позволяющие получать решения задач с квадратичными целевыми функциями приблизительно за N шагов при условии использования недесятичных дробей, будем называть квадратично сходящимися. Среди таких методов можно выделить класс алгоритмов, в основе которых лежит построение сопряженных направлений. Выше было сформулировано условие сопряженности для системы направлений s (k) , k = 1, 2, 3,…, r ≤ N, и симметрической матрицы С порядка N N. Была также установлена связь между построением указанных направлений и преобразованием произвольной квадратичной функции к виду суммы полных

1) Задачи такого типа возникают, например, в регрессионном анализе - Прим. перев.

квадратов; сделан вывод о том, что последовательный поиск вдоль каждого из N направлений, обладающих свойством С -сопряженности, позволяет найти точку минимума квадратичной функции N переменных. Рассмотрена процедура определения системы сопряженных направлений с использованием только значений целевой функции. Ниже для получения сопряженных направлений применяются квадратичная аппроксимация f(x) и значения компонент градиента. Кроме того, потребуем, чтобы рассматриваемые методы обеспечивали убывание целевой функции при переходе от итерации к итерации.

Пусть в пространстве управляемых переменных заданы две произвольные несовпадающие точки x (0) и x (1) . Градиент квадратичной функции равен

f(x) = q(x) = C x + b = g(x) (3.60)

Обозначение g(x) введено здесь для удобства записи формулы. Таким образом,

g (x (0)) = C x (0) + b ,

g (x (1)) = C x (1) + b .

Запишем изменение градиента при переходе от точки х (0) к точке х (1) :

g(x) = g (x (1)) – g (x (0)) = C (x (1) - x (0)), (3.61)

g(x) = C x

Равенство (3.61) выражает свойство квадратичных функций, которое будет использовано ниже.

В 1952 г. Эстенс и Штифель предложили эффективный итерационный алгоритм для решения систем линейных уравнений, который по существу представлял собой метод сопряженных градиентов. Они рассматривали левые части линейных уравнений как компоненты градиента квадратичной функции и решали задачу минимизации этой функции. Позже Флетчер и Ривс обосновали квадратичную сходимость метода и обобщили его для случая неквадратичных функций. Фрид и Метцлер продемонстрировали (допустив, однако, некоторые неточности) возможности использования метода для решения линейных систем с разреженной матрицей коэффициентов. (Определение разреженной матрицы см. в приложении А.) Они подчеркнули простоту реализации метода по сравнению с другими, более общими алгоритмами, что является особенно важной характеристикой с позиций нашего изложения.

Рассмотрение метода будем проводить в предположении, что "целевая функция является квадратичной:

f(x) = q(x) = a + b T x + ½ x T C x ,

аитерации проводятся по формуле (3.42), т.е.

x = x + α s(x ) .

Направления поиска на каждой итерации определяются с помощью следующих формул:

s (k) = – g (k) + (3.62)

s (0) = –g (0) , (3.63)

где g (k) = f (x ). Так как после определения системы направлений проводится последовательный поиск вдоль каждого из направлений, полезно напомнить, что в качестве критерия окончания одномерного поиска обычно используется условие

f (x ) T s (k) = 0 (3.64)

Значения , i = 1, 2, 3,...,k - 1 ,выбираются таким образом, чтобы направление s (k) было С -сопряжено со всеми построенными ранее направлениями поиска. Рассмотрим первое направление

s (1) = –g (1) + γ (0) s (0) = –g (1) –γ (0) g (0)

и наложим условие его сопряженности с s (0)

s (1)T C s (0) = 0,

откуда T C s (0) = 0.

На начальной итерации

s (0) = ;

следовательно,

T C = 0

Используя свойство квадратичных функций (3.61), получаем

T g = 0, (3.65)

γ (0) = – ( g T g (1))/( g T g (0)). (3.66)

Из уравнения (3.65) следует, что

g (1)T g (1) + γ (0) g (0)T g (1) – g (1) T g (0) – γ (0) g (0)T g (0) = 0.

При соответствующем выборе α (0) и с учетом формулы (3.64) имеем

g (1) T g (0) = 0.

Таким образом,

s (2) = –g (2) + γ (0) s (0) + γ (1) s (1) .

и выберем γ (0) γ (1) таким образом, чтобы выполнялись условия

s (2) T C s (0) = 0 и s (2) C s (1) = 0,

т. е. условия С -сопряженности направления s (2) с направлениями s (0) и s (1) . С помощью формул (3.61) и (3.64) можно показать (это предоставляется читателю в качестве упражнения), что здесь γ (0) = 0, а в общем случае γ ( i ) = 0, i = 0, 1, 2,...,k-2, при любом значении k. Отсюда следует, что общая формула для направлений поиска может быть записана в виде, предложенном Флетчером и Ривсом:

s (k ) = –g (k ) + s (3.68)

Если f(x) - квадратичная функция, для нахождения точки минимума требуется определить N -1 таких направлений и провести N поисков вдоль прямой (при отсутствии ошибок округления). Если же функция f(х) не является квадратичной, количество направлений и соответствующих поисков возрастает.

Некоторые исследователи на основе опыта проведения вычислительных экспериментов предлагают после реализации каждой серии из N или N + 1 шагов возвращаться к начальной итерации алгоритма, положив s(x) = -g(x). Это предложение остается предметом для изучения, поскольку при минимизации функций общего вида в ряде случаев влечет за собой замедление сходимости. С другой стороны, циклические переходы к начальной итерации повы­шают надежность алгоритма, так как вероятность построения линейно зависимых направлений уменьшается. Если полученное на­правление оказывается линейной комбинацией одного или нескольких полученных ранее направлений, то метод может не привести к получению решения, поскольку поиск в соответствующем подпространстве R N уже проводился. Однако следует отметить, что на практике такие ситуации встречаются достаточно редко. Метод оказывается весьма эффективным при решении практических задач, характеризуется простотой однопараметрической вычислительной схемы и небольшим объемом памяти ЭВМ, необходимым для проведения поиска. Относительно невысокий уровень требований к объему памяти ЭВМ делает метод Флетчера - Ривса (ФР) и его модификации особенно полезным при решении задач большой размерности.

Пример 3.9. Метод Флетчера - Ривса

Найти точку минимума функции

f(x) = 4x + 3x – 4x x + x

если x = T .

Шаг 1. f(x) = T ,

s (0) = – f (x (0)) = – T .

Шаг 2.Поиск вдоль прямой:

x = x – α f (x (0)) → α = ⅛,

x = T – ⅛ T = [– ⅛, 0] T

Шаг 3.k = 1.

s (1) = – T – [¼, 1] T = [– ¼, ½] T .

Шаг 4.Поиск вдоль прямой:

x = x + α s (1) → α = ¼,

x = [– ⅛, 0] T – ¼ [¼, ½] T = [– , ] T ,

f (x (2)) = T .

Таким образом, x = х*. Решение получено в результате проведения двух одномерных поисков, поскольку целевая функция квадратичная, а ошибки округления отсутствуют.

Миль и Кентрелл обобщили подход Флетчера и Ривса, предложив формулу

x = x + α {– f (x (k )) + γ s (x )} (3.69)

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

Напомним, что рассмотрение метода Флетчера - Ривса прово­дилось в предположении квадратичности целевой функции и отсутствия ошибок округления в процессе поиска вдоль прямой. Од­нако при реализации ряда методов сопряженные направления определяются без учета одного из указанных предположений (или даже обоих предположений). Среди таких методов наибольшего внима­ния, по-видимому, заслуживает метод, разработанный Ползком и Рибьером в 1969 г. Метод основан на точной процедуре проведения поиска вдоль прямой и на более общем предположении об аппроксимации целевой функции. При этом

γ = , (3.70)

где, как и прежде,

g (x )= g (x ) – g (x ). (3.71)

Если α - параметр, значение которого определяется в результате поиска вдоль прямой, и γ - параметр, значение которого вычисляется по формуле (3.70), то метод сопряженных градиентов Полака - Рибьера реализуется с помощью следующих соотношений:

x = x + α s (x ),

s (x ) = – f (x ) + γ s (x ). (3.72)

Легко видеть, что единственное различие между методами Полака - Рибьера и Флетчера - Ривса заключается в способах выбора параметра γ.

Известны и другие подобные методы, которые также основаны на проведении точных вычислений при одномерном поиске и на более общей (чем квадратичная) аппроксимации целевой функции (см., например, ). Краудер и Вульф в 1972 г., а затем Пауэлл доказали, что методы сопряженных градиентов обладают линейной скоростью сходимости при отсутствии периодического возврата к начальной итерации. Указанные возвраты осуществляются в соответствии со специальной процедурой, которая прерывает процесс построения векторов направлений поиска, а далее вычисления продолжаются по аналогии с построением s (x (0)). Выше отмечалось, что по ряду причин наличие процедуры возврата к начальной итерации повышает устойчивость работы алгоритма, так как позволяет избежать построения линейно зависимых векторов направлений поиска. Пауэлл доказал, что метод Полака - Рибьера также характеризуется линейной скоростью сходимости при отсутствии возвратов к начальной итерации, однако имеет несомненное преимущество перед методом Флетчера - Ривса при решении задач с целевыми функциями общего вида и обладает менее высокой чувствительностью к ошибкам округления при проведении одномерных поисков.

Вопросы разработки эффективных процедур и методов, обеспечивающих возврат к начальной итерации и при этом обладающих малой чувствительностью к ошибкам округления, остаются предметом активных исследований. Бил предложил метод сопряженных градиентов, аналогичный стандартной процедуре Флетчера - Ривса, но вместе с тем не использующий направление вдоль градиента при возвращении к начальной итерации. Он показал, как на основе анализа направления, полученного непосредственно перед возвращением к начальной итерации, можно уменьшить объем необходимых вычислений при решении задач, требующих нескольких возвратов. Пауэлл исследовал стратегию Била, а также другие стратегии возврата к начальной итерации и предложил использовать процедуру возврата либо после проведения каждой серии из N шагов, либо при выполнении неравенства

| g (x ) T g (x ) | ≥ 0.2 ||g (x )|| . (3.73)

Он продемонстрировал, что стратегию Била, дополненную услови­ем (3.73), можно успешно применять как вместе с формулой Флет­чера - Ривса, так и с формулой Полака - Рибьера, и провел ряд вычислительных экспериментов, результаты которых подтверж­дают превосходство метода Полака - Рибьера (с возвратом). Шэнно исследовал влияние ошибок округления при проведе­нии поисков вдоль прямой и различных стратегий возврата на эффективность методов сопряженных градиентов. Он показал, что стратегия Била (с использованием соответствующей двухпараметрической формулы), дополненная предложенным Пауэллом усло­вием возврата, приводит к существенному уменьшению требуемой точности одномерного поиска и, следовательно, к значительному повышению эффективности полной вычислительной схемы метода сопряженных градиентов. Шэнно также представил численные ре­зультаты, которые указывают на преимущество метода Полака - Рибьера с использованием процедур возврата и округления при поисках вдоль прямой. В работе продемонстрирована ведущая роль методов сопряженных градиентов при решении задач нелиней­ного программирования большой размерности.

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

Эти методы подобны методам, рассмотренным в подразд. 3.3.5, поскольку также основаны на свойствах квадратичных функций. В соответствии с изложенными выше методами поиск решения осуществляется по системе сопряженных направлений, тогда как квазиньютоновские методы обладают положительными чертами метода Ньютона, однако используют только первые производные. Во всех методах указанного класса построение векторов направлений поиска осуществляется с помощью формулы (3.42), в которой s (x (k))записывается в виде

s (x ) = –A f (x ), (3.74)

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

x = C g . (3.75)

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

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

A = A + A (3.76)

где A - корректирующая матрица. Матрица A будет использоваться в формулах (3.74) и (3.42). Задача заключается в том, чтобы построить матрицу A таким образом, чтобы последовательность А (0) , А (1) , А (2) ,...,A (k +1) давала приближение к Н -1 = f (x *) -1 ; при этом для получения решения х* требуется один дополнительный поиск вдоль прямой, если f(x) - квадратичная функция. Как неоднократно подчеркивалось выше, имеются определенные основания полагать, что метод, обеспечивающий нахождение оптимумов квадратичных функций, может привести к успеху при решении задач с нелинейными целевыми функциями общего вида.

Вернемся к важному свойству квадратичных функций (3.75) и предположим, что матрица С -1 аппроксимируется по формуле

С -1 = βA , (3.77)

где р - скалярная величина. Наиболее предпочтительным является приближение, удовлетворяющее (3.75), т. е.

x = A g . (3.78)

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

x = x – x , (3.79)

g = g (x ) – g (x ). (3.80)

С другой стороны, можно потребовать, чтобы новое приближение удовлетворяло формуле (3.75):

x = βA g . (3.81)

Подставляя выражение (3.76) в (3.81), получаем

A g = x – A g . (3.82)

С помощью непосредственной подстановки можно убедиться, что матрица

A = – (3.83)

является решением этого уравнения. Здесь у и z - произвольные векторы, т. е. формула (3.83) определяет некоторое семейство решений. Если положить

y = x и z =A g , (3.84)

то получим формулу, реализующую известный и широко приме­няемый метод Дэвидона - Флетчера - Пауэлла (ДФП) :

A = A + – . (3.85)

Можно показать, что эта рекуррентная формула сохраняет свойства симметрии и положительной определенности матриц. Поэтому если А (0) симметрическая положительно определенная матрица,то матрицы А (1) , А (2) , ... также оказываются симметрическими и положительно определенными при отсутствии ошибок округления; обычно удобно выбирать А (0) = I .

Первая вариация f(x) равна

f(x) = f (x ) x . (3.86)

Используя формулы (3.42) и (3.74), получаем

f(x) = f (x ) α A f (x ), (3.87)

f(x) = – α f (x ) A f (x ), (3.88)

и неравенство f (x (k +1) ) < f (x k ) выполняется при любых значениях α > 0, если A - положительно определенная матрица. Таким образом, алгоритм обеспечивает убывание целевой функции при переходе от итерации к итерации. Метод Дэвидона - Флетчера - Пауэлла в течение ряда лет продолжает оставаться наиболее широко используемым градиентным методом. Он отличается устойчивостью и успешно применяется при решении самых различных задач, возникающих на практике. Основным недостатком методов такого типа является необходимость хранить в памяти ЭВМ матрицу А порядка N N.

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

1. Понятие о методах сопряженных направлений.

Рассмотрим задачу минимизации квадратичной функции

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

Будем говорить, что ненулевые векторы являются взаимно сопряженными (относительно матрицы А), если для всех

Под методом сопряженных направлений для минимизации квадратичной функции (10.33) будем понимать метод

в котором направления взаимно сопряжены, а шаги

получаются как решение задач одномерной минимизации:

Теорема 10.4. Метод сопряженных направлений позволяет найти точку минимума квадратичной функции (10 33) не более чем за шагов.

Методы сопряженных направлений отличаются один от другого способом построения сопряженных направлений. Наиболее известным среди них является метод сопряженных градиентов

2. Метод сопряженных градиентов.

В этом методе направления строят правилу

Так как то первый шаг этого метода совпадает с шагом метода наискорейшего спуска. Можно показать (мы этого делать не будем), что направления (10.34) действительно являются

сопряженными относительно матрицы А. Более того, градиенты оказываются взаимно ортогональными.

Пример 10.5. Применим метод сопряженных градиентов для минимизации квадратичной функции - из примера 10.1. Запишем виде где

Возьмем начальное приближение

1-й шаг метода совпадает с первым шагом метода наискорейшего спуска. Поэтому (см. пример 10.1)

2-й шаг. Вычислим

Так как то и решение оказалось найденным за два шага.

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

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

или к виду

Преимущество формул (10 36), (10.37) в том, что они не содержат явным образом матрицу А.

Минимизацию функции методом сопряженных градиентов производят в соответствии с формулами

Коэффициенты вычисляют по одной из формул (10.36), (10.37).

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

Решение задач одномерной минимизации (10.40) приходится осуществлять численно. Отметим также то, что часто в методе сопряженных градиентов при коэффициент не вычисляют по формулам (10.36), (10.37), а полагают равным нулю. При этом очередной шаг производят фактически методом наискорейшего спуска. Такое "обновление" метода позволяет уменьшить влияние вычислительной погрешности.

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

f (x )

		f (xk ) = f(x0 ) + ∑ α i Api .
		i= 1
	обе части	этого равенства скалярно на p k		учитывая
исчерпывающего спуска по направлению p k :(f (x k ), p k ) = 0			и A −ортогональность
векторов, получаем
		(f (x 0 ),p k )+ α k (Ap k ,p k )= 0.
		A положительно определена,	квадратичная
(Ap k , p k ) > 0 и для величины шагаα k получаем выражение (5.17).
	Последовательный исчерпывающий спуск			A –ортогональным

направлениям (5.16) приводит к точке минимума квадратичной формы не более чем за n шагов.

□ Доказать самостоятельно. Предположить, что существуют u k ≠ α k , и

получить, что они совпадают. ■

Вопрос о нахождении базиса из A –ортогональных векторов в пространствеE n

решается неоднозначно. В качестве такого базиса можно, например, взять ортогональный базис из собственных векторов матрицы A . Однако их поиск приn > 2 представляет собой самостоятельную и довольно сложную задачу.

и без предварительного построения векторов p 1 , ..., p n , последовательно находя их в процессе минимизации, как это было сделано выше в примере с минимизацией функции двух переменных. И в этом случае для квадратичной функции с положительно определенной матрицейA для нахождения минимума достаточно конечное число шагов. Если не является квадратичной функцией или

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

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

5.6. Метод сопряженных градиентов

При использовании методов градиентного и наискорейшего спуска в итерационной процедуре

антиградиента: p k = − f (x k ). Однако такой выбор направления убывания не всегда бывает удачным. В частности, для плохо обусловленных задач минимизации направление антиградиента в точкеx k может значительно отличаться от направления к точке минимумаx . В результате траектория приближения к точке минимума имеет зигзагообразный характер. Воспользуемся другим подходом, идея которого была изложена при построении метода сопряженных направлений. Будем определять направления спускаp k не только через вектор антиградиента− f (x k ) ,

в котором величина шага α k находится из условия исчерпывающего спуска по

направлению p k . Далее,		после вычисления очередной точки x k + 1 ,					k = 0, 1, ..., новое
направление поиска p k + 1		находится по формуле, отличной от антиградиента:
	pk + 1 = − f(xk + 1 ) + β k pk ,			k = 0, 1, ...,
где коэффициенты		выбираются так, чтобы при минимизации квадратичной
функции f (x ) с	положительно определенной					матрицей	A получалась
последовательность	A −ортогональных			векторов		p 0 ,p 1 , ....	Из условия
(Ap k + 1 ,p k )= 0имеем:
		β k=	(A f (x k + 1 ),p k )
			(Ap k ,p k )

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

для квадратичной функции шаг исчерпывающего спуска по направлению p k равен

α k = −	(f (x k ),p k )
	(Ap k ,p k )
Утверждение . Итерационный процесс			(5.19)−(5.22) минимизации

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

f (x )

A дает точки	x 0 , ...,x k	и векторы p 0 , ..., p k такие, что если				f (x i )≠ 0при
0 ≤i	то векторы		p 0 , ...,		A −ортогональны,		градиенты
f (x 0 ), ...,f (x i )	взаимно ортогональны.
Так как направления				являются A −ортогональными,

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

С учетом взаимной ортогональности градиентов f (x i ) и условий

исчерпывающего спуска по направлениям p k можно упростить выражения (5.21) и

(5.22) для α k

и β k . В результате получим,

что итерационный процесс метода

сопряженных градиентов описывается соотношениями

x k+ 1

X k +α k

p k ,k = 0, 1, ...;

x0 En ,

p0 = − f(x0 ) ,

f (x k + α k p k )= minf (x k

+ αp k ),

k = 0, 1, ...,

α> 0

p k+ 1

= − f (x k + 1 ) +β k

p k ,k = 0, 1, ...,

β k=

f (xk + 1 )

k = 1, 2, ...

f (xk )

Следует отметить, что выражение для коэффициента β k не содержит в явном виде матрицуA квадратичной формы. Поэтому метод сопряженных градиентов может применяться для минимизации неквадратичных функций.

Итерационный процесс (5.23)−(5.26) может не приводить к точке минимума неквадратичной функции за конечное число итераций. Более того, точное

определение α k из условия (5.22) возможно лишь в редких случаях, а вектораp k

на образуют, вообще говоря, A −ортогональную систему относительно какой-либо матрицыA . Поэтому реализация каждой итерации метода будет сопровождаться неизбежными погрешностями. Эти погрешности, накапливаясь, могут привести к

тому, что векторы p k перестанут указывать направление убывания функции и

сходимость метода может нарушаться. Поэтому в методе сопряженных градиентов применяется практический прием − через каждые N шагов производят обновление метода, полагаяβ m N = 0, m = 1, 2, ... . Номераm N называют моментами

обновления метода, или рестарта . Часто полагаютN = n − размерности пространстваE n . ЕслиN = 1 , то получается частный случай метода сопряженных градиентов − метод наискорейшего спуска.

Вблизи точки минимума дважды дифференцируемая функция с положительно определенной матрицей Гессе H (x ) , как правило, достаточно хорошо

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

Пример 5.7. Методом сопряженных градиентов найти точку минимума

функции f (x ) = 4 x 2 + 3 x 2 − 4 x x

из начальной точки x 0 = (0, 0) T .

□ Итерация 1.

Шаг 1. Положим ε = 0,01,

= (0, 0)T ,

и найдем f (x 0 ) = (1, 0) T . Перейдем к

Шаг 2. Положим k = 0,

= − f (x 0 ) = (− 1, 0) T . Перейдем к шагу 3.

f (x0

+ α p 0 )→ min.Получим

α 0 = 1/ 8 . – Здесь применили формулуα 0 = −

(f (x 0 ),p 0 )

(Ax 0 + b ,p 0 )

Перейдем

(Ap 0 ,p 0 )

(Ap 0 ,p 0 )

Шаг 4. Найдем

x 1= x 0

+ α 0 p 0 = (− 1/ 8,

и f (x 1 ) = (0, 1/ 2) T . Точность не

достигнута, прейдем к шагу 5.

Шаг 5. Условие k + 1 = n не выполняется (нет рестарта), перейдем к шагу 6.

Шаг 6. Найдем коэффициент β 0 = 1/ 4 и новое направление спуска

p 1 = − f (x 1 ) + β 0 p 0 = (− 1/ 4, − 1/ 2) T . Перейдем к следующей итерации.

Поскольку x 1 , f (x 1 ) иp 1

= − f (x 1 ) +β 0

уже вычислены на итерации 1, то

итерацию 2 начинаем с шага 3.

Итерация 2.

Шаг 3. Решим задачу одномерной минимизации

f (x 1 + α p 1 ) → min . Получим

α = 1/ 4 . Перейдем к шагу 4.

Шаг 4. Найдем x 2

X 1 +α 1

p 1 = (− 3 /16,− 1/ 8)T и f (x 2 )= (0, 0)T − задача решена

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

Описание алгоритма

Шаг 0 . Выбирается точка начального приближения , параметр длины шага , точность решения и вычисляется начальное направление поиска .

Шаг k . На k -м шаге находится минимум (максимум) целевой функции на прямой, проведенной из точки по направлению . Найденная точка минимума (максимума) определяет очередное k -е приближение , после чего определяется направление поиска

Формула (5.4) может быть переписана в эквивалентном виде

.

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

Метод Ньютона

Метод предназначен для решения задачи (5.1) и принадлежит классу методов второго порядка. В основе метода лежит разложение Тейлора целевой функции и то, что в точке экстремума градиент функции равен нулю, то есть .

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

. (5.5)

Формулу (5.5) перепишем в матричной форме, учитывая при этом, что :

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

Предположим, что матрица Гессе невырождена. Тогда она имеет обратную матрицу . Умножая обе части уравнения (5.6) на слева, получим , откуда

. (5.7)

Формула (5.7) определяет алгоритм метода Ньютона: пересчет приближений на k

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

,

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

Метод Ньютона-Рафсона

Метод является методом первого порядка и предназначен для решения систем n нелинейных уравнений c n неизвестными:

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

Пусть точка есть решение системы (5.9), а точка расположена вблизи . Разлагая функцию в точке по формуле Тейлора, имеем

откуда (по условию ) вытекает

, (5.11)

где матрица Якоби вектор-функции . Предположим, что матрица Якоби невырождена. Тогда она имеет обратную матрицу . Умножая обе части уравнения (5.11) на слева, получим , откуда

. (5.12)

Формула (5.12) определяет алгоритм метода Ньютона-Рафсона: пересчет приближений на k -й итерации выполняется в соответствии с формулой

В случае одной переменной, когда система (5.9) вырождается в единственное уравнение , формула (5.13) принимает вид

, (5.14)

где значение производной функции в точке .

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

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

Замечание 5.2. Методы Ньютона и Ньютона-Рафсона требуют большого объема вычислений (надо на каждом шаге вычислять и обращать матрицы Гессе и Якоби).

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

ЛИТЕРАТУРА

1. Афанасьев М.Ю. , Суворов Б.П. Исследование операций в экономике: Учебное пособие. – М.: Экономический факультет МГУ, ТЕИС, 2003 – 312 с.

2. Базара М, Шетти К. Нелинейное программирование. Теория и алгоритмы: Пер. с англ. – М.: Мир, 1982 – 583 с.

3. Берман Г .Н . Сборник задач по курсу математического анализа: Учебное пособие для вузов. – СПб: «Специальная Литература», 1998. – 446 с.

4. Вагнер Г. Основы исследования операций: В 3-х томах. Пер. с англ. – М.: Мир, 1972. – 336 с.

5. Вентцель Е. С. Исследование операций. Задачи, принципы, методология – М.: Наука, 1988. – 208 с.

6. Демидович Б.П. Сборник задач и упражнений по математическому анализу. – М.: Наука, 1977. – 528 с.

7. Дегтярев Ю.И. Исследование операций. – М.: Высш. шк., 1986. – 320 с.

8. Нуреев Р.М. Сборник задач по микроэкономике. – М.: НОРМА, 2006. – 432 с.

9. Солодовников А. С., Бабайцев В.А., Браилов А.В. Математика в экономике: Учебник: В 2-х ч. – М.:Финансы и статистика, 1999. – 224 с.

10. Таха Х. Введение в исследование операций, 6-е изд.: Пер. с англ. – М.: Издательский дом «Вильямс», 2001. – 912 с.

11. Химмельблау Д. Прикладное нелинейное программирование: Пер. с англ. – М.: Мир, 1975 – 534 с.

12. Шикин Е. В., Шикина Г.Е. Исследование операций: Учебник – М.: ТК Велби, Изд-во Проспект, 2006. – 280 с.

13. Исследование операций в экономике : Учебн. пособие для вузов/ Н.Ш.Кремер, Б.А.Путко, И.М.Тришин, М.Н.Фридман; Под ред. проф. Н.Ш.Кремера. – М.: Банки и биржи, ЮНИТИ, 1997. – 407 с.

14. Матрицы и векторы : Учебн. пособие/ Рюмкин В.И. – Томск: ТГУ, 1999. – 40 с.

15. Системы линейных уравнений : Учебн. пособие / Рюмкин В.И. – Томск: ТГУ, 2000. – 45 с.

ВВЕДЕНИЕ……………………………………...................................
1. ОСНОВЫ МАТЕМАТИЧЕСКОГО ПРОГРАММИРОВАНИЯ………………...
1.1. Постановка задачи математического программирования...............................
1.2. Разновидности ЗМП…………….…………..........................................
1.3. Базовые понятия математического программирования................................
1.4. Производная по направлению. Градиент………….........................................
1.5. Касательные гиперплоскости и нормали…………..........................................
1.6. Разложение Тейлора……………………………...............................................
1.7. ЗНЛП и условия существования ее решения...................................................
1.8. Задачи ……………..……...................................................................................
2. РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ БЕЗ ОГРАНИЧЕНИЙ................................................................................................................
2.1. Необходимые условия решения ЗНЛП без ограничений...............................
2.2. Достаточные условия решения ЗНЛП без ограничений.................................
2.3. Классический метод решения ЗНЛП без ограничений...................................
2.4. Задачи……………..............................................................................................
3. РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ ПРИ ОГРАНИЧЕНИЯХ-РАВЕНСТВАХ.................................................................................
3.1. Метод множителей Лагранжа…………………………...................................
3.1.1. Назначение и обоснование метода множителей Лагранжа……………
3.1.2. Схема реализации метода множителей Лагранжа……………………...
3.1.3. Интерпретация множителей Лагранжа…………………………………
3.2. Метод подстановки…………………………….................................................
3.3. Задачи…………………………..........................................................................
4. РЕШЕНИЕ ЗАДАЧИ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ ПРИ ОГРАНИЧЕНИЯХ-НЕРАВЕНСТВАХ………………………………………………..
4.1. Обобщенный метод множителей Лагранжа…………………………………
4.2. Условия Куна-Таккера…………………………..............................................
4.2.1. Необходимость условий Куна-Таккера…………………………………
4.2.2. Достаточность условий Куна-Таккера…………………………………..
4.2.3. Метод Куна-Таккера………………………...............................................
4.3. Задачи…………………………..........................................................................
5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ …………………………...……………………………………
5.1. Понятие алгоритма…………………………....................................................
5.2. Классификация численных методов…………………………………………
5.3. Алгоритмы численных методов……………………………………………...
5.3.1. Метод наискорейшего спуска (подъема)…………………………………
5.3.2. Метод сопряженных градиентов………………………….........................
5.3.3. Метод Ньютона………………………….....................................................
5.3.4. Метод Ньютона-Рафсона………………………………………………...
ЛИТЕРАТУРА………………………………..............................................................

Определения линейной и нелинейной функций см. в разделе 1.2

Введение

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

Постановка задачи оптимизации

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

Если , то задача оптимизации называется безусловной (unconstrained ). Если , то задача оптимизации называется условной (constrained ).

Метод сопряжённых градиентов для квадратичного функционала

Изложение метода

Рассмотрим следующую задачу оптимизации:

Здесь - симметричная положительно определённая матрица размера . Такая задача оптимизации называется квадратичной. Заметим, что . Условие экстремума функции эквивалентно системе Функция достигает своей нижней грани в единственной точке , определяемой уравнением . Таким образом, данная задача оптимизации сводится к решению системы линейных уравнений
Идея метода сопряжённых градиентов состоит в следующем:
Пусть - базис в . Тогда для любой точки вектор раскладывается по базису Таким образом, представимо в виде

Каждое следующее приближение вычисляется по формуле:

Определение. Два вектора и называются сопряжёнными относительно симметричной матрицы B, если

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

Базисные вектора вычисляются по формулам:

Коэффициенты выбираются так, чтобы векторы и были сопряжёнными относительно А.

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

Анализ метода

Для метода сопряжённых градиентов справедлива следующая теорема:
Теорема Пусть , где - симметричная положительно определённая матрица размера . Тогда метод сопряжённых градиентов сходится не более чем за шагов и справедливы следующие соотношения:

Сходимость метода

Если все вычисления точные, и исходные данные точны то метод сходится к решению системы не более чем за итераций, где - размерность системы. Более тонкий анализ показывает, что число итераций не превышает , где - число различных собственных значений матрицы A. Для оценки скорости сходимости верна следующая (довольно грубая) оценка:

, где

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

.

Вычислительная сложность

На каждой итерации метода выполняется операций. Такое количество операций требуется для вычисления произведения - это самая трудоёмкая процедура на каждой итерации. Отальные вычисления требуют O(n) операций. Суммарная вычислительная сложность метода не превышает - так как число итераций не больше n.

Численный пример

Применим метод сопряжённых градиентов для решения системы , где

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

Заключение

Метод сопряжённых градиентов - один из наиболее эффективных методов решения СЛАУ с положительно определённой матрицей. Метод гарантирует сходимость за конечное число шагов, а нужная точность может быть достигнута значительно раньше. Основная проблема заключается в том, что из-за накопления погрешностей может нарушаться ортогональность базисных веторов , что ухудшает сходимость

Метод сопряжённых градиентов в общем случае

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

.

Непрерывно дифференцируемая в функция. Чтобы модифицировать метод сопряжённых градиентов для решения этой задачи необходимо получить для формулы, в которые не входит матрица А:

Можно вычислять по одной из трёх формул:

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

Анализ метода

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

Сходимость метода

Для метода Флетчера - Ривса существует теорема о сходимости, накладывающая не слишком жёсткие условия на минимизируемую функцию :
Теорема .
Пусть и выполняются следующие условия:

множества M: .
Тогда

Для метода Полака-Райбера доказана сходимость в предположении, что - строго выпуклая функция. В общем случае доказать сходимость метода Полака - Райбера невозможно. Напоротив, верна следующая теорема:
Теорема.
Предположим, что в методе Полака-Райбера значения на каждом шаге вычисляются точно. Тогда существует функция , и начальное приближение , такие что %200,%20%5Cforall%20k%20=%200,%201,%202,%20...%20%5Cquad%20%7C%7Cf(x_k)%7C%7C%20>%20%5Cdelta" alt="\exists \delta > 0, \forall k = 0, 1, 2, ... \quad ||f(x_k)|| > \delta">.

Тем не менее, на практике метод Полака-Райбера работает лучше.
Наиболее распространённые критерии останова на практике: Норма градиента становится меньше некоторого порога
Значение функции в течении m последовательных итераций почти не изменилось

Вычислительная сложность

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

Числовой пример

Будем искать методом сопряжённых градиентов минимум функции . Минимум этой функции равен 1 и достигается в точке (5, 4). Сравним на примере этой функции методы Полака-Райбера и Флетчера-Ривса. Итерации в обоих методах прекращаются, когда на текущем шаге квадрат нормы градиента становится меньше . Для выбора используется метод золотого сечения

	Метод Флетчера - Ривса			Метод Полака - Райбера
	Число итераций	Найденное решение	Значение функции	Число итераций	Найденное решение	Значение функции
0.01	18	(5.01382198,3.9697932)	1.00110367	15	(5.03942877,4.00003512)	1.00155463
0.001	20	(5.01056482,3.99018026)	1.00020805	18	(4.9915894,3.99999044)	1.00007074
0.0001	24	(4.9979991,4.00186173)	1.00000747	20	(5.00336181,4.0000018)	1.0000113
0.00001	25	(4.99898277,4.00094645)	1.00000193	22	(4.99846808,3.99999918)	1.00000235
0.00001	29	(4.99974658,4.0002358)	1.00000012	26	(4.99955034,3.99999976)	1.0000002

Реализовано два варианта метода сопряжённых градиентов: для минимизации квадратичного функционала, и для минимизации произвольной функции. В первом случае метод реализуется функцией
vector FindSolution(matrix A, vector b)
Здесь A и b - матрица и вектор, фигурирющие в определении квадратичной задачи оптимизации.
Для минимизации произвольного функционала можно использовать одну из двух функций:
vector FletcherRievesMethod(int spaceSize, Function F, vector (*GradF) (vector))
vector PolakRibiereMethod(int spaceSize, Function F, vector (*GradF) (vector))
Параметры для обеих функций совпадают и имеют следующий смысл:
spaceSize - размерность пространства(число переменных, от которых зависит минимизируемый функционал)
F - указатель на минимизируемую функцию. Функция должна иметь вид double <имя функции>(vector)
GradF - указатель на функцию, вычисляющую градиент минимизируемого функционала
Оба метода используют вспомогательную функцию для решения задачи одномерной оптимизации. В программе реализована одномерная оптимизация методом золотого сечения.

Заключение

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

См. также

Список литературы

• Васильев Ф. П. Методы оптимизации - Издательство «Факториал Пресс», 2002
• Nocedal J., Wright S.J. Numerical Optimization ,Springer, 1999