Алгоритм метода Нелдера–Мида

1. Проецируют точку x [ h, k ] через центр тяжести:

x [ n + 3, k ] = x [ n + 2, k ] + a (x [ n + 2, k ] – x [ h, k ]),

где a – некоторая константа (обычно а = 1).

2. Растягивают вектор x [ n + 3, k ] – x [ n + 2, k ]:

x [ n + 4, k ] = x [ n + 2, k ] + с (x [ n + 3, k ] – x [ n + 2, k ]),

где с > 1 – коэффициент растяжения (обычно 2,8 ¸ 3).

Если f (x [ n + 4, k ]) < f (x [ l, k ]), то x [ h, k ] заменяют на x [ n + 4, k ] и продолжают с п. 1 при k = k + 1.

Иначе x [ h, k ] заменяют на x [ n + 3, k ] и продолжают с п. 1 при k = k + 1.

3. Если f (x [ n + 3, k ]) > f (x [ i, k ]) для всех i ¹ h, то сжимают вектор x [ h, k ] – x [ n + 2, k ]:

x [ n + 5, k ] = x [ n + 2, k ] + b (x [ h, k ] – x [ n + 2, k ]),

где b > 0 – коэффициент сжатия (обычно 0,4 ¸ 0,6).

Затем точку x [ h, k ] заменяют на x [ n + 5, k ] и продолжают с п. 1 при k = k + 1.

4. Если f (x [ n + 3, k ]) > f (x [ h, k ]), то все векторы x [ i, k ] – x [ l, k ], I = 1, 2, …, n + 1 уменьшают в два раза:

x [ i, k ] = x [ l, k ] + 0,5(x [ i, k ] – x [ l, k ]).

Затем продолжают с п. 1 при k = k + 1.

Выход из процесса происходит при предельном сжатии многогранника:

где e = (e₁, …, e_n) – заданный вектор.

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

Метод вращающихся координат (метод Розенброка)

Суть метода состоит во вращении системы координат в соответствии с изменением скорости убывания целевой функции.

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

Идея метода состоит в следующем (рис. 9.6).

Из начальной точки x [0] осуществляют спуск в точку x [1] по направлениям, параллельным координатным осям.

На следующей итерации одна из осей должна проходить в направлении y ₁ = х [1] – х [0], а другая – в направлении, перпендикулярном к y ₁.

Спуск вдоль этих осей приводит в точку x [2], что дает возможность построить новый вектор х [2] – х [1] и на его базе новую систему направлений поиска.

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

Рис. 9.6. Геометрическая интерпретация метода Розенброка

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

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

Метод параллельных касательных (метод Пауэлла)

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

Суть метода: выбирается начальная точка х [0] и выполняется одномерный поиск вдоль произвольного направления, приводящий в точку х [1].

Затем выбирается точка х [2], не лежащая на прямой х [0] – х [1], и проводится одномерный поиск вдоль прямой, параллельной х [0] – х [1].

Полученная в результате точка х [3] вместе с точкой х [1] определяет направление х [1] – х [3] одномерного поиска, дающее точку минимума х *.

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

 

Рис. 9.7. Геометрическая интерпретация метода Пауэлла

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

9.3 Методы безусловной оптимизации первого порядка

Градиентом дифференцируемой функции f (x) в точке х [0] называется n-мерный вектор f ' (х [0]), компоненты которого являются частными производными функции f (x), вычисленными в точке х [0], т. е.

^.

Этот вектор перпендикулярен к плоскости, проведенной через точку х [0], и касательной к поверхности уровня функции f (x), проходящей через точку х [0].

В каждой точке такой поверхности функция f (x) принимает одинаковое значение.

Приравнивая функцию различным постоянным величинам С₀, С₁,..., получим серию поверхностей, характеризующих ее топологию (рис. 9.8).

Рис. 9.8. Градиент функции

Вектор-градиент направлен в сторону наискорейшего возрастания функции в данной точке.

Вектор, противоположный градиенту (– f ' (x [0])), называется антиградиентом и направлен в сторону наискорейшего убывания функции.

В точке минимума градиент функции равен нулю.

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

Выбирая в качестве направления спуска p [ k ] антиградиент в точке x [ k ], получаем итерационный процесс вида

x [ k + 1] = x [ k ] – a _k , a _k > 0; k = 0, 1, 2, ….

Вектор-градиент направлен в сторону наискорейшего возрастания функции в данной точке.

В координатной форме этот процесс записывается следующим образом:

Критерий останова – либо условие малости приращения аргумента,

|| x [ k + 1] – x [ k ] || < e,

либо выполнение условия малости градиента

|| || < g.

Здесь e, g – заданные малые величины.

Возможен и комбинированный критерий.

В методе с постоянным шагом для всех итераций выбирается некоторая постоянная величина шага. Достаточно малый шаг a_k обеспечит убывание функции.

Цена – (возможно) очень большое количество итераций для достижения точки минимума.

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

Методы с постоянным шагом применяются на практике редко.

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

Метод наискорейшего спуска

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

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

С математической точки зрения на каждой итерации необходимо решать задачу одномерной минимизации по α функции φ (α) = f (x [ k ] – α f ' (x [ k ])).

В рассматриваемом методе направление движения из точки х [ k ] касается линии уровня в точке
x [ k + 1] (рис. 9.9).

Рис. 9.9. Геометрическая интерпретация метода наискорейшего спуска

Траектория спуска – зигзаг, соседние звенья зигзага ортогональны друг другу. Действительно, шаг a _k выбирается путем минимизации по a функции

j (a) = f (x [ k ] – a f ' (x [ k ])).

Необходимое условие минимума функции dj (a)/ d a = 0.

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

dj (a)/ d a = – f ' (x [ k +1 ]) f ' (x [ k ]) = 0.

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

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

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

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

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

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

Затем применяют другие, более эффективные методы.

Рис. 9.10. Овражная функция

 

9.4 Методы оптимизации второго порядка

Методы безусловной оптимизации второго порядка используют вторые частные производные минимизируемой функции f (х). Эти методы называют иногда квази - ньютоновскими.

Идея методов

Необходимым условием экстремума функции многих переменных f (x)в точке х* является равенство нулю ее градиента в этой точке:

Разложение в окрестности точки x [ k ] в ряд Тейлора с точностью до членов первого порядка позволяет переписать предыдущее уравнение в виде

Здесь – матрица вторых производных (матрица Гессе) минимизируемой функции:

(9.2)

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

х [ k +1 ] = х [ k ] – H ^–1 (x [ k ]) f ' (x [ k ]),

где H ^–1 (x [ k ]) – обратная матрица для матрицы Гессе, а – H ^–1 (x [ k ]) f ' (x [ k ]) = р [ k ] – направление спуска.

Полученный метод минимизации называют методом Ньютона.

Очевидно, что в данном методе величина шага вдоль направления р [ k ] полагается равной единице.

Последовательность точек { x [ k ]}, получаемая в результате применения итерационного процесса, при определенных предположениях сходится к некоторой стационарной точке x * функции f (x).

Если матрица Гессе H (x *)положительно определена, точка x *, будет точкой строгого локального минимума функции f (x).

Последовательность { x [ k ]} сходится к точке x *, только если матрица Гессе целевой функции положительно определена на каждой итерации.

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

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

При минимизации «овражных» функций скорость сходимости метода Ньютона более высока по сравнению с градиентными методами.

Существенный недостаток метода Ньютона – зависимость сходимости для невыпуклых функций от начального приближения x [0].

Если x [0] далека от точки минимума, метод может расходиться, т. е. каждая следующая точка будет более удаленной от точки минимума, чем предыдущая.

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

p [ k ] = – H ^–1 (x [ k ]) f ' (x [ k ]),

но и величины шага a вдоль этого направления.

Соответствующий алгоритм называют методом Ньютона с регулировкой шага. Употребляется также термин «метод с переменной метрикой ».

Итерационный процесс в таком случае определяется выражением:

x [ k + 1] = x [ k ] – a _k H ^–1(x [ k ]) f ' (x [ k ]).

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

Вследствие накопления ошибок в процессе счета матрица Гессе на некоторой итерации может оказаться отрицательно определенной или ее нельзя будет обратить. В таких случаях в подпрограммах оптимизации полагается H ^–1(x [ k ]) = Е, где Е – единичная матрица. Очевидно, что итерация при этом осуществляется по методу наискорейшего спуска.

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

Причина – необходимость вычисления и обращения матрицы вторых производных целевой функции.

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

Он обладает сверхлинейной, или квадратичной, скоростью сходимости в зависимости от требований, которым удовлетворяет минимизируемая функция f (x).

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

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

В начале процесса минимизации, когда точка x [0] находится далеко от точки экстремума x *, можно применять какой-либо вариант градиентных методов. Далее, при уменьшении скорости сходимости градиентного метода можно перейти к методу Ньютона.

Существуют и широко применяются две основные модификации алгоритма метода с переменной метрикой:

• алгоритм Дэвидона–Флетчера–Пауэлла (Davidon–Fletcher–Powell, DFP алгоритм);

• алгоритм Бройдена–Флетчера–Гольдфарба–Шанно (Broyden–Fletcher–Goldfarb–Shanno, BFGS алгоритм).

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

Тема 10

Метод Монте–Карло

10.0 Метод Монте–Карло. Некоторые задачи

Задача Бюффона

Название «метод Монте–Карло» для методов, систематически использующих случайную величину, появилось во второй половине 40-х годов ХХ века, когда фон Нейман и Улам использовали случайные числа для моделирования поведения нейтронов.

Сама идея восходит, по меньшей мере, к французскому естествоиспытателю Бюффону (XVIII век).

Формулировка задачи: на листе бумаги имеется набор параллельных прямых на расстоянии друг от друга. С некоторой высоты на лист роняют иголку длиной L < d. Упав, иголка случайным образом располагается на бумаге, пересекая, либо не пересекая линии. Из количества пересечений и общего числа всех результатов надо оценить число p (рис. 10.1).

Рис. 10.1. К задаче Бюффона

Пусть d = 1, а длина иголки L < d.

Когда иголку бросают, центр ее может упасть от прямой на любом расстоянии от 0 до 1/2. Обозначим эту переменную через х.

Угол наклона иголки к прямым обозначим через j (рис. 10.1). Переменные х и j можно считать случайными и независимыми.

Условие пересечения: .

На рис. 10.2 заштрихована область, ограниченная кривой x = L∙ cos(j/2), в которую должна попадать точка с координатами (j, х), чтобы иголка пересекала одну из линий.

Рис. 10.2. К задаче Бюффона. Условие пересечения одной из линий

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

Можно рассматривать рассуждения Бюффона как определение значения p методом Монте–Карло. Если бы относительно p было известно только, что оно лежит между 1 и 10, это был бы прекрасный способ обнаружить, что p примерно равно 3. Сделав много бросаний, мы могли бы получить p примерно 3,1.

Ошибка уменьшается , где N – число испытаний.

Метод Монте–Карло называют также методом статистических испытаний.

Площадь круга

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

Из 15 нанесенных точек внутрь попали 11. Отсюда площадь четверти круга 11/15 ≈ 0,73. Т.к. эта площадь равна p/4, получаем оценку для p ≈ 2,9.

Взяв из той же таблицы 150 точек, получим оценку для p ≈ 3,4 (!!).

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

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

В рассматриваемом примере можно принять, что ошибка при определении некоторого числа A (число испытаний = N) равна

 

Рис. 10.3. Вычисление площади круга

Т. о. при N =150 мы получили Δp ≈ 0,28, так что наше p должно быть записано в виде p ≈ 3,4 ± 0,28.

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

Это число равно N = A ²/(Δ A)².

Для получения p с точностью до одной сотой нужно около 100 000 испытаний.

10.1 Проблема получения случайных чисел

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

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

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

Применялись два основных подхода.

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

Достоинства: высокая «степень случайности» (при соблюдении некоторых дополнительных условий).

Недостатки: порой слишком медленно; невоспроизводимость последовательностей (т. е. невозможно повторить расчет с тем же набором случайных чисел).

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

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

Недостаток: последовательность всегда одна и та же (от библиотечной функции типа RANDOM).

10.2 Общая схема метода

Допустим, что нам требуется вычислить какую-то неизвестную величину m.

Попытаемся предложить такую случайную величину x, чтобы М x = m. Пусть при этом дисперсия D x = b ².

Рассмотрим N случайных независимых величин x₁, x₂, …, x _N, распределения которых совпадают с распределением x.

Если N достаточно велико, то согласно предельной центральной теореме распределение суммы r _N = x₁+ x₂ +…+ x _N будет приблизительно нормальным с параметрами a = Nm, s² = Nb².

Из «правила трех сигм» следует, что

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

Перепишем это соотношение в несколько ином виде:

(10.1)

Соотношение (10.1) – ключевое для метода Монте–Карло. Оно дает и метод расчета m, и оценку погрешности.

В самом деле, найдем N значений случайной величины x.

Важно: безразлично, находить ли один раз по одному значению каждой из величин x₁, x₂, …, x _N или найти N значений одной величины x, т. к. все эти случайные величины одинаковы (имеют одно и то же распределение).

Из (10.1) видно, что среднее арифметическое этих значений приближенно равно m.

С большой вероятностью ошибка этого приближения не превосходит и стремится к нулю с ростом N.

10.3 Метод Монте–Карло. Пример