SciPy (čte se jako saj paj) je balíček aplikovaných matematických procedur založený na rozšíření Numpy Python. Se SciPy se interaktivní sezení Pythonu mění v plnohodnotné prostředí zpracování dat a prototypování složitých systémů jako MATLAB, IDL, Octave, R-Lab a SciLab. Dnes bych chtěl stručně přiblížit, jak by se měly uplatňovat některé známé algoritmy optimalizace v balíčku scipy.optimize. Podrobnější a aktuální dokumentaci k používání funkcí lze vždy získat pomocí příkazu help() nebo kombinací Shift+Tab.
Úvod
Abych se zbavil sebe a čtenářů hledání a čtení původních zdrojů, odkazy na popisy metod budou většinou na Wikipedii. Obecně lze říci, že tyto informace jsou dostatečné pro pochopení metod v obecných obrysech a podmínek jejich použití. Pro pochopení podstaty matematických metod odkazujeme na autoritativnější publikace, které lze najít na konci každého článku nebo v oblíbeném vyhledávači.
Modul scipy.optimize zahrnuje implementaci následujících procedur:
- Podmíněné a nepodmíněné minimalizace skalárních funkcí více 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 fittingu křivek nelineární MNK (curve_fit)
- Minimalizace skalární funkce jedné proměnné (minim_scalar) a hledání kořenů (root_scalar)
- Vícedimenzionální řešitelé systémů rovnic (root) s použ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.
Nepodmíněná minimalizace skalární funkce více proměnných
Funkce minim z balíčku scipy.optimize poskytuje obecné rozhraní pro řešení úloh podmíněné a nepodmíněné minimalizace skalárních funkcí více proměnných. Abychom ukázali její fungování, potřebujeme vhodnou funkci více proměnných, kterou budeme různě minimalizovat.
Pro tyto účely se skvěle hodí Rosenbrockova funkce s N proměnnými, která má tvar:
I když je Rosenbrockova funkce a její Jacobian a Hessian (první a druhá derivace) již definována v balíčku scipy.optimize, definujme ji sami.
import numpy as np
def rosen(x):
"""Rosenbrockova funkce"""
return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)Abychom to zjednodušili, vykreslíme 3D hodnoty Rosenbrockovy funkce v závislosti na dvou proměnných.
Kód pro vykreslení
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 úhel pohledu
ax.view_init(45, 30)
# Vytváříme data pro 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]))
# Vykreslíme plochu
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
S vědomím, že minimum je rovno 0 při 
, se podíváme na příklady toho, jak určit minimální hodnotu Rosenbrockovy funkce pomocí různých procedur scipy.optimize.
Metoda simplex Neldera-Meada (Nelder-Mead)
Předpokládejme, že máme počáteční bod x0 v 5-dimenzionálním prostoru. Najděme nejbližší bod minima Rosenbrockovy funkce pomocí algoritmu (algoritmus uvedený jako hodnota parametru 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)Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 339
Vyhodnocení funkce: 571
[1. 1. 1. 1. 1.]Metoda simplex je nejjednodušší způsob, jak minimalizovat explicitně definovanou a poměrně hladkou funkci. Nevyžaduje výpočet derivací funkce, stačí zadat pouze její hodnoty. Metoda Neldera-Meada je dobrou volbou pro jednoduché úlohy minimalizace. Nicméně protože nevyužívá odhady gradientu, může to trvat déle, než se minimum nalézá.
Metoda Powella
Dalším optimalizačním algoritmem, který počítá pouze hodnoty funkcí, je Chcete-li ji použít, je třeba nastavit method = 'powell' ve funkci 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)Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 19
Vyhodnocení funkce: 1622
[1. 1. 1. 1. 1.]Algoritmus Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Pro rychlejší konvergenci k řešení, procedura používá gradient cílové funkce. Gradient může být definován jako funkce nebo se vypočítává pomocí prvního ordinálu. V každém případě metoda BFGS obvykle vyžaduje méně volání funkcí než simplexní metoda.
Najdeme analytickou derivaci od Rosenbrockovy funkce:
Toto vyjádření platí pro derivace všech proměnných kromě první a poslední, které jsou definovány jako:
Podívejme se na funkci v Pythonu, která tento gradient počítá:
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 derFunkce pro výpočet gradientu se zadává jako hodnota parametru jac funkce minim, jak je ukázáno níže.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 25
Vyhodnocení funkcí: 30
Vyhodnocení gradientu: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Algoritmus Conjugate Gradient (Newton)
Algoritmus je modifikovanou metodou Newtona.
Newtonova metoda vychází z aproximace funkce v místní oblasti polynomem druhého stupně:
kde 
je matice druhých derivací (Hesseova matice, hessian).
Pokud je hessian pozitivně definitní, místní minimum této funkce lze nalézt vyrovnáním nulového gradientu kvadratické formy na nulu. Výsledkem bude vyjádření:
Inverzní hessian se počítá pomocí metody konjugovaných gradientů. Příklad použití této metody pro minimalizaci Rosenbrockovy funkce je uveden níže. Pro použití metody Newton-CG je třeba zadat funkci, která počítá hessian.
Hessian funkce Rosenbrock v analytické podobě je:
kde 
a 
, určuje matici 
.
Ostatní nenulové prvky matice jsou:
Například v pětirozměrném prostoru N = 5 má Hesseova matice pro Rosenbrockovu funkci pásový tvar:
Kód, který počítá tento hessian spolu s kódem pro minimalizaci Rosenbrockovy funkce pomocí metody konjugovaných gradientů (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)Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 24
Vyhodnocení funkce: 33
Vyhodnocení gradientu: 56
Vyhodnocení Hessovy matice: 24
[1. 1. 1. 0.99999999 0.99999999]Příklad s definicí funkce pro součin Hessovy matice a libovolného vektoru
V reálných úlohách může výpočet a uchovávání celé Hessovy matice vyžadovat značné časové a paměťové zdroje. Ve skutečnosti však není nutné definovat samotnou Hessovu matici, neboť k procesu minimalizace je potřeba pouze vektor, který je roven součinu Hessovy matice s jiným libovolným vektorem. Z pohledu výpočetní efektivity je tedy mnohem výhodnější ihned definovat funkci, která vrací výsledek součinu Hessovy matice s libovolným vektorem.
Zvažme funkci hess, která přijímá vektor minimalizace jako první argument, a libovolný vektor jako druhý argument (spolu s dalšími argumenty minimalizované funkce). V našem případě není příliš složité spočítat součin Hessovy matice funkce Rosenbrok s libovolným vektorem. Jestliže p — libovolný vektor, pak součin 
má tvar:
Funkce, která počítá součin Hessovy matice a libovolného vektoru, je předána jako hodnota argumentu hessp funkce 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})
Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 24
Vyhodnocení funkce: 33
Vyhodnocení gradientu: 56
Vyhodnocení Hessovy matice: 66Algoritmus s důvěryhodnou oblastí (trust region) pro konjugované gradienty (Newtonovy)
Špatná kondicionalita Hessovy matice a nesprávné směry hledání mohou vést k neefektivnosti algoritmu konjugovaných gradientů Newtona. V takových případech se upřednostňuje (trust-region) konjugovaných gradientů Newtona.
Příklad s definicí Hessovy matice:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 20
Vyhodnocení funkce: 21
Vyhodnocení gradientu: 20
Vyhodnocení hessiana: 19
[1. 1. 1. 1. 1.]Příklad s funkcí produktu hessiana a libovolného vektoru:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 20
Vyhodnocení funkce: 21
Vyhodnocení gradientu: 20
Vyhodnocení hessiana: 0
[1. 1. 1. 1. 1.]Metody typu Krylov
Podobně jako metoda trust-ncg, metody typu Krylov jsou vhodné pro řešení rozsáhlých úloh, protože využívají pouze matice-vektorové součiny. Jejich podstatou je řešení úlohy v důvěrné oblasti, omezené na zúžený Krylovův prostor. Pro neurčité úlohy je lepší tento metod, protože vyžaduje méně nelineárních iterací díky menšímu počtu matice-vektorových součinů na jednu podúlohu ve srovnání s metodou trust-ncg. Navíc, řešení kvadratické podúlohy je přesnější než u metody trust-ncg.
Příklad s definicí Hessovy matice:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 19
Vyhodnocení funkce: 20
Vyhodnocení gradientu: 20
Vyhodnocení hessiana: 18
print(res.x)
[1. 1. 1. 1. 1.]
Příklad s funkcí produktu hessiana a libovolného vektoru:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Optimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 19
Vyhodnocení funkce: 20
Vyhodnocení gradientu: 20
Vyhodnocení hessiana: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algoritmus pro aproximaci v důvěrné oblasti
Všechny metody (Newton-CG, trust-ncg a trust-krylov) jsou velmi vhodné pro řešení rozsáhlých úloh (s tisíci proměnnými). To je způsobeno tím, že podkladový algoritmus konjugovaných gradientů předpokládá přibližné určení inverzní matice hessiana. Řešení se nachází iterativně, bez explicitní dekompozice hessiana. Jelikož je třeba určit pouze funkci pro součin hessiana a libovolného vektoru, je tento algoritmus zvlášť výhodný pro práci se sparse (páskovými diagonálními) maticemi. To zajišťuje nízké náklady na paměť a velkou úsporu času.
U středně velkých úlohách náklady na ukládání a faktorizaci Hessovy matice nejsou rozhodující. To znamená, že lze dosáhnout řešení za méně iterací tím, že se podúlohy důvěryhodnosti vyřeší téměř přesně. K tomu se některé nelineární rovnice iterativně řeší pro každou kvadratickou podúlohu. Takové řešení obvykle vyžaduje 3 nebo 4 rozklady Choleského Hessovy matice. Výsledkem je, že metoda konverguje za méně iterací a vyžaduje méně výpočtů cílové funkce než jiné realizované metody důvěryhodné oblasti. Tento algoritmus předpokládá pouze určení plné Hessovy matice a nepodporuje možnost použití funkce součinu hessianu a libovolného vektoru.
Příklad s minimalizací Rosenbrockovy funkce:
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimalizace byla úspěšně ukončena.
Aktuální hodnota funkce: 0.000000
Iterace: 13
Vyhodnocení funkce: 14
Vyhodnocení gradientu: 13
Vyhodnocení hessians: 14
array([1., 1., 1., 1., 1.])Tímto se asi zastavíme. V dalším článku se pokusím říci to nejzajímavější o podmínkové minimalizaci, aplikacích minimalizace při řešení úloh aproximace, minimalizaci funkce jedné proměnné, libovolných minimalizátorech a hledání kořenů soustavy rovnic pomocí balíku scipy.optimize.
Zdroj:
Zdroj: habr.com
