Симплексный метод

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

Симплекс обладает замечательным свойством: если взять только одну точку вне симплекса в качестве новой вершины, то, соединив ее с вершинами прилежащей грани, получим новый симплекс. На рис. 8.23 слева из исходного симплекса с вершинами 1,2,3 построено 2 других симплекса, а справа новый симплекс построен на грани 2,3,4 исходного и новой вершины 5.

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

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

Последовательное отражение наихудших вершин перемещает симплекс в направлении уменьшения значения функции, что в конечном итоге приводит к отысканию минимума с заданной точностью. Характер движения симплекса представлен на рис. 8.25. Отражение вершин показано пунктирными линиями. Сначала отражается вершина 1, затем в полученном симплексе отражается вершина 3, симплекс 2,4,5 заменяется на 4,5,6 и так далее. Очевидно, что в процессе движения, особенно вблизи экстремума, новая вершина может оказаться не лучше отраженной. Для предотвращения зацикливания либо отражается другая вершина, либо уменьшается размер нового симплекса.

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

где i – номер вершины, j – координата вершины (индекс переменной), n – число переменных (размерность пространства), а приращения d ₁ и d ₂ зависят от n и выбранного масштабного множителя m:

При m =1 все ребра имеют единичную длину.

Вычисление координат новой вершины, получающейся в результате отражения одной из вершин последнего симплекса, производится следующим образом. Пусть отражается вершина k (например, вершина 1 на рис. 8.24). Тогда сначала определяется центр противолежащей грани C

а затем координаты новой вершины s:

x ^(s) = x ^(k) + l ×(x ^C – x ^(k))

или

x ^(s) = x ^С + q ×(x ^C – x ^(k)), (8.42)

где l >0 и q – параметры отражения.

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

Значения коэффициентов отражения были подобраны экспериментально: a = 1; b = 0,5; g = 2. Предложены и другие наборы значений, однако они дают лучшие результаты только в частных случаях.

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

Сначала производится нормальное отражение. Получаемая в результате точка D считается временной (рис. 8.26). Далее идут проверки:

· если f^L < f^D < f^B, нормальное отражение приемлемо и временная точка фиксируется;

· если f^B £ f^D < f^A, то выполняется сжатие с q =b и фиксируется точка D _b;

· если f^D ³ f^А, симплекс сжимается с q =–b и фиксируется точка D_{– b} ;

· если f^D < f^L, то производится растяжение (q =g), дающее точку D _g, которая фиксируется, если она не хуже временной; иначе она не учитывается и фиксируется временная точка.

В условии останова поиска может использоваться показатель размера симплекса, например, максимальная длина ребра, и/или разность значений целевой функции f^А – f^L.

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

Интересен вариант метода, в котором на всех итерациях симплекс остается регулярным. В нем используется коэффициент уменьшения длины ребра b ³ 0,85. Симплекс сжимается, стягиваясь к лучшей вершине, как показано на рис. 8.27, и после этого происходит отражение. Такое значение b обусловлено тем, что на вершины уменьшенного симплекса переносятся значения функции соответствующих вершин до сжатия. Тем самым исключаются дополнительные вычисления функции. Возможно также увеличение симплекса за счет удлинения всех ребер с коэффициентом 1/ b при сохранении положения наихудшей вершины.

8.8.4. Градиентные методы

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

Скорость изменения функции f (X) в произвольном направлении l определяется производной

. (8.43)

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

.

Поверхность уровня целевой функции, описываемая уравнением f (X) = const, имеет размерность n – 1.

Зададим координаты следующим образом. Проведем через рассматриваемую точку n – 1 взаимно ортогональные касательные к поверхности уровня, приняв их за оси координат. В качестве последней оси возьмем нормаль к поверхности уровня (рис. 8.29).

В такой системе координат производные функции по всем x_j равны нулю. Поэтому из (8.43) имеем

.

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

Обозначается градиент как grad f (X) или Ñ f (X). Он полностью определяется своими проекциями – производными функции по переменным:

В задачах на минимум используется противоположное направление – антиградиент.

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

,

.

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

Если переменные имеют разные единицы измерения, можно перейти к относительным переменным y_i, используя минимально и максимально возможные значения переменных x_i:

.

Очевидно, что значения y_i лежат в диапазоне [0, 1].

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

В градиентном методе поиск минимума производится согласно рекуррентным формулам:

, (8.44)

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

Формула (8.44) в записи по отдельным переменным принимает вид:

. (8.45)

Алгоритм градиентного метода.

1) Задать начальную точку, начальное значение шага и точность по величине градиента e, вычислить f (X⁰) и положить k = 0.

2) В текущей точке X ^k вычислить градиент (производные и его длину.

3) Проверить: если закончить поиск.

4) Определить новую точку X ^{k +1} согласно (8.45) и вычислить в ней значение функции.

5) Проверить: если f (X ^{k+ 1})£ f (X ^k), увеличить k на единицу и перейти на шаг 2; если иначе, уменьшить h^k в два раза и перейти на шаг 4.▲

При гладкой производной траектория поиска в идеале (при непрерывном вычислении градиента) тоже будет гладкой и ортогональной линиям уровня. При конечной величине шага траектория становится кусочно-линейной. Качественный характер поиска показан на рис. 8.30, а на рис. 8.31 приведена реализация алгоритма с двух начальных точек применительно к функции f (X)=(x ₁- x ₂)²+(x ₂-2)⁴. Обе траектории приводят практически в одну точку.

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

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

(8.46)

В результате решения задачи (8.46) определяется оптимальный шаг h ^*, по которому находится следующая точка:

(8.47)

Очевидно, что при таком определении новой точки, значение функции в ней будет лучше, чем в X ^k.

Алгоритм наискорейшего спуска.

1) Задать начальную точку и точность по величине градиента e, положить k = 0.

2) В текущей точке X ^k вычислить градиент (производные и его длину.

3) Проверить: если закончить поиск.

4) Провести одномерный поиск (8.46).

5) Определить новую точку X ^{k +1} согласно (8.47), увеличить k на единицу и перейти на шаг 2. ▲

На рис. 8.32 показаны две траектории движения к минимуму функции f (X)=(x ₁- x ₂)²+(x ₂-2)⁴, полученные алгоритмом. Минимум на направлении антиградиента достигается в точке касания с линией уровня, а градиент в этой точке ортоганален ей. Поэтому каждое последующее направление ортогонально непосредственно предшест­вую­щему. Из рис. 8.32 видно, что с приближением к экстремуму частота вычисления градиента увеличивается, и вблизи минимума метод наискорейшего спуска вырождается в градиентный (см. рис. 8.31 и 8.32).

Градиентные методы плохо работают в условиях оврага: при попадании на дно оврага резко падает скорость движения и возможна остановка поиска до достижения окрестности минимума (рис. 8.33, f =(2- x ₁)²+3(x ₁²- x ₂)²).

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

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

Квадратичная аппроксимация дифферен­цируемой функции f в точке X ^k записывается в виде

, (8.48)

где Н – матрица вторых производных функции f (матрица Гессе). В стационарной точке градиент функции равен нулю, что применительно к (8.48) дает

. (8.49)

Полагая матрицу H неособенной (существует обратная к ней матрица), из (8.49) получаем рекуррентный способ генерирования последовательности точек

. (8.50)

Исходя из вывода формулы (8.50) ясно, что для квадратичной функции цели X ^{k +1} является точкой минимума. Следовательно, минимум такой функции из любой начальной точки достигается за один шаг (8.50). Для нелинейной унимодальной функции, отличной от квадратичной, точки последовательности (8.50) асимтотически приближаются к минимуму. Условие окончания поиска такое же, как и в градиентных методах: . (В случае линейной аппроксимации матрица Н становится единичной и поиск вырождается в градиентный).

Необходимым условием сходимости метода является существование обратной матрицы во всех генерируемых точках. Доказательство сходимости метода получено при условии, что начальная точка достаточно близка к точке минимума. При этом метод имеет высокую скорость сходимости. Если поиск начинается с удаленной точки и функция существенно отличается от квадратичной, метод может не сходиться или давать слабую сходимость. Причина такого поведения в том, что значение функции в точке X ^{k+ 1} не обязательно меньше, чем в точке X ^k. Так на рис. 8.34 траектория поиска имеет зигзагообразный характер из-за того, что после 2-й итерации значение функции ухудшилось.

Для устранения отмеченных недостатков предложены модификации метода.

Чтобы обеспечить регуляризацию матрицы Н (невырожденность), она деформируется добавлением к ней единичной матрицы с неотрицательным скалярным множителем e: H^¢= Н+ e ×Е. Значение e зависит от X и подбирается так, чтобы все собственные значения матрицы H^¢ были положительными. Тогда эта матрица положительно определена и существует обратнгая матрица H^¢-1, которая также положительно определена.

Возможное ухудшение функции в очередной точке исключается введением в формулу (8.50) изменяемого параметра h. Модифицированная формула принимает вид

. (8.51)

При этом возможен вариант с дискретным шагом h, как в градиентном методе, либо с определением оптимального h^* с помощью одномерной минимизации (наискорейший спуск):

Пример поиска минимума функции
f =(1- x ₁)²+5(x ₁²- x ₂)² методом Ньютона с регулируемым шагом приведен на рис. 8.35. Если сравнить траектории поиска на рис. 8.34 и 8.35, то легко заметить преимущество модифицированного варианта метода.

8.8.6. Методы сопряженных направлений

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

Пусть дана матрица Н _{n ´ n}. Направления d ₁, d ₂,..., d _k (k £ n) называются сопряженными или Н-сопряженными, если они линейно независимы и

. (8.52)

Эти векторы определяют сопряженные направления. Для квадратичной функции двух переменных сопряженные направления получаются следующим образом. Возьмем произвольное направление d ₁ и на нем найдем минимум, двигаясь из точки X¹. Повторим поиск минимума на d ₁ из точки X²¹ X¹ (рис. 8.36). Направление d ₂, определяемое прямой, проходящей через найденные минимумы, является сопряженным с направлением d ₁. При этом направление d ₂ проходит через искомый минимум функции f. Следовательно, при любой начальной точке минимум квадратичной функции двух переменных достигается за два одномерных поиска вдоль сопряженных направлений.

Пример 8.6. Используя сопряженные направления, найти минимум функции (точка минимума X^*=(2,4)).

Запишем матрицу Гессе .

За первое направление возьмем . Компоненты d ₂ найдем из условия (8.52):

.

Положив а = 1, получаем b = 2 и . Возьмем начальную точку X⁰=(-1;1). Найдем минимум на направлении d ₁. Для этого подставим в функцию X = X⁰+ h d ₁, то есть x ₁= x ₁⁰+ h = -1+ h, x ₂= x ₂⁰=1. Тогда f = h ²-3 h -3 и минимум по h будет при h ^*=1,5. Следовательно, минимум на d ₁ достигается в точке X¹=(0,5;1). Приняв ее за начальную для поиска вдоль d ₂ и подставляя в функцию x ₁= 0,5+ h, x ₂=1+2 h, получаем f = 3 h ²-9 h -5,25. Находим h ^*=1,5 и соответствующую новую точку X²=(2;4). Как видим, второй одномерный поиск привел в точку искомого минимума f (рис. 8.37).▲

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

Методы, основанные на концепции сопряженных направлений, различаются способами построения таких направлений. Ряд из них относятся к прямым методам, например, метод Пауэлла и его модификация – метод Зангвилла. Другие используют первые производные, например, метод сопряженных градиентов (метод Флетчера-Ривса). Одним из самых эффективных является метод Дэвидона-Флетчера-Пауэлла. В нем генерируются сопряженные направления с использованием градиента и матрицы D, аппроксимирующей обратную матрицу H^-1. Поэтому его относят также к квазиньютоновским методам. Рассмотрение этих методов выходит за рамки настоящего пособия.

8.8.7. Методы случайного поиска

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

Случайный механизм выбора направления реализуется с помощью датчика случайных чисел b, равномерно распределенных на интервале [- b, b ]. Направление задается случайным вектором

X = (x ₁, x ₂, x ₃,..., x _n),

компоненты которого вычисляются по формуле

,

где n случайных чисел b_i генерируются датчиком. Очевидно, что такой случайный вектор имеет единичную длину и определяет только направление. При этом все направления равновероятны.

Приведем несколько простых алгоритмов случайного поиска.