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

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

въведение

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

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

1. Условно и безусловно минимизиране на скаларни функции на няколко променливи (минимум) с помощта на различни алгоритми (симплекс на Nelder-Mead, BFGS, конюгирани градиенти на Нютон, КОБИЛА и SLSQP)
2. Глобална оптимизация (например: басейнови хълмове, diff_evolution)
3. Минимизиране на остатъците MNC (least_squares) и нелинейни алгоритми за напасване на кривата на най-малките квадрати (curve_fit)
4. Минимизиране на скаларни функции на една променлива (minim_scalar) и намиране на корени (root_scalar)
5. Решаватели на системи с многовариантни уравнения (root), използващи различни алгоритми (Hybrid Powell, Левенберг-Марквард или широкомащабни методи като напр Нютон-Крилов).

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

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

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

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

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

import numpy as np

def rosen(x):
    """The Rosenbrock function"""
    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 при , нека да разгледаме примери как да определим минималната стойност на функцията Rosenbrock с помощта на различни scipy.optimize процедури.

Симплексен метод на Nelder-Mead (Nelder-Mead)

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

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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 339
         Function evaluations: 571
[1. 1. 1. 1. 1.]

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

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

Друг алгоритъм за оптимизация, при който се изчисляват само стойностите на функцията, е Метод на Пауъл. За да го използвате, трябва да зададете 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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 25
         Function evaluations: 30
         Gradient evaluations: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

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

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

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

Обратният хесиан се изчислява с помощта на метода на спрегнатия градиент. Пример за използване на този метод за минимизиране на функцията на Розенброк е даден по-долу. За да използвате метода Newton-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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 24
[1.         1.         1.         0.99999999 0.99999999]

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

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

Да разгледаме функцията hess, която приема вектор за минимизиране като първи аргумент и произволен вектор като втори аргумент (заедно с други аргументи на функцията, която трябва да бъде минимизирана). В нашия случай не е много трудно да се изчисли произведението на хесиана на функцията на Розенброк с произволен вектор. Ако p е произволен вектор, тогава продуктът има формата:

Функцията, която изчислява произведението на Хесиан и произволен вектор, се предава като стойност на аргумента hessp към функцията за минимизиране:

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})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 24
         Function evaluations: 33
         Gradient evaluations: 56
         Hessian evaluations: 66

Алгоритъм за доверен регион за конюгирани градиенти (Нютон)

Лошата условност на матрицата на Хесиан и неправилните посоки на търсене могат да доведат до неефективността на алгоритъма на конюгирания градиент на Нютон. В такива случаи се дава предимство метод на доверен регион (trust-region) Нютонови спрегнати градиенти.

Пример с дефиницията на матрицата на Хесиан:

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

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 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)

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 21
         Gradient evaluations: 20
         Hessian evaluations: 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})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 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})

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 19
         Function evaluations: 20
         Gradient evaluations: 20
         Hessian evaluations: 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

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 13
         Function evaluations: 14
         Gradient evaluations: 13
         Hessian evaluations: 14

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

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

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

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