SciPy (izgovorjeno sai pie) je matematični aplikacijski paket, ki temelji na razširitvi Numpy Python. S SciPy vaša interaktivna seja Python postane isto popolno podatkovno znanost in okolje za izdelavo prototipov kompleksnega sistema kot MATLAB, IDL, Octave, R-Lab in SciLab. Danes želim na kratko spregovoriti o tem, kako uporabiti nekatere dobro znane optimizacijske algoritme v paketu scipy.optimize. Podrobnejšo in posodobljeno pomoč pri uporabi funkcij lahko vedno dobite z ukazom help() ali s kombinacijo Shift+Tab.
Predstavitev
Da bi sebi in bralcem prihranili iskanje in branje primarnih virov, bodo povezave do opisov metod v glavnem na Wikipediji. Ti podatki praviloma zadostujejo za razumevanje metod na splošno in pogojev za njihovo uporabo. Če želite razumeti bistvo matematičnih metod, sledite povezavam do bolj verodostojnih publikacij, ki jih najdete na koncu vsakega članka ali v vašem priljubljenem iskalniku.
Modul scipy.optimize torej vključuje izvedbo naslednjih postopkov:
- Pogojna in brezpogojna minimizacija skalarnih funkcij več spremenljivk (minim) z uporabo različnih algoritmov (Nelder-Mead simplex, BFGS, Newton konjugirani gradienti, и )
- Globalna optimizacija (na primer: , )
- Minimiziranje ostankov (least_squares) in algoritmi za prilagajanje krivulj z uporabo nelinearnih najmanjših kvadratov (curve_fit)
- Minimiziranje skalarnih funkcij ene spremenljivke (minim_scalar) in iskanje korenin (root_scalar)
- Večdimenzionalni reševalci sistema enačb (koren) z različnimi algoritmi (hibridni Powell, ali metode velikega obsega, kot je npr ).
V tem članku bomo obravnavali samo prvi element s tega celotnega seznama.
Brezpogojna minimizacija skalarne funkcije več spremenljivk
Funkcija minim iz paketa scipy.optimize nudi splošen vmesnik za reševanje problemov pogojne in brezpogojne minimizacije skalarnih funkcij več spremenljivk. Za prikaz, kako deluje, bomo potrebovali ustrezno funkcijo več spremenljivk, ki jih bomo minimizirali na različne načine.
Za te namene je Rosenbrockova funkcija N spremenljivk popolna, ki ima obliko:
Kljub temu, da so funkcija Rosenbrock in njeni Jacobijevi in Hessianovi matriki (prvi oziroma drugi odvod) že definirani v paketu scipy.optimize, jo bomo definirali 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)Zaradi jasnosti narišimo v 3D vrednosti Rosenbrockove funkcije dveh spremenljivk.
Risalna koda
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()
Vnaprej vedeti, da je minimum 0 pri , poglejmo primere, kako določiti najmanjšo vrednost Rosenbrockove funkcije z različnimi postopki scipy.optimize.
Simpleksna metoda Nelder-Mead
Naj obstaja začetna točka x0 v 5-dimenzionalnem prostoru. Z algoritmom poiščimo točko minimuma Rosenbrockove funkcije, ki ji je najbližja (algoritem je določen kot vrednost parametra metode):
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.]Simpleksna metoda je najpreprostejši način minimiziranja eksplicitno definirane in dokaj gladke funkcije. Ne zahteva izračuna odvodov funkcije, dovolj je, da navedete samo njene vrednosti. Metoda Nelder-Mead je dobra izbira za preproste probleme minimizacije. Ker pa ne uporablja ocen gradienta, lahko iskanje minimuma traja dlje.
Powellova metoda
Drug optimizacijski algoritem, v katerem se izračunajo samo vrednosti funkcij, je . Če ga želite uporabiti, morate v funkciji minim nastaviti method = '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.]Algoritem Broyden-Fletcher-Goldfarb-Shanno (BFGS).
Za hitrejšo konvergenco k rešitvi je postopek uporablja gradient ciljne funkcije. Gradient je mogoče določiti kot funkcijo ali izračunati z uporabo razlik prvega reda. V vsakem primeru metoda BFGS običajno zahteva manj klicev funkcij kot metoda simpleksa.
Poiščimo odvod Rosenbrockove funkcije v analitični obliki:
Ta izraz je veljaven za izpeljanke vseh spremenljivk razen prve in zadnje, ki sta definirani kot:
Poglejmo funkcijo Python, ki izračuna ta 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 derFunkcija izračuna gradienta je navedena kot vrednost parametra jac funkcije minim, kot je prikazano spodaj.
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]Algoritem konjugiranega gradienta (Newton)
Algoritem je modificirana Newtonova metoda.
Newtonova metoda temelji na aproksimaciji funkcije v lokalnem območju s polinomom druge stopnje:
če je matrika drugih odvodov (Hessian matrix, Hessian).
Če je Hessian pozitivno določen, potem lahko lokalni minimum te funkcije najdemo z enačenjem ničelnega gradienta kvadratne oblike z nič. Rezultat bo izraz:
Inverzni Hessian se izračuna z metodo konjugiranega gradienta. Spodaj je podan primer uporabe te metode za minimiziranje Rosenbrockove funkcije. Če želite uporabiti metodo Newton-CG, morate določiti funkcijo, ki izračuna Hessian.
Hessian Rosenbrockove funkcije v analitični obliki je enak:
če и , definirajte matriko .
Preostali neničelni elementi matrike so enaki:
Na primer, v petdimenzionalnem prostoru N = 5 ima Hessova matrika za Rosenbrockovo funkcijo obliko pasu:
Koda, ki izračuna ta Hessian skupaj s kodo za minimiziranje Rosenbrockove funkcije z uporabo metode konjugiranega gradienta (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]Primer z definicijo funkcije zmnožka hessenovega in poljubnega vektorja
Pri težavah v resničnem svetu lahko računanje in shranjevanje celotne Hessove matrike zahteva precej časa in pomnilniških virov. V tem primeru pravzaprav ni potrebe po določanju same Hessove matrike, ker postopek minimizacije zahteva samo vektor, ki je enak zmnožku Hessiana z drugim poljubnim vektorjem. Tako je z računskega vidika veliko bolje takoj definirati funkcijo, ki vrne rezultat zmnožka Hessana s poljubnim vektorjem.
Razmislite o funkciji Hess, ki vzame vektor minimizacije kot prvi argument in poljuben vektor kot drugi argument (skupaj z drugimi argumenti funkcije, ki jo je treba minimizirati). V našem primeru izračun produkta Hessiana Rosenbrockove funkcije s poljubnim vektorjem ni zelo težaven. če p je poljuben vektor, potem produkt izgleda kot:
Funkcija, ki izračuna zmnožek hessenovega in poljubnega vektorja, se kot vrednost argumenta hessp posreduje funkciji minimiziranja:
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: 66Algoritem regije zaupanja konjugiranega gradienta (Newton)
Slabo pogojevanje Hessove matrike in nepravilne smeri iskanja lahko povzročijo, da je Newtonov konjugirani gradientni algoritem neučinkovit. V takih primerih se daje prednost (območje zaupanja) konjugirani Newtonovi gradienti.
Primer z definicijo Hessove matrike:
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.]Primer z zmnožkom Hessove funkcije in poljubnega vektorja:
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.]Metode tipa Krylov
Tako kot metoda zaupanja ncg so metode tipa Krylov zelo primerne za reševanje problemov velikega obsega, ker uporabljajo samo produkte matričnih vektorjev. Njihovo bistvo je rešiti problem v območju zaupanja, omejenem z okrnjenim Krylovovim podprostorom. Za negotove probleme je bolje uporabiti to metodo, saj uporablja manjše število nelinearnih iteracij zaradi manjšega števila matrično-vektorskih produktov na podproblem v primerjavi z metodo trust-ncg. Poleg tega je rešitev kvadratnega podproblema najdena natančneje kot z uporabo metode trust-ncg.
Primer z definicijo Hessove matrike:
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.]
Primer z zmnožkom Hessove funkcije in poljubnega vektorja:
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.]
Algoritem za približno rešitev v območju zaupanja
Vse metode (Newton-CG, trust-ncg in trust-krylov) so zelo primerne za reševanje velikih problemov (s tisočimi spremenljivkami). To je posledica dejstva, da osnovni algoritem konjugiranega gradienta implicira približno določitev inverzne Hessove matrike. Rešitev se najde iterativno, brez eksplicitne razširitve Hessiana. Ker morate definirati samo funkcijo za zmnožek Hessovega in poljubnega vektorja, je ta algoritem še posebej dober za delo z redkimi (pasovnimi diagonalnimi) matricami. To zagotavlja nizke stroške pomnilnika in znatne prihranke časa.
Za srednje velike težave stroški shranjevanja in faktoriziranja Hessian niso kritični. To pomeni, da je možno dobiti rešitev v manj ponovitvah, s čimer skoraj natančno rešimo podprobleme regije zaupanja. Da bi to naredili, se nekatere nelinearne enačbe rešujejo iterativno za vsak kvadratni podproblem. Takšna rešitev običajno zahteva 3 ali 4 Choleskyjeve razgradnje Hessove matrike. Posledično metoda konvergira v manj ponovitvah in zahteva manj izračunov ciljne funkcije kot druge implementirane metode regije zaupanja. Ta algoritem vključuje le določanje celotne Hessove matrike in ne podpira zmožnosti uporabe produktne funkcije Hessovega in poljubnega vektorja.
Primer z minimizacijo Rosenbrockove funkcije:
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.])Verjetno se bomo tam ustavili. V naslednjem članku bom poskušal povedati najbolj zanimive stvari o pogojni minimizaciji, uporabi minimizacije pri reševanju aproksimacijskih problemov, minimizaciji funkcije ene spremenljivke, poljubnih minimizatorjih in iskanju korenin sistema enačb z uporabo scipy.optimize paket.
Vir:
Vir: www.habr.com