SciPy (vyslovuje se jako saj paj) je balíček aplikačních matematických procedur založený na rozšíření Numpy Python. Se SciPy se interaktivní relace Pythonu mění na plnohodnotné prostředí pro zpracování dat a prototypování složitých systémů, jako je MATLAB, IDL, Octave, R-Lab a SciLab. Dnes bych chtěl krátce přiblížit, jak by se měly používat některé známé optimalizační algoritmy v balíčku scipy.optimize. Podrobnější a aktuální dokumentaci o použití funkcí můžete vždy získat pomocí příkazu help() nebo pomocí Shift+Tab.
Úvod
Abychom ušetřili sebe i čtenáře od hledání a čtení primárních zdrojů, odkazy na popisy metod budou převážně na Wikipedii. Obvykle tyto informace stačí k pochopení metod v obecných rysech a podmínek jejich použití. Pro pochopení podstaty matematických metod pokračujeme odkazy na autoritativnější publikace, které najdete na konci každého článku nebo ve své oblíbené vyhledávačce.
Modul scipy.optimize obsahuje implementaci následujících procedur:
- Podmíněné a bezpodmínečné minimalizace skalárních funkcí vícero proměnných (minim) pomocí různých algoritmů (simplex Neldera-Meada, BFGS, konjugovaných gradientů Newtona, a )
- Globální optimalizace (například: , )
- Minimalizace zbytků (least_squares) a algoritmy pro přizpůsobení křivek nelineárním MNK (curve_fit)
- Minimalizace skalárních funkcí jedné proměnné (minim_scalar) a hledání kořenů (root_scalar)
- Více dimenzionální řešitelé systémů rovnic (root) s využitím různých algoritmů (hybridní Powell, nebo velkoplošné metody, jako je ).
V tomto článku se zaměříme pouze na první bod z tohoto seznamu.
Bezpodmínečná minimalizace skalární funkce vícero proměnných
Funkce minim z balíčku scipy.optimize poskytuje obecné rozhraní pro řešení úloh podmíněné a bezpodmínečné minimalizace skalárních funkcí vícero proměnných. K demonstraci jejího fungování budeme potřebovat vhodnou funkci vícero proměnných, kterou budeme různě minimalizovat.
K těmto účelům se skvěle hodí Rosenbrockova funkce N proměnných, která má tvar:
Hoci je funkcia Rosenbrock a jej Jacobiho a Hessova matica (prvá a druhá derivácia) už definovaná v balíku scipy.optimize, definujme ju samostatne.
import numpy as np
def rosen(x):
"""Funkcia Rosenbrock"""
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 funkcie Rosenbrock od dvoch premenných.
Kód na vykreslenie
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# Nastavujeme 3D graf
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Nastavujeme uhol pohľadu
ax.view_init(45, 30)
# Vytvárame dáta pre graf
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]))
# Vykresľujeme plochu
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
Keď už vieme, že minimum je 0 pri 
, zvážime príklady toho, ako určiť minimum funkcie Rosenbrock pomocou rôznych procedúr scipy.optimize.
Simpliciálny metóda Neldera-Meada (Nelder-Mead)
Nech máme počiatočný bod x0 v 5-rozmernom priestore. Nájdeme najbližší bod minima funkcie Rosenbrock pomocou algoritmu (algoritmus uvedený ako hodnota parametra 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)Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 339
Vyhodnotenia funkcie: 571
[1. 1. 1. 1. 1.]Simpliciálny metóda je najjednoduchší spôsob, ako minimalizovať explicitne definovanú a celkom hladkú funkciu. Nevyžaduje výpočty derivácií funkcie, stačí zadať len jej hodnoty. Metóda Neldera-Meada je dobrá voľba pre jednoduché úlohy minimalizácie. Avšak, pretože nepoužíva odhady gradientu, nájdenie minima môže trvať dlhšie.
Metóda Powell
Ďalším algoritmom optimalizácie, pri ktorom sa počítajú iba hodnoty funkcií, je Na jej použitie je potrebné nastaviť method = ‘powell’ vo funkcii 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)Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 19
Vyhodnotenia funkcie: 1622
[1. 1. 1. 1. 1.]Algoritmus Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Na dosiahnutie rýchlejšej konvergencie k riešeniu, procedúra používa gradient cieľovej funkcie. Gradient môže byť daný ako funkcia alebo sa môže vypočítať pomocou prvého poriadku. V každom prípade metóda BFGS zvyčajne vyžaduje menej volaní funkcií ako simplexová metóda.
Nájdeme deriváciu funkcie Rosenbrock v analytickom tvare:
Tento výraz platí pre derivácie všetkých premenných, okrem prvej a poslednej, ktoré sú definované ako:
Pozrime sa na funkciu Python, 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 na výpočet gradientu sa uvádza ako hodnota parametru jac funkcie minim, ako je znázornené nižšie.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 25
Vyhodnotenia funkcie: 30
Vyhodnotenia gradientu: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Algoritmus konjugovaných gradientov (Newtona)
Algoritmus je modifikovaná metóda Newtona.
Metóda Newtona je založená na aproximácii funkcie v lokálnej oblasti polynomom druhého stupňa:
kde 
je matica druhých derivácií (matica Hessiana, hessian).
Ak je hessian pozitívne definitný, lokálny minimum tejto funkcie možno nájsť, keď sa nulový gradient kvadratickej formy rovná nule. Výsledkom je výraz:
Obratený hessian sa vypočítava pomocou metódy konjugovaných gradientov. Príklad použitia tejto metódy na minimalizáciu funkcie Rosenbrock je uvedený nižšie. Na použitie metódy Newton-CG je potrebné zadať funkciu, ktorá vypočítava hessian.
Hessian funkcie Rosenbrock v analytickom tvare sa rovná:
kde 
a 
, určujú maticu 
.
Ostatné nenulové prvky matice sú rovné:
Napríklad, v päťrozmernom priestore N = 5, matica Hessiana pre funkciu Rosenbrock má pásový tvar:
Kód, ktorý vypočítava tento hessian spolu s kódom na minimalizáciu funkcie Rosenbrock pomocou metódy konjugovaných gradientov (Newtona):
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)Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 24
Vyhodnotenia funkcie: 33
Vyhodnotenia gradientu: 56
Vyhodnotenia hessianu: 24
[1. 1. 1. 0.99999999 0.99999999]Príklad s definíciou funkcie súčinu hessiana a ľubovoľného vektora
Pri reálnych úlohách môže výpočet a uchovávanie celej matice Hessian vyžadovať značné množstvo času a pamäte. V skutočnosti nie je potrebné zadávať samotnú maticu Hessian, pretože na minimalizačný postup potrebujeme iba vektor, ktorý je súčinom hessiana s iným ľubovoľným vektorom. Tým pádom je z výpočtového hľadiska oveľa vhodnejšie definovať funkciu, ktorá vracia výsledok súčinu hessiana s ľubovoľným vektorom.
Zvážime funkciu hess, ktorá prijíma vektor minimalizácie ako prvý argument a ľubovoľný vektor ako druhý argument (spolu s ďalšími argumentami minimalizovanej funkcie). V našom prípade nie je ťažké vypočítať súčin hessiana funkcie Rosenbrockovej s ľubovoľným vektorom. Ak p je ľubovoľný vektor, tak súčin 
má tvar:
Funkcia, ktorá vypočítava súčin hessiana a ľubovoľného vektora, sa predáva ako hodnota argumentu hessp pre funkciu 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})
Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 24
Vyhodnotenia funkcie: 33
Vyhodnotenia gradientu: 56
Vyhodnotenia hessianu: 66Algoritmus dôveryhodného regiónu (trust region) spojených gradientov (Newtona)
Zlá podmienka matice Hessianu a nesprávne smerovania hľadania môžu spôsobiť, že algoritmus spojených gradientov Newtona môže byť neefektívny. V takýchto prípadoch sa uprednostňuje (trust-region) spojených gradientov Newtona.
Príklad s definíciou matice Hessian:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 20
Vyhodnotenia funkcie: 21
Vyhodnotenia gradientu: 20
Vyhodnotenia Hessovej matice: 19
[1. 1. 1. 1. 1.]Príklad s funkciou produktu Hessovej matice a ľubovoľného vektora:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 20
Vyhodnotenia funkcie: 21
Vyhodnotenia gradientu: 20
Vyhodnotenia Hessovej matice: 0
[1. 1. 1. 1. 1.]Metódy Krylovovho typu
Podobne ako metóda trust-ncg, metódy Krylovovho typu sú vhodné na riešenie veľkoplošných úloh, pretože používajú iba maticovo-vektorové súčiny. Ich podstata spočíva v riešení úlohy v dôveryhodnej oblasti, obmedzenej skráteným Krylovovým podpriestorom. Pre neurčité úlohy je lepšie použiť túto metódu, pretože vyžaduje menej nelineárnych iterácií vďaka menšiemu počtu maticovo-vektorových súčinov pre jednu podúlohu, v porovnaní s metódou trust-ncg. Okrem toho sa kvadratická podúloha nachádza presnejšie ako v metóde trust-ncg.
Príklad s definíciou matice Hessian:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 19
Vyhodnotenia funkcie: 20
Vyhodnotenia gradientu: 20
Vyhodnotenia Hessovej matice: 18
print(res.x)
[1. 1. 1. 1. 1.]
Príklad s funkciou produktu Hessovej matice a ľubovoľného vektora:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Optimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 19
Vyhodnotenia funkcie: 20
Vyhodnotenia gradientu: 20
Vyhodnotenia Hessovej matice: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algoritmus aproximovaného riešenia v dôveryhodnej oblasti
Všetky metódy (Newton-CG, trust-ncg a trust-krylov) sú vhodné na riešenie veľkoplošných úloh (s tisíckami premenných). To vyplýva z toho, že základný algoritmus konjugovaných gradientov predpokladá aproximované nájdenie inverznej Hessovej matice. Riešenie sa nachádza iteratívne, bez explicitného rozkladu Hessovej matice. Keďže je potrebné určiť iba funkciu pre súčin Hessovej matice a ľubovoľného vektora, tento algoritmus je obzvlášť dobrý na prácu s riedkymi (pásikovými diagonálnymi) maticami. To zabezpečuje nízke náklady na pamäť a značnú úsporu času.
Pri úlohách strednej veľkosti náklady na uchovávanie a faktorovanie Hessovej matice nemajú rozhodujúci význam. To znamená, že je možné získať riešenie za menej iterácií, ak sa podúlohy dôveryhodnej oblasti vyriešia takmer presne. Na tento účel sa niektoré nelineárne rovnice riešia iteratívne pre každú kvadratickú podúlohu. Takéto riešenie zvyčajne vyžaduje 3 alebo 4 rozklady Choleského Hessovej matice. Výsledkom je, že metóda konverguje za menej iterácií a vyžaduje menej výpočtov cieľovej funkcie ako iné realizované metódy dôveryhodnej oblasti. Tento algoritmus predpokladá iba určenie celej Hessovej matice a nepodporuje možnosť používať funkciu súčinu Hessovej matice a ľubovoľného vektora.
Príklad s minimalizáciou funkcie Rosenbrock:
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimalizácia bola úspešne ukončená.
Aktuálna hodnota funkcie: 0.000000
Iterácie: 13
Vyhodnotenia funkcie: 14
Vyhodnotenia gradientu: 13
Vyhodnotenia Hessovej matice: 14
array([1., 1., 1., 1., 1.])Na tomto sa asi zastavíme. V nasledujúcom článku sa pokúsim povedať najzaujímavejšie o podmienkovej minimalizácii, aplikácii minimalizácie pri riešení úloh aproximácie, minimalizácii funkcie jednej premennej, ľubovoľných minimalizátoroch a hľadaní koreňov systému rovníc pomocou balíka scipy.optimize.
Zdroj:
Zdroj: habr.com
