SciPy (vyslovuje sa sai koláč) je matematický balík aplikácií založený na rozšírení Numpy Python. So SciPy sa vaša interaktívna Python relácia stane rovnakým prostredím pre kompletnú vedu o údajoch a komplexné systémové prototypové prostredie ako MATLAB, IDL, Octave, R-Lab a SciLab. Dnes chcem stručne hovoriť o tom, ako používať niektoré známe optimalizačné algoritmy v balíku scipy.optimize. Podrobnejšiu a najaktuálnejšiu pomoc o používaní funkcií môžete vždy získať pomocou príkazu help() alebo pomocou Shift+Tab.
Úvod
Aby ste seba a čitateľov ušetrili pred hľadaním a čítaním primárnych zdrojov, odkazy na popisy metód budú hlavne na Wikipédii. Tieto informácie spravidla postačujú na pochopenie metód vo všeobecných podmienkach a podmienok ich aplikácie. Aby ste pochopili podstatu matematických metód, postupujte podľa odkazov na autoritatívnejšie publikácie, ktoré nájdete na konci každého článku alebo vo svojom obľúbenom vyhľadávači.
Modul scipy.optimize teda zahŕňa implementáciu nasledujúcich postupov:
- Podmienená a nepodmienená minimalizácia skalárnych funkcií viacerých premenných (minimálne) pomocou rôznych algoritmov (Nelder-Mead simplex, BFGS, Newtonove konjugované gradienty, и )
- Globálna optimalizácia (napríklad: , )
- Minimalizácia zvyškov (najmenší_štvorce) a algoritmy prekladania kriviek pomocou nelineárnych najmenších štvorcov (prispôsobenie krivke)
- Minimalizácia skalárnych funkcií jednej premennej (minim_scalar) a hľadanie koreňov (root_scalar)
- Viacrozmerní riešitelia sústavy rovníc (koreň) využívajúci rôzne algoritmy (hybridný Powell, alebo metódy vo veľkom meradle, ako napr ).
V tomto článku sa budeme zaoberať iba prvou položkou z celého tohto zoznamu.
Bezpodmienečná minimalizácia skalárnej funkcie viacerých premenných
Funkcia minim z balíka scipy.optimize poskytuje všeobecné rozhranie na riešenie problémov podmienenej a nepodmienenej minimalizácie skalárnych funkcií viacerých premenných. Na demonštráciu jej fungovania budeme potrebovať vhodnú funkciu viacerých premenných, ktoré budeme rôznymi spôsobmi minimalizovať.
Na tieto účely je dokonalá Rosenbrockova funkcia N premenných, ktorá má tvar:
Napriek tomu, že Rosenbrockova funkcia a jej Jacobiho a Hessova matica (prvá a druhá derivácia) sú už definované v balíku scipy.optimize, zadefinujeme si ju sami.
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)Pre prehľadnosť nakreslíme v 3D hodnoty Rosenbrockovej funkcie dvoch premenných.
Kód kreslenia
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()
Vopred vediac, že minimum je 0 at , pozrime sa na príklady, ako určiť minimálnu hodnotu Rosenbrockovej funkcie pomocou rôznych postupov scipy.optimize.
Nelder-Mead simplexná metóda
Nech je počiatočný bod x0 v 5-rozmernom priestore. Pomocou algoritmu nájdime minimálny bod Rosenbrockovej funkcie, ktorý je k nej najbližšie (algoritmus je špecifikovaný ako hodnota parametra metódy):
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.]Simplexová metóda je najjednoduchší spôsob, ako minimalizovať explicitne definovanú a pomerne hladkú funkciu. Nevyžaduje výpočet derivácií funkcie, stačí zadať len jej hodnoty. Metóda Nelder-Mead je dobrou voľbou pre jednoduché problémy s minimalizáciou. Keďže však nepoužíva odhady gradientu, nájdenie minima môže trvať dlhšie.
Powellova metóda
Ďalším optimalizačným algoritmom, v ktorom sa počítajú iba funkčné hodnoty, je . Ak ju chcete použiť, musíte vo funkcii minim nastaviť metódu = '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.]Algoritmus Broyden-Fletcher-Goldfarb-Shanno (BFGS).
Ak chcete získať rýchlejšiu konvergenciu k riešeniu, postup využíva gradient účelovej funkcie. Gradient môže byť špecifikovaný ako funkcia alebo vypočítaný pomocou rozdielov prvého rádu. V každom prípade metóda BFGS zvyčajne vyžaduje menej volaní funkcií ako simplexná metóda.
Nájdite deriváciu Rosenbrockovej funkcie v analytickej forme:
Tento výraz je platný pre deriváty všetkých premenných okrem prvej a poslednej, ktoré sú definované ako:
Pozrime sa na funkciu Pythonu, ktorá vypočítava tento gradient:
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 derFunkcia výpočtu gradientu je špecifikovaná ako hodnota parametra jac funkcie minima, ako je uvedené nižšie.
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]Algoritmus konjugovaného gradientu (Newton)
Algoritmus je modifikovaná Newtonova metóda.
Newtonova metóda je založená na aproximácii funkcie v lokálnej oblasti polynómom druhého stupňa:
kde je matica druhých derivácií (Hessova matica, Hessova).
Ak je Hessián kladne definitný, potom lokálne minimum tejto funkcie možno nájsť prirovnaním nulového gradientu kvadratickej formy k nule. Výsledkom bude výraz:
Inverzný Hessián sa vypočíta pomocou metódy konjugovaného gradientu. Príklad použitia tejto metódy na minimalizáciu Rosenbrockovej funkcie je uvedený nižšie. Ak chcete použiť metódu Newton-CG, musíte zadať funkciu, ktorá vypočíta Hessian.
Hessián Rosenbrockovej funkcie v analytickej forme sa rovná:
kde и , definujte maticu .
Zostávajúce nenulové prvky matice sa rovnajú:
Napríklad v päťrozmernom priestore N = 5 má Hessova matica pre Rosenbrockovu funkciu tvar pásu:
Kód, ktorý vypočíta tento Hessian spolu s kódom na minimalizáciu Rosenbrockovej funkcie pomocou metódy konjugovaného gradientu (Newton):
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]Príklad s definíciou funkcie súčinu Hessian a ľubovoľného vektora
Pri problémoch v reálnom svete môže výpočet a ukladanie celej Hessovej matice vyžadovať značné časové a pamäťové zdroje. V tomto prípade vlastne nie je potrebné špecifikovať samotnú Hessovu maticu, pretože minimalizačná procedúra vyžaduje len vektor rovný súčinu Hessianu s iným ľubovoľným vektorom. Z výpočtového hľadiska je teda oveľa vhodnejšie okamžite definovať funkciu, ktorá vráti výsledok súčinu Hessian s ľubovoľným vektorom.
Zvážte funkciu hess, ktorá berie minimalizačný vektor ako prvý argument a ľubovoľný vektor ako druhý argument (spolu s ďalšími argumentmi funkcie, ktorá sa má minimalizovať). V našom prípade nie je veľmi ťažké vypočítať súčin Hessiánskej Rosenbrockovej funkcie s ľubovoľným vektorom. Ak p je ľubovoľný vektor, potom súčin vyzerá ako:
Funkcia, ktorá vypočítava súčin Hessian a ľubovoľného vektora, sa odovzdá ako hodnota argumentu hessp do funkcie minimalizácie:
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: 66Algoritmus konjugovanej oblasti dôveryhodnosti gradientu (Newton)
Slabá úprava Hessovej matice a nesprávne smery vyhľadávania môžu spôsobiť, že Newtonov algoritmus konjugovaného gradientu bude neúčinný. V takýchto prípadoch sa dáva prednosť (trust-region) konjugovať Newtonove gradienty.
Príklad s definíciou Hessovej matice:
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.]Príklad s funkciou súčinu Hessian a ľubovoľným vektorom:
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.]Metódy typu Krylov
Podobne ako metóda trust-ncg, metódy Krylovovho typu sú vhodné na riešenie rozsiahlych problémov, pretože používajú iba produkty maticového vektora. Ich podstatou je vyriešiť problém v oblasti dôvery obmedzenej skráteným Krylovovým podpriestorom. Pre neisté problémy je lepšie použiť túto metódu, pretože používa menší počet nelineárnych iterácií v dôsledku menšieho počtu produktov maticového vektora na podproblém v porovnaní s metódou trust-ncg. Okrem toho sa riešenie kvadratického podproblému nájde presnejšie ako pomocou metódy trust-ncg.
Príklad s definíciou Hessovej matice:
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.]
Príklad s funkciou súčinu Hessian a ľubovoľným vektorom:
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.]
Algoritmus na približné riešenie v oblasti spoľahlivosti
Všetky metódy (Newton-CG, trust-ncg a trust-krylov) sú vhodné na riešenie rozsiahlych problémov (s tisíckami premenných). Je to spôsobené skutočnosťou, že základný algoritmus konjugovaného gradientu predpokladá približné určenie inverznej Hessovej matice. Riešenie sa nachádza iteračne, bez explicitného rozšírenia hesiančiny. Keďže potrebujete definovať funkciu iba pre súčin hessiánskeho a ľubovoľného vektora, tento algoritmus je vhodný najmä na prácu s riedkymi (pásmovými diagonálnymi) maticami. To poskytuje nízke náklady na pamäť a výraznú úsporu času.
Pri stredne veľkých problémoch nie sú náklady na skladovanie a faktoring Hessian kritické. To znamená, že je možné získať riešenie v menšom počte opakovaní, čím sa riešia podproblémy dôveryhodnej oblasti takmer presne. Na tento účel sa niektoré nelineárne rovnice riešia iteračne pre každý kvadratický podproblém. Takéto riešenie zvyčajne vyžaduje 3 alebo 4 Choleského rozklady Hessovej matice. Výsledkom je, že metóda konverguje v menšom počte opakovaní a vyžaduje menej výpočtov objektívnej funkcie ako iné implementované metódy oblasti spoľahlivosti. Tento algoritmus zahŕňa iba určenie úplnej Hessovej matice a nepodporuje možnosť použiť funkciu súčinu Hessovej matice a ľubovoľného vektora.
Príklad s minimalizáciou Rosenbrockovej funkcie:
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.])Asi sa tam zastavíme. V ďalšom článku sa pokúsim povedať to najzaujímavejšie o podmienenej minimalizácii, aplikácii minimalizácie pri riešení aproximačných úloh, minimalizácii funkcie jednej premennej, ľubovoľných minimalizéroch a hľadaní koreňov sústavy rovníc pomocou scipy.optimize balík.
Zdroj:
Zdroj: hab.com