SciPy, оптимизация

SciPy (произнасян като сай пай) е пакет от приложни математически процедури, основан на разширението Numpy Python. С SciPy интерактивната сесия Python се превръща в пълноценна среда за обработка на данни и прототипиране на сложни системи, подобна на MATLAB, IDL, Octave, R-Lab и SciLab. Днес искам накратко да разкажа как да прилагаме някои известни алгоритми за оптимизация в пакета scipy.optimize. По-подробна и актуална справка за приложението на функциите винаги може да се получи с командата help() или с Shift+Tab.

Въведение

За да спестя време на себе си и на читателите при търсене и четене на първоначалните източници, линковете за описания на методите ще бъдат основно към Wikipedia. В повечето случаи тази информация е достатъчна за разбиране на методите в общи линии и условията за тяхното приложение. За по-добро разбиране на същността на математическите методи се насочваме към линковете към по-авторитетни публикации, които могат да се намерят в края на всяка статия или в любимата търсачка.

И така, модулът scipy.optimize включва реализация на следните процедури:

1. Условна и безусловна минимизация на скаларни функции на няколко променливи (minim) с помощта на различни алгоритми (симплекс на Нелдер-Мийд, BFGS, сопряжени градиенти на Нютон, COBYLA и SLSQP)
2. Глобална оптимизация (например: basinhopping, diff_evolution)
3. Минимизиране на остатъците МНК (least_squares) и алгоритми за подгонка на криви с нелинейни МНК (curve_fit)
4. Минимизиране на скаларни функции на една променлива (minim_scalar) и търсене на корени (root_scalar)
5. Многомерни решатели на системи уравнения (root) с помощта на различни алгоритми (хибриден Паул, Левенберг-Маркварт или мащабни методи, като Нютон-Крилов).

В тази статия ще разгледаме само първата точка от списъка.

Безусловна минимизация на скалярната функция на няколко променливи

Функцията minim от пакета scipy.optimize предлага общ интерфейс за решаване на задачи за условна и безусловна минимизация на скаларни функции на няколко променливи. За да демонстрираме работата ѝ, ще ни е нужна подходяща функция на няколко променливи, която ще минимизираме по различни начини.

За тези цели функцията на Розенброк от N променливи е идеална, която има вида:

Въпреки че функцията на Розенброк и матриците на Якоби и Хес (първи и втори производни съответно) вече са дефинирани в пакета scipy.optimize, ще я дефинираме сами.

import numpy as np

def rosen(x):
    """Функцията на Розенброк"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

За по-добра визуализация, ще изобразим в 3D стойностите на функцията на Розенброк от две променливи.

Код за визуализация

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Настройка на 3D графика
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Задаване на ъгъла на виждане
ax.view_init(45, 30)

# Създаване на данни за графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))

# Изображение на повърхността
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

Предварително знаейки, че минимумът е равен на 0 при , ще разгледаме примери за това как да определим минималната стойност на функцията на Розенброк с помощта на различни процедури scipy.optimize.

Метод на симплекс на Нелдър-Мийд (Nelder-Mead)

Нека имаме начална точка x0 в 5-измерно пространство. Ще намерим най-близката до нея точка на минимума на функцията на Розенброк с помощта на алгоритъра симплекса Нелдер-Мийд (алгоритъмът е указан като стойност на параметъра method):

from scipy.optimize import minimize
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Оптимизацията завърши успешно.
         Текуща стойност на функцията: 0.000000
         Итерации: 339
         Оценки на функцията: 571
[1. 1. 1. 1. 1.]

Методът на симплекс е най-простият начин за минимизиране на явно определена и сравнително гладка функция. Не изисква изчисляване на производни на функцията, достатъчно е да зададете само нейните стойности. Методът на Нелдер-Мийд е добър избор за прости задачи по минимизация. Въпреки това, тъй като не използва оценки на градиента, откритията на минимума могат да отнемат повече време.

Метод на Пауъл

Друг алгоритъм за оптимизация, при който се изчисляват само стойности на функции, е методът на Пауъл.За да го използвате, трябва да зададете method = 'powell' в функцията minim.

x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='powell',
    options={'xtol': 1e-8, 'disp': True})
print(res.x)

Оптимизацията завърши успешно.
         Текуща стойност на функцията: 0.000000
         Итерации: 19
         Оценки на функцията: 1622
[1. 1. 1. 1. 1.]

Алгоритъм на Бройдена-Флетчера-Голдфарба-Шанно (BFGS)

За постигане на по-бързо сходство към решението, процедурата BFGS използва градиент на целевата функция. Градиентът може да бъде зададен като функция или да се изчисли с помощта на разлики от първи ред. Във всеки случай методът BFGS обикновено изисква по-малко извиквания на функции, отколкото симплекс метод.

Нека намерим производната на функцията на Розенброк в аналитичен вид:

Това изразение е валидно за производните на всички променливи, с изключение на първата и последната, които се определят като:

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

def rosen_der (x):
    xm = x [1: -1]
    xm_m1 = x [: - 2]
    xm_p1 = x [2:]
    der = np.zeros_like (x)
    der [1: -1] = 200 * (xm-xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1-xm)
    der [0] = -400 * x [0] * (x [1] -x [0] ** 2) - 2 * (1-x [0])
    der [-1] = 200 * (x [-1] -x [-2] ** 2)
    return der

Функцията за изчисление на градиента се задава като стойност на параметъра jac на функцията minim, както е показано по-долу.

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Оптимизацията успешно приключи.
         Текуща стойност на функцията: 0.000000
         Итерации: 25
         Оценки на функцията: 30
         Оценки на градиента: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

Алгоритъм на съпружените градиенти (Нютон)

Алгоритъм съпружените градиенти на Нютон е модифициран метод на Нютон.
Методът на Нютон е основан на аппроксимацията на функцията в локалната област с полином от втори ред:

където е матрицата на вторите производни (матрица на Хесен, хесиан).
Ако хесиянът е положително определен, локалният минимум на тази функция може да бъде намерен, като се приравни нулевия градиент на квадратичната форма на нула. В резултат ще се получи израз:

Обратният хесиян се изчислява с помощта на метода на съпружените градиенти. Пример за използването на този метод за минимизиране на функцията на Розенброк е показан по-долу. За да се използва методът Нютон-CG, е необходимо да се зададе функция, която изчислява хесиян.
Хесиянът на функцията на Розенброк в аналитичен вид е равен на:

където и , определят матрицата .

Останалите ненулеви елементи на матрицата са:

Например, в петмерно пространство N = 5, матрицата на Хесен за функцията на Розенброк има лентов вид:

Код, който изчислява този хесиян заедно с кода за минимизиране на функцията на Розенброк с помощта на метода на съпружените градиенти (Нютон):

def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

res = minimize(rosen, x0, method='Newton-CG', 
               jac=rosen_der, hess=rosen_hess,
               options={'xtol': 1e-8, 'disp': True})
print(res.x)

Оптимизация завершена успешно.
         Текущее значение функции: 0.000000
         Итерации: 24
         Оценки функции: 33
         Оценки градиента: 56
         Оценки гессиана: 24
[1.         1.         1.         0.99999999 0.99999999]

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

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

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

Функция, вычисляющая произведение гессиана и произвольного вектора, передается как значение аргумента hessp функции minimize:

def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] 
    -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'xtol': 1e-8, 'disp': True})

Оптимизация завершена успешно.
         Текущее значение функции: 0.000000
         Итерации: 24
         Оценки функции: 33
         Оценки градиента: 56
         Оценки гессиана: 66

Алгоритм доверительной области (trust region) сопряженных градиентов (Ньютона)

Плохая обусловленность матрицы Гессе и неверные направления поиска могут привести к тому, что алгоритм сопряженных градиентов Ньютона может быть неэффективным. В таких случаях предпочтение отдается методу доверительной области (trust-region) сопряженных градиентов Ньютона.

Пример с определением матрицы Гессе:

res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
print(res.x)

Оптимизация успешно завершена.
         Текущая функция значения: 0.000000
         Итерации: 20
         Оценки функции: 21
         Оценки градиента: 20
         Оценки гессиана: 19
[1. 1. 1. 1. 1.]

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

res = minimize(rosen, x0, method='trust-ncg', 
                jac=rosen_der, hessp=rosen_hess_p, 
                options={'gtol': 1e-8, 'disp': True})
print(res.x)

Оптимизация успешно завершена.
         Текущая функция значения: 0.000000
         Итерации: 20
         Оценки функции: 21
         Оценки градиента: 20
         Оценки гессиана: 0
[1. 1. 1. 1. 1.]

Методи от типа Крилов

Подобно на метода trust-ncg, методите от типа Крилов добре подхождат за решаване на мащабни задачи, тъй като използват само матрично-векторни произведения. Същността им е да решат проблем в доверителната област, ограничена до отсеченото подпространство на Крилов. За неопределени задачи, този метод е предпочитан, тъй като изисква по-малко нелинейни итерации за сметка на по-малкото матрично-векторни произведения за една подзадача, в сравнение с метода trust-ncg. Освен това, решението на квадратичната подзадача е по-точно отколкото при метода trust-ncg.
Пример с определением матрицы Гессе:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

Оптимизация успешно завършена.
         Текуща стойност на функцията: 0.000000
         Итерации: 19
         Оценки на функцията: 20
         Оценки на градиента: 20
         Оценки на гессиана: 18

print(res.x)

    [1. 1. 1. 1. 1.]

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

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Оптимизация успешно завършена.
         Текуща стойност на функцията: 0.000000
         Итерации: 19
         Оценки на функцията: 20
         Оценки на градиента: 20
         Оценки на гессиана: 0

print(res.x)

    [1. 1. 1. 1. 1.]

Алгоритъм за приближенное решение в доверителната област

Всички методи (Newton-CG, trust-ncg и trust-krylov) са подходящи за решаване на големи задачи (с хиляди променливи). Това е свързано с основния алгоритъм на свързаните градиенти, който предполага приблизителното намиране на обратно на гессианна матрица. Решението се намира итеративно, без явното разлагане на гессиана. Тъй като изисква само функция за произведение на гессианна матрица с произволен вектор, този алгоритъм е особено подходящ за работа с разредени (ленточни диагонални) матрици. Това осигурява ниски разходи по памет и значителна икономия на време.

В задачите със среден размер разходите за съхранение и факторизация на хесиевата матрица не са решаващи. Това означава, че може да се получи решение с по-малко итерации, като се разрешат подзадачите в областта на доверие почти точно. За тази цел някои нелинейни уравнения се решават итеративно за всяка квадратична подзадача. Такова решение обикновено изисква 3 или 4 разлагания на Холеца на хесиевата матрица. В резултат методът конвергира с по-малко итерации и изисква по-малко изчисления на целевата функция, отколкото другите реализирани методи на доверителната област. Този алгоритъм предполага единствено определяне на пълната хесиева матрица и не поддържа възможността да използва функцията за произведение на хесиева матрица и произволен вектор.

Пример за минимизиране на функцията на Розенброк:

res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Оптимизацията беше успешно прекратена.
         Текуща стойност на функцията: 0.000000
         Итерации: 13
         Оценки на функцията: 14
         Оценки на градиента: 13
         Оценки на хесиана: 14

array([1., 1., 1., 1., 1.])

На това, вероятно, ще се спрем. В следващата статия ще се опитам да разкажа най-интересното за условното минимизиране, приложението на минимизацията при решаване на задачи за аппроксимация, минимизация на функция на една променлива, произволни минимизатори и търсене на корени на система от уравнения с помощта на пакета scipy.optimize.

Източник: https://docs.scipy.org/doc/scipy/reference/

Източник: habr.com