Метод нелдера мида python

Тестирование метода Хука-Дживса на функции Розенброка показало его чувствительность к выбору начальной точки поиска экстремума функции.
В этой работе, опять же на функции Розенброка, реализуется и тестируется метод Нелдера-Мида [1].
Приводится Python-реализации поиска минимума функции методом Нелдера-Мида.
При желании можно познакомиться и с Фортран-, и C#-реализациями метода Нелдера-Мида.
Замечание. Для поиска максимума нужно умножить на -1 результат, получаемый при вычислении значения оптимизируемой функции.
В качестве тестовой берется функция Розенброка:

Глобальный минимум функции равен 0.0 и находится в точке (x1, x2) = (1.0, 1.0).
При работе с функцией Розенброка для проверки используемого метода в качестве начальной часто берут точку
(-5, 10)
или точку
(-2.048, 2.048).
Приводимая программа может работать с произвольной оптимизируемой функцией.
Код, реализующий метод Нелдера-Мида, предваряет Python-программа построения графика функции Розенброка.

Вывод графика функции Розенброка

Выводится изометрическая проекция 3d-поверхности. Используется библиотека Matplotlib [2].
График показан на рис. 1.

Рис. 1. График функции Розенброка

Для его вывода употреблен следующий код:

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np

# Формирование сетки
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1.5, 3, 0.1)
X, Y = np.meshgrid(X, Y)
# Функция Розенброка
Z = (1.0 – X)**2 + 100.0 * (Y – X * X)**2
#
fig = plt.figure()
# Будем выводить 3d-проекцию графика функции
ax = fig.gca(projection = ‘3d’)
# Вывод поверхности
surf = ax.plot_surface(X, Y, Z, cmap = cm.Spectral, linew > # Метки осей координат
ax.set_xlabel(‘X’)
ax.set_ylabel(‘Y’)
ax.set_zlabel(‘Z’)
# Настройка оси X
for label in ax.xaxis.get_ticklabels():
label.set_color(‘black’)
label.set_rotation(-45)
label.set_fontsize(9)
# Настройка оси Y
for label in ax.yaxis.get_ticklabels():
label.set_fontsize(9)
# Настройка оси Z
for label in ax.zaxis.get_ticklabels():
label.set_fontsize(9)
# Изометрия
ax.view_init(elev = 30, azim = 45)
# Шкала цветов
fig.colorbar(surf, shrink = 0.5, aspect = 5)
# Отображение результата (рис. 1)
plt.show()

Python-реализация метода Нелдера-Мида

Используется библиотека SciPy [3]. Реализуются два варианта вызова библиотечной процедуры оптимизации:

1. Задается начальная точка поиска минимума функции.
2. Задается начальный симплекс поиска минимума функции.

# Вариант 1
import numpy as np
import math # Для sqrt()
import scipy.optimize as opt
# Функция Розенброка
def Rosenbrock(X):
return (1.0 – X[0])**2 + 100.0_8 * (X[1] – X[0] * X[0] )**2
#
n = 2
x0 = np.zeros(2, dtype = float) # Вектор с двумя элементами типа float
# Начальная точка поиска минимума функции
x0[0] = -5.0
x0[1] = 10.0
xtol = 1.0e-5 # Точность поиска экстремума
# Находим минимум функции
res = opt.minimize(Rosenbrock, x0, method = ‘Nelder-Mead’, options = <‘xtol’: xtol, ‘disp’: True>)
print(res)

Результат:
final_simplex: (array([[1.00000132, 1.0000028 ],
[1.00000014, 0.99999997],
[0.99999681, 0.99999355]]), array([4.38559817e-12, 9.00569749e-12, 1.05977059e-11]))
fun: 4.385598172677925e-12
message: ‘Optimization terminated successfully.’
nfev: 269 (число оценок функции)
nit: 143 (число итераций)
status: 0
success: True
x: array([1.00000132, 1.0000028 ])

# Вариант 2
import numpy as np
import math # Для sqrt()
import scipy.optimize as opt
# Функция Розенброка
def Rosenbrock(X):
return (1.0 – X[0])**2 + 100.0_8 * (X[1] – X[0] * X[0] )**2
#
# Процедура формирования начального симплекса
def makeInitialSimplex(X, L, n, initialSimplex):
qn = math.sqrt(1.0 + n) – 1.0
q2 = L / math.sqrt(2.0) * n
r1 = q2 * (qn + n)
r2 = q2 * qn
initialSimplex[0, :] = X
for j in range(n):
initialSimplex[j + 1, :] = X + r2
for i in range(n):
initialSimplex[i + 1, i] += (r1 – r2)
#
n = 2
x0 = np.zeros(2, dtype = float) # Вектор с двумя элементами типа float
# Начальная точка поиска минимума функции
x0[0] = -5.0
x0[1] = 10.0
xtol = 1.0e-5 # Точность поиска экстремума
# Начальная симплекс поиска минимума функции
initialSimplex = np.zeros((n + 1, n), dtype = float)
L = 0.4 # Длина ребра начального симплекса
# Формируем начальный симплекс
makeInitialSimplex(x0, L, n, initialSimplex)
# Находим минимум функции
res = opt.minimize(Rosenbrock, x0, method = ‘Nelder-Mead’, options = <‘xtol’: xtol, ‘disp’: True, ‘initial_simplex’: initialSimplex>)
print(res)

Результат:
final_simplex: (array([[1.00000168, 1.00000305],
[0.99999702, 0.99999376],
[0.99999679, 0.99999416]]), array([1.18711392e-11, 1.70297440e-11, 4.41663303e-11]))
fun: 1.1871139195622139e-11
message: ‘Optimization terminated successfully.’
nfev: 255 (число оценок функции)
nit: 138 (число итераций)
status: 0
success: True
x: array([1.00000168, 1.00000305])

Материал из MachineLearning.

Не путать с «симплекс-методом» из линейного программирования — методом оптимизации линейной системы с ограничениями.

Итоговое количество итераций 79. Вычисленный минимум функционала равен .

Рекомендации программисту

Метод Нелдера-Мида прост в реализации, в качестве параметров алгоритма как правило используются следующие:

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

Ниже приведен пример функции на языке C++, реализующей алгоритм.

Заключение

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

Последовательные симплексы в методе Нелдера-Мида для функции Розенброка (вверху) и функции Химмельблау (внизу)
Не путать с «симплекс-методом» из линейного программирования — методом оптимизации линейной системы с ограничениями.

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

Суть метода заключается в последовательном перемещении и деформировании симплекса вокруг точки экстремума.

Метод находит локальный экстремум и может «застрять» в одном из них. Если всё же требуется найти глобальный экстремум, можно пробовать выбирать другой начальный симплекс. Более развитый подход к исключению локальных экстремумов предлагается в алгоритмах, основанных на методе Монте-Карло, а также в эволюционных алгоритмах.

Алгоритм [ править | править код ]

Пусть требуется найти безусловный минимум функции n переменных f ( x ( 1 ) , x ( 2 ) , … , x ( n ) ) <displaystyle fleft(x^<(1)>,x^<(2)>,ldots ,x^<(n)>
ight)> . Предполагается, что серьёзных ограничений на область определения функции нет, то есть функция определена во всех встречающихся точках.

Параметрами метода являются:

• коэффициент отражения 0>"> α > 0 <displaystyle alpha >0>0"/> , обычно выбирается равным 1 <displaystyle 1>.
• коэффициент сжатия 0>"> β > 0 <displaystyle eta >0>0"/> , обычно выбирается равным 0 , 5 <displaystyle 0<,>5>.
• коэффициент растяжения 0>"> γ > 0 <displaystyle gamma >0>0"/> , обычно выбирается равным 2 <displaystyle 2>.

1. «Подготовка». Вначале выбирается n + 1 <displaystyle n+1>точка x i = ( x i ( 1 ) , x i ( 2 ) , … , x i ( n ) ) , i = 1.. n + 1 <displaystyle x_=left(x_^<(1)>,x_^<(2)>,ldots ,x_^<(n)>
   ight),i=1..n+1>, образующие симплекс n-мерного пространства. В этих точках вычисляются значения функции: f 1 = f ( x 1 ) , f 2 = f ( x 2 ) , … , f n + 1 = f ( x n + 1 ) <displaystyle f_<1>=f(x_<1>),f_<2>=f(x_<2>),ldots ,f_=f(x_)>.
2. «Сортировка». Из вершин симплекса выбираем три точки: x h <displaystyle x_>с наибольшим (из выбранных) значением функции f h <displaystyle f_>, x g <displaystyle x_>со следующим по величине значением f g <displaystyle f_>и x l <displaystyle x_>с наименьшим значением функции f l <displaystyle f_>. Целью дальнейших манипуляций будет уменьшение по крайней мере f h <displaystyle f_>.
3. Найдём центр тяжести всех точек, за исключением x h <displaystyle x_>: x c = 1 n ∑ i ≠ h x i <displaystyle x_=<frac <1>>sum limits _x_>. Вычислять f c = f ( x c ) <displaystyle f_=f(x_)>не обязательно.
4. «Отражение». Отразим точку x h <displaystyle x_>относительно x c <displaystyle x_>с коэффициентом α <displaystyle alpha >(при α = 1 <displaystyle alpha =1>это будет центральная симметрия, в общем случае — гомотетия), получим точку x r <displaystyle x_>и вычислим в ней функцию: f r = f ( x r ) <displaystyle f_=f(x_)>. Координаты новой точки вычисляются по формуле: x r = ( 1 + α ) x c − α x h <displaystyle x_=(1+alpha )x_-alpha x_>.
5. Далее смотрим, насколько нам удалось уменьшить функцию, ищем место f r <displaystyle f_>в ряду f h , f g , f l <displaystyle f_,f_,f_>. Если f r f l <displaystyle f_, то направление выбрано удачное и можно попробовать увеличить шаг. Производим «растяжение». Новая точка x e = ( 1 − γ ) x c + γ x r <displaystyle x_=(1-gamma )x_+gamma x_>и значение функции f e = f ( x e ) <displaystyle f_=f(x_)>. Если f e f r <displaystyle f_, то можно расширить симплекс до этой точки: присваиваем точке x h <displaystyle x_>значение x e <displaystyle x_>и заканчиваем итерацию (на шаг 9). Если f r f e <displaystyle f_, то переместились слишком далеко: присваиваем точке x h <displaystyle x_>значение x r <displaystyle x_>и заканчиваем итерацию (на шаг 9). Если f l f r f g <displaystyle f_, то выбор точки неплохой (новая лучше двух прежних). Присваиваем точке x h <displaystyle x_>значение x r <displaystyle x_>и переходим на шаг 9. Если f g f r f h <displaystyle f_, то меняем местами значения x r <displaystyle x_>и x h <displaystyle x_>. Также нужно поменять местами значения f r <displaystyle f_>и f h <displaystyle f_>. После этого идём на шаг 6. Если f h f r <displaystyle f_, то просто идём на следующий шаг 6. В результате (возможно, после переобозначения) f l f g f h f r <displaystyle f_.
6. «Сжатие». Строим точку x s = β x h + ( 1 − β ) x c <displaystyle x_=eta x_+(1-eta )x_>и вычисляем в ней значение f s = f ( x s ) <displaystyle f_=f(x_)>.
7. Если f s f h <displaystyle f_, то присваиваем точке x h <displaystyle x_>значение x s <displaystyle x_>и идём на шаг 9.
8. Если f_>"> f s > f h <displaystyle f_>f_>f_"/> , то первоначальные точки оказались самыми удачными. Делаем «глобальное сжатие» симплекса — гомотетию к точке с наименьшим значением x l <displaystyle x_>: x i ← x l + ( x i − x l ) / 2 <displaystyle x_gets x_+(x_-x_)/2>, i ≠ l <displaystyle i
   eq l>.
9. Последний шаг — проверка сходимости. Может выполняться по-разному, например, оценкой дисперсии набора точек. Суть проверки заключается в том, чтобы проверить взаимную близость полученных вершин симплекса, что предполагает и близость их к искомому минимуму. Если требуемая точность ещё не достигнута, можно продолжить итерации с шага 2.