SciPy (вымаўляецца як сай пай) - гэта пакет прыкладных матэматычных працэдур, заснаваны на пашырэнні Numpy Python. З SciPy інтэрактыўны сеанс Python ператвараецца ў такое ж паўнавартаснае асяроддзе апрацоўкі дадзеных і прататыпавання складаных сістэм, як MATLAB, IDL, Octave, R-Lab і SciLab. Сёння я хачу сцісла распавесці аб тым, як варта ўжываць некаторыя вядомыя алгарытмы аптымізацыі ў пакеце scipy.optimize. Больш падрабязную і актуальную даведку па ўжыванні функцый заўсёды можна атрымаць з дапамогай каманды help() ці з дапамогай Shift+Tab.
Увядзенне
Каб пазбавіць самога сябе і чытачоў ад пошуку і чытання першакрыніц, спасылкі на апісанні метадаў будуць у асноўным на вікіпедыю. Як правіла, гэтай інфармацыі дастаткова для разумення метадаў у агульных рысах і ўмоў іх прымянення. Для разумення сутнасці матэматычных метадаў ідзем па спасылках на больш аўтарытэтныя публікацыі, якія можна знайсці ў канцы кожнага артыкула ці ў каханай пошукавай сістэме.
Такім чынам, модуль scipy.optimize уключае ў сябе рэалізацыю наступных працэдур:
- Умоўнай і безумоўнай мінімізацыі скалярных функцый некалькіх зменных (minim) з дапамогай розных алгарытмаў (сімплекс Нэлдэра-Міда, BFGS, спалучаных градыентаў Ньютана, и )
- Глабальнай аптымізацыі (напрыклад: , )
- Мінімізацыя рэшткаў (least_squares) і алгарытмы падганяння крывых нелінейным МНК (curve_fit)
- Мінімізацыі скалярнай функцый адной зменнай (minim_scalar) і пошуку каранёў (root_scalar)
- Шматмерныя решатели сістэмы раўнанняў (root) з выкарыстаннем розных алгарытмаў (гібрыдны Паўэла, ці буйнамаштабныя метады, такія як ).
У гэтым артыкуле мы разгледзім толькі першы пункт з усяго гэтага сьпісу.
Безумоўная мінімізацыя скалярнай функцыі некалькіх зменных
Функцыя 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 пры , разгледзім прыклады таго, як вызначыць мінімальнае значэнне функцыі Розенброка з дапамогай розных працэдур 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 339
Function evaluations: 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 1622
[1. 1. 1. 1. 1.]Алгарытм Бройдэна-Флетчара-Голдфарба-Шанно (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 функцыі 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})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 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)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.xOptimization 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.
Крыніца:
Крыніца: habr.com