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

21.09.2019

Вычислительные эксперименты для оценки эффективности параллельного варианта метода верхней релаксации проводились при условиях, указанных во введении. С целью формирования симметричной положительно определенной матрицы элементы подматрицы генерировались в диапа-зоне от 0 до 1, значения элементов подматрицы получались из симметрии матриц и , а элементы на главной диагонали (подматрица ) генерировались в диапазоне от до , где – размер матрицы.

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

Таблица 7.20. Результаты экспериментов (метод верхней релаксации)

n	1 поток	Параллельный алгоритм
		2 потока		4 потока		6 потоков		8 потоков
		T	S	T	S	T	S	T	S
2500	0,73	0,47	1,57	0,30	2,48	0,25	2,93	0,22	3,35
5000	3,25	2,11	1,54	1,22	2,67	0,98	3,30	0,80	4,08
7500	7,72	5,05	1,53	3,18	2,43	2,36	3,28	1,84	4,19
10000	14,60	9,77	1,50	5,94	2,46	4,52	3,23	3,56	4,10
12500	25,54	17,63	1,45	10,44	2,45	7,35	3,48	5,79	4,41
15000	38,64	26,36	1,47	15,32	2,52	10,84	3,56	8,50	4,54

Рис. 7.50.

Эксперименты демонстрируют неплохое ускорение (порядка 4 на 8-и потоках).

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

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

Обращается в ноль. Таким образом, решение системы (7.49) можно искать как решение задачи безусловной минимизации (7.56).

Последовательный алгоритм

С целью решения задачи минимизации (7.56) организуется следующий итерационный процесс.

Подготовительный шаг () определяется формулами

Где – произвольное начальное приближение; а коэффициент вычисляется как

Основные шаги (

) определяются формулами

Здесь – невязка -го приближения, коэффициент находят из условия сопряженности

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

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

Таким образом, выполнение итераций метода потребует

(7.58)

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

Организация параллельных вычислений

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

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

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

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

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

Где – длина вектора, – число потоков, – накладные расходы на созда- ние и закрытие параллельной секции.

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

Где – число итераций метода.

Результаты вычислительных экспериментов

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

Результаты вычислительных экспериментов приведены в таблице 7.21 (время работы алгоритмов указано в секундах).

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

Несколько слов об обозначениях, используемых далее.

Скалярное произведение двух векторов записывается x T y и представляет сумму скаляров: . Заметим, что x T y = y T x. Если x и y ортогональны, то x T y = 0. В общем, выражения, которые преобразуются к матрице 1х1, такие как x T y и x T Ax, рассматриваются как скалярные величины.

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

Ax = b (1)

где x – неизвестный вектор, b – известный вектор, а A – известная, квадратная, симметричная, положительно–определенная матрица. Решение этой системы эквивалентно нахождению минимума соответствующей квадратичной формы.
Квадратичная форма – это просто скаляр, квадратичная функция некого вектора x следующего вида:

f(x) = (1/2)x T Ax-b T x+c (2)

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

x T Ax > 0 (3)

На рисунке 1 изображено как выглядят квадратичные формы соответственно для положительно-определенной матрицы (а), отрицательно-определенной матрицы (b), положительно-неопределенной матрицы (с), неопределенной матрицы (d).

Рис. 1. Квадратичные формы для положительно-определенной матрицы, отрицательно-определенной матрицы, положительно-неопределенной матрицы, неопределенной матрицы.

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

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

в начальной точке x(0) вычисляется градиент, и движение осуществляется в направлении антиградиента до тех пор, пока уменьшается целевая функция;
в точке, где функция перестает уменьшаться, опять вычисляется градиент, и спуск продолжается в новом направлении;
процесс повторяется до достижения точки минимума.

Рис. 2. Траектория движения в точку минимума методом наискорейшего спуска.

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

Рис. 3. Траектория движения в точку минимума при использовании метода сопряженных градиентов

Определение сопряженности формулируется следующим образом: два вектора x и y называют А-сопряженными (или сопряженными по отношению к матрице А) или А–ортогональными, если скалярное произведение x и Ay равно нулю, то есть:

x T Ay = 0 (4)

Сопряженность можно считать обобщением понятия ортогональности. Действительно, когда матрица А – единичная матрица, в соответствии с равенством 4, векторы x и y – ортогональны. Можно и иначе продемонстрировать взаимосвязь понятий ортогональности и сопряженности: мысленно растяните рисунок 3 таким образом, чтобы линии равного уровня из эллипсов превратились в окружности, при этом сопряженные направления станут просто ортогональными.

Остается выяснить, каким образом вычислять сопряженные направления. Один из возможных способов – использовать методы линейной алгебры, в частности, процесс ортогонализации Грамма–Шмидта. Но для этого необходимо знать матрицу А, поэтому для большинства задач (например, обучение многослойных нейросетей) этот метод не годится. К счастью, существуют другие, итеративные способы вычисления сопряженного направления, самый известный – формула Флетчера-Ривса:

(6)

Формула 5 означает, что новое сопряженное направление получается сложением антиградиента в точке поворота и предыдущего направления движения, умноженного на коэффициент, вычисленный по формуле 6. Направления, вычисленные по формуле 5, оказываются сопряженными, если минимизируемая функция задана в форме 2. То есть для квадратичных функций метод сопряженных градиентов находит минимум за n шагов (n – размерность пространства поиска). Для функций общего вида алгоритм перестает быть конечным и становится итеративным. При этом, Флетчер и Ривс предлагают возобновлять алгоритмическую процедуру через каждые n + 1 шагов.

Можно привести еще одну формулу для определения сопряженного направления, формула Полака–Райбера (Polak-Ribiere):

(7)

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

Из приведенного алгоритма следует, что на шаге 2 осуществляется одномерная минимизация функции. Для этого, в частности, можно воспользоваться методом Фибоначчи, методом золотого сечения или методом бисекций. Более быструю сходимость обеспечивает метод Ньютона–Рафсона, но для этого необходимо иметь возможность вычисления матрицы Гессе. В последнем случае, переменная, по которой осуществляется оптимизация, вычисляется на каждом шаге итерации по формуле:

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

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

Вариант метода сопряженных направлений, использующий формулу Флетчера-Ривса для расчета сопряженных направлений.

K:= 0
r:= -f"(x) // антиградиент целевой функции
d:= r // начальное направление спуска совпадает с антиградиентом
Sigma new: = r T * r // квадрат модуля антиградиента
Sigma 0: = Sigma new

// Цикл поиска (выход по счетчику или ошибке)
while i < i max and Sigma new > Eps 2 * Sigma 0
begin
j: = 0
Sigma d: = d T * d

// Цикл одномерной минимизации (спуск по направлению d)
repeat
a: =
x: = x + a
j: = j + 1
until (j >= j max) or (a 2 * Sigma d <= Eps 2)

R: = -f"(x) // антиградиент целевой функции в новой точке
Sigma old: = Sigma new
Sigma new: = r T * r
beta: = Sigma new / Sigma old
d: = r + beta * d // Вычисление сопряженного направления
k: = k + 1

If (k = n) or (r T * d <= 0) then // Рестарт алгоритма
begin
d: = r
k: = 0
end

I: = i + 1
end

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

Литература

Н.Н.Моисеев, Ю.П.Иванилов, Е.М.Столярова "Методы оптимизации", М. Наука, 1978
А.Фиакко, Г.Мак-Кормик "Нелинейное программирование", М. Мир, 1972
У.И.Зангвилл "Нелинейное программирование", М. Советское радио, 1973
Jonathan Richard Shewchuk "Second order gradients methods", School of Computer Science Carnegie Mellon University Pittsburg, 1994

В предыдущих подразделах рассматривались методы Коши и Ньютона. Отмечалось, что метод Коши эффективен при поиске на значительных расстояниях от точки минимума х* и плохо «работает» в окрестности этой точки, тогда как метод Ньютона не отличается высокой надежностью при поиске х* из удаленной точки, однако оказывается весьма эффективным в тех случаях, когда 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.

Метод предназначен для решения задачи (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

Несколько слов об обозначениях, используемых далее.

Скалярное произведение двух векторов записывается $x^Ty$ и представляет сумму скаляров: $\sum_{i=1}^n\, x_i\,y_i$. Заметим, что $x^Ty = y^Tx$. Если x и y ортогональны, то $x^Ty = 0$. В общем, выражения, которые преобразуются к матрице 1х1, такие как $x^Ty$ и $x^TA_x$, рассматриваются как скалярные величины.

$f\,(x) = (\frac{1}{2})\,x^TA_x\,-\,b^Tx\,+\,c$, (2)

$x^TA_x\,>\,0$, (3)

в начальной точке x(0) вычисляется градиент, и движение осуществляется в направлении антиградиента до тех пор, пока уменьшается целевая функция;
в точке, где функция перестает уменьшаться, опять вычисляется градиент, и спуск продолжается в новом направлении;
процесс повторяется до достижения точки минимума.

$x^TA_y\,=\,0$, (4)

$d_{(i+1)} = d_{(i+1)}\,+\,\beta_{(i+1)}\,d_i$ , (5)

$\beta_{(i+1)} = \frac{r_{(i+1)}^T}{r_{i}^T}\,\frac{r_{(i+1)}}{r_{(i)}}$, (6)

Можно привести еще одну формулу для определения сопряженного направления, формула Полака–Райбера (Polak-Ribiere):

$\beta_{(i+1)} = \frac{r_{(i+1)}^T\,(r_{(i+1)}\,-\,r_{(i)})}{r_{i}^T\,r_{(i)}}$, (7)

Метод Флетчера-Ривса сходится, если начальная точка достаточно близка к требуемому минимуму, тогда как метод Полака-Райбера может в редких случаях бесконечно циклиться. Однако последний часто сходится быстрее первого метода. К счастью, сходимость метода Полака-Райбера может быть гарантирована выбором $\beta = \max \{\beta;\,0\}$. Это эквивалентно рестарту алгорима по условию $\beta \leq 0$. Рестарт алгоритмической процедуры необходим, чтобы забыть последнее направление поиска и стартовать алгоритм заново в направлении скорейшего спуска.

Вычисляется антиградиент в произвольной точке x (0) .
$d_{(0)} = r_{(0)} = -\,f"(x_{(0)})$
Осуществляется спуск в вычисленном направлении пока функция уменьшается, иными словами, поиск a (i) , который минимизирует
$f\,(x_{(i)}\,+\,a_{(i)}\,d_{(i)})$
Переход в точку, найденную в предыдущем пункте
$x_{(i+1)} = x_{(i)}\,+\,a_{(i)}\,d_{(i)}$
Вычисление антиградиента в этой точке
$r_{(i+1)} = -\,f"(x_{(i+1)})$
Вычисления по формуле 6 или 7. Чтобы осуществить рестарт алгоритма, то есть забыть последнее направление поиска и стартовать алгоритм заново в направлении скорейшего спуска, для формулы Флетчера–Ривса присваивается 0 через каждые n + 1 шагов, для формулы Полака-Райбера – $\beta_{(i+1)} = \max \{\beta_{(i+1)},\,0\}$
Вычисление нового сопряженного направления
$d_{(i+1)} = r_{(i+1)}\,+\,\beta_{(i+1)}\,d_{(i)}$
Переход на пункт 2.

$$a = -\,\frac{{f"}^T\,(x)\,d}{d^T\,f""(x)\,d}$$

$f""(x)\,= \begin{pmatrix} \frac{\partial^2\,f}{\partial x_1\,\partial x_1}&\frac{\partial^2\,f}{\partial x_1\,\partial x_2}&\cdots&\frac{\partial^2\,f}{\partial x_1\,\partial x_n}& \\ \frac{\partial^2\,f}{\partial x_2\,\partial x_1}&\frac{\partial^2\,f}{\partial x_2\,\partial x_2}& \cdots&\frac{\partial^2\,f}{\partial x_2\,\partial x_n}& \\ \vdots&\vdots&\ddots&\vdots &\\ \frac{\partial^2\,f}{\partial x_n\,\partial x_1}& \frac{\partial^2\,f}{\partial x_n\,\partial x_2}&\cdots&\frac{\partial^2\,f}{\partial x_n\,\partial x_n} \end{pmatrix}$
Матрица Гессе

r:= -f"(x) // антиградиент целевой функции

d:= r // начальное направление спуска совпадает с антиградиентом

Sigma new: = r T * r // квадрат модуля антиградиента

Sigma 0: = Sigma new

// Цикл поиска (выход по счетчику или ошибке)
while i < i max and Sigma new > Eps 2 * Sigma 0
begin
j: = 0
Sigma d: = d T * d

// Цикл одномерной минимизации (спуск по направлению d)
repeat
a: =
x: = x + a
j: = j + 1
until (j >= j max) or (a 2 * Sigma d <= Eps 2)

If (k = n) or (r T * d <= 0) then // Рестарт алгоритма
begin
d: = r
k: = 0
end

I: = i + 1
end

Литература

Н.Н.Моисеев, Ю.П.Иванилов, Е.М.Столярова "Методы оптимизации", М. Наука, 1978
А.Фиакко, Г.Мак-Кормик "Нелинейное программирование", М. Мир, 1972
У.И.Зангвилл "Нелинейное программирование", М. Советское радио, 1973
Jonathan Richard Shewchuk "Second order gradients methods", School of Computer Science Carnegie Mellon University Pittsburg, 1994

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

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

Последовательный алгоритм

Организация параллельных вычислений

Результаты вычислительных экспериментов

Литература

Интересные статьи

Интересные статьи