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

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

Вовед

Со цел да се спасите себеси и читателите од пребарување и читање примарни извори, врските до описите на методите ќе бидат главно на Википедија. Како по правило, овие информации се доволни за да се разберат методите во општи услови и условите за нивна примена. За да ја разберете суштината на математичките методи, следете ги врските до повеќе авторитетни публикации, кои може да се најдат на крајот од секоја статија или во вашиот омилен пребарувач.

Значи, модулот scipy.optimize вклучува имплементација на следниве постапки:

1. Условно и безусловно минимизирање на скаларните функции на неколку променливи (минимум) со користење на различни алгоритми (Nelder-Mead simplex, BFGS, Њутн конјугирани градиенти, КОБИЛА и SLSQP)
2. Глобална оптимизација (на пример: басенот, diff_evolution)
3. Минимизирање на остатоците МНЦ (најмалку_квадрати) и алгоритми за фитинг на крива користејќи нелинеарни најмали квадрати (крива_фит)
4. Минимизирање на скаларните функции на една променлива (minim_scalar) и пребарување на корени (root_scalar)
5. Повеќедимензионални решавачи на систем на равенки (корен) со користење на различни алгоритми (хибриден Пауел, Левенберг-Маркарт или методи од големи размери како што се Њутн-Крилов).

Во оваа статија ќе ја разгледаме само првата ставка од целата оваа листа.

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

Функцијата minim од пакетот scipy.optimize обезбедува општ интерфејс за решавање на проблеми со условно и безусловно минимизирање на скаларните функции на неколку променливи. За да покажеме како функционира, ќе ни треба соодветна функција од неколку променливи, кои ќе ги минимизираме на различни начини.

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

И покрај фактот дека функцијата Rosenbrock и нејзините Jacobi и Hessian матрици (првиот и вториот дериват, соодветно) се веќе дефинирани во пакетот 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 во , ајде да погледнеме примери за тоа како да се одреди минималната вредност на функцијата Розенброк користејќи различни процедури scipy.optimize.

Нелдер-Мид симплекс метод

Нека има почетна точка x0 во 5-димензионален простор. Ајде да ја најдеме минималната точка на функцијата Розенброк најблиску до неа користејќи го алгоритмот Нелдер-Мид симплекс (алгоритмот е наведен како вредност на параметарот на методот):

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 е добар избор за едноставни проблеми со минимизирање. Меѓутоа, бидејќи не користи проценки на градиент, може да биде потребно подолго време да се најде минимумот.

Пауел метод

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

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.]

Алгоритам Broyden-Fletcher-Goldfarb-Shanno (BFGS).

За да се добие побрзо конвергенција до решение, постапката BFGS користи градиент на функцијата цел. Градиентот може да се одреди како функција или да се пресмета со користење на разлики од прв ред. Во секој случај, методот BFGS обично бара помалку повици на функции од методот симплекс.

Да го најдеме дериватот на функцијата Розенброк во аналитичка форма:

Овој израз е валиден за изводите на сите променливи освен првата и последната, кои се дефинирани како:

Ајде да ја погледнеме функцијата Пајтон што го пресметува овој градиент:

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 на функцијата минимална, како што е прикажано подолу.

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

Алгоритам за доверлив регион на конјугиран градиент (Њутн)

Лошата кондиција на Хесијанската матрица и неточните насоки за пребарување може да предизвикаат Њутновиот коњугатен градиент алгоритам да биде неефикасен. Во такви случаи, предност се дава на метод на доверлив регион (доверба-регион) конјугирани Њутнови градиенти.

Пример со дефиниција на Хесијанската матрица:

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