SciPy (vyslovováno sai pie) je matematický aplikační balíček založený na rozšíření Numpy Python. Se SciPy se vaše interaktivní relace v Pythonu stane stejným prostředím pro kompletní datovou vědu a komplexním systémovým prototypováním jako MATLAB, IDL, Octave, R-Lab a SciLab. Dnes chci krátce pohovořit o tom, jak používat některé známé optimalizační algoritmy v balíčku scipy.optimize. Podrobnější a aktuálnější nápovědu k používání funkcí lze vždy získat pomocí příkazu help() nebo pomocí Shift+Tab.
úvod
Abyste sebe i čtenáře ušetřili hledání a čtení primárních zdrojů, budou odkazy na popisy metod hlavně na Wikipedii. Tyto informace zpravidla postačují k pochopení metod obecně a podmínek jejich aplikace. Abyste pochopili podstatu matematických metod, postupujte podle odkazů na autoritativnější publikace, které najdete na konci každého článku nebo ve svém oblíbeném vyhledávači.
Modul scipy.optimize tedy zahrnuje implementaci následujících procedur:
- Podmíněná a nepodmíněná minimalizace skalárních funkcí několika proměnných (minim) pomocí různých algoritmů (Nelder-Mead simplex, BFGS, Newtonovy konjugované gradienty, и )
- Globální optimalizace (například: , )
- Minimalizace zbytků (nejmenší_čtverce) a algoritmy prokládání křivek pomocí nelineárních nejmenších čtverců (curve_fit)
- Minimalizace skalárních funkcí jedné proměnné (minim_scalar) a hledání kořenů (root_scalar)
- Vícerozměrné řešiče soustav rovnic (kořen) využívající různé algoritmy (hybridní Powell, nebo metody ve velkém měřítku, jako např ).
V tomto článku se budeme zabývat pouze první položkou z celého tohoto seznamu.
Bezpodmínečná minimalizace skalární funkce více proměnných
Funkce minim z balíku scipy.optimize poskytuje obecné rozhraní pro řešení problémů podmíněné a nepodmíněné minimalizace skalárních funkcí několika proměnných. Abychom demonstrovali, jak to funguje, budeme potřebovat vhodnou funkci několika proměnných, které budeme různými způsoby minimalizovat.
Pro tyto účely je perfektní Rosenbrockova funkce N proměnných, která má tvar:
Navzdory tomu, že Rosenbrockova funkce a její Jacobiho a Hessova matice (první a druhá derivace) jsou již definovány v balíčku scipy.optimize, nadefinujeme si ji 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)Pro přehlednost nakreslíme ve 3D hodnoty Rosenbrockovy funkce dvou proměnných.
Kód výkresu
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()
Předem vědět, že minimum je 0 at 
, podívejme se na příklady, jak určit minimální hodnotu Rosenbrockovy funkce pomocí různých postupů scipy.optimize.
Nelder-Mead simplexová metoda
Nechť je počáteční bod x0 v 5-rozměrném prostoru. Pomocí algoritmu najdeme minimální bod Rosenbrockovy funkce, který je k ní nejblíže (algoritmus je uveden jako hodnota parametru metody):
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á metoda 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 Nelder-Mead je dobrou volbou pro jednoduché problémy s minimalizací. Protože však nepoužívá odhady gradientu, může nalezení minima trvat déle.
Powellova metoda
Dalším optimalizačním algoritmem, ve kterém se počítají pouze funkční hodnoty, je . Chcete-li ji použít, musíte ve funkci minim nastavit metodu = '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).
Chcete-li získat rychlejší konvergenci k řešení, postup používá gradient účelové funkce. Gradient lze zadat jako funkci nebo vypočítat pomocí rozdílů prvního řádu. V každém případě metoda BFGS obvykle vyžaduje méně volání funkcí než simplexní metoda.
Najdeme derivaci Rosenbrockovy funkce v analytické podobě:
Tento výraz je platný pro derivace všech proměnných kromě první a poslední, které jsou definovány jako:
Podívejme se na funkci Pythonu, která vypočítá 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 derFunkce výpočtu přechodu je určena jako hodnota parametru jac funkce minima, jak je uvedeno níže.
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 metoda.
Newtonova metoda je založena na aproximaci funkce v lokální oblasti polynomem druhého stupně:
kde 
je matice druhých derivací (Hessova matice, Hessian).
Je-li Hessian kladně definitní, pak lze lokální minimum této funkce nalézt přirovnáním nulového gradientu kvadratické formy k nule. Výsledkem bude výraz:
Inverzní Hessian se vypočítá pomocí metody konjugovaného gradientu. Níže je uveden příklad použití této metody k minimalizaci funkce Rosenbrock. Chcete-li použít metodu Newton-CG, musíte zadat funkci, která vypočítá Hessian.
Hessián Rosenbrockovy funkce v analytické formě se rovná:
kde 
и 
, definujte matici 
.
Zbývající nenulové prvky matice se rovnají:
Například v pětirozměrném prostoru N = 5 má Hessova matice pro Rosenbrockovu funkci tvar pásma:
Kód, který vypočítá tento Hessian spolu s kódem pro minimalizaci Rosenbrockovy funkce pomocí metody 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]Příklad s definicí součinové funkce Hessianu a libovolného vektoru
V problémech reálného světa může výpočet a ukládání celé Hessovy matice vyžadovat značné časové a paměťové zdroje. V tomto případě vlastně není potřeba specifikovat samotnou Hessovu matici, protože minimalizační procedura vyžaduje pouze vektor rovný součinu Hessianu s jiným libovolným vektorem. Z výpočetního hlediska je tedy mnohem vhodnější okamžitě definovat funkci, která vrací výsledek součinu Hessianu s libovolným vektorem.
Uvažujme funkci hess, která bere minimalizační vektor jako první argument, a libovolný vektor jako druhý argument (spolu s dalšími argumenty funkce, která má být minimalizována). V našem případě není výpočet součinu Hessianovy Rosenbrockovy funkce s libovolným vektorem příliš obtížný. Li p je libovolný vektor, pak součin 
má tvar:
Funkce, která vypočítá součin Hessianu a libovolného vektoru, je předána jako hodnota argumentu hessp funkci minimalizace:
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 konjugované oblasti důvěryhodnosti gradientu (Newton)
Špatná úprava Hessovy matice a nesprávné směry hledání mohou způsobit, že Newtonův algoritmus sdruženého gradientu bude neúčinný. V takových případech se dává přednost (trust-region) konjugovat Newtonovy gradienty.
Příklad s definicí hessovské 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.]Příklad s funkcí součinu Hessianu 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 0
[1. 1. 1. 1. 1.]Metody Krylovova typu
Stejně jako metoda trust-ncg jsou metody Krylovova typu dobře vhodné pro řešení rozsáhlých problémů, protože používají pouze produkty maticových vektorů. Jejich podstatou je vyřešit problém v oblasti důvěry omezené zkráceným Krylovovým podprostorem. Pro nejisté problémy je lepší použít tuto metodu, protože používá menší počet nelineárních iterací kvůli menšímu počtu matic-vektorových produktů na podproblém ve srovnání s metodou trust-ncg. Řešení kvadratického dílčího problému je navíc nalezeno přesněji než pomocí metody trust-ncg.
Příklad s definicí hessovské 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.]
Příklad s funkcí součinu Hessianu a libovolného vektoru:
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 pro přibližné řešení v oblasti spolehlivosti
Všechny metody (Newton-CG, trust-ncg a trust-krylov) jsou vhodné pro řešení rozsáhlých problémů (s tisíci proměnnými). To je způsobeno skutečností, že základní algoritmus konjugovaného gradientu implikuje přibližné určení inverzní Hessovy matice. Řešení se nalézá iterativně, bez explicitního rozšíření hessia. Vzhledem k tomu, že potřebujete definovat funkci pouze pro součin Hessianu a libovolného vektoru, je tento algoritmus obzvláště vhodný pro práci s řídkými (pásmovými diagonálními) maticemi. To poskytuje nízké náklady na paměť a významnou úsporu času.
U středně velkých problémů nejsou náklady na ukládání a faktoring Hessian kritické. To znamená, že je možné získat řešení v menším počtu iterací, které řeší dílčí problémy důvěryhodné oblasti téměř přesně. K tomu jsou některé nelineární rovnice řešeny iterativně pro každý kvadratický podproblém. Takové řešení obvykle vyžaduje 3 nebo 4 Choleského rozklady Hessovy matice. Výsledkem je, že metoda konverguje v menším počtu iterací a vyžaduje méně výpočtů objektivních funkcí než jiné implementované metody oblasti spolehlivosti. Tento algoritmus zahrnuje pouze určení kompletní hessovské matice a nepodporuje možnost použití součinové funkce hessovské matice a libovolného vektoru.
Příklad s minimalizací funkce Rosenbrock:
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.])Tam se asi zastavíme. V příštím článku se pokusím říci to nejzajímavější o podmíněné minimalizaci, aplikaci minimalizace při řešení aproximačních úloh, minimalizaci funkce jedné proměnné, libovolných minimalizátorech a hledání kořenů soustavy rovnic pomocí scipy.optimize balík.
Zdroj:
Zdroj: www.habr.com
