SciPy (izgovara se kao saj paj) je paket primenjenih matematičkih procedura, zasnovan na ekstenziji Numpy Python. Sa SciPy, interaktivna Python sesija postaje potpuno funkcionalno okruženje za obradu podataka i prototipiranje složenih sistema, kao što su MATLAB, IDL, Octave, R-Lab i SciLab. Danas ću ukratko govoriti o tome kako primeniti neka poznata algoritma optimizacije u paketu scipy.optimize. Detaljniju i aktuelnu pomoć u vezi sa korišćenjem funkcija uvek možete dobiti pomoću komande help() ili pritiskom na Shift+Tab.
Uvod
Da bih sebi i čitaocima olakšao pretragu i čitanje izvora, linkovi ka opisima metoda biće uglavnom na Vikipediji. Obično je ova informacija dovoljna za opšte razumevanje metoda i uslova njihovog korišćenja. Da bismo razumeli suštinu matematičkih metoda, upućujemo vas na autoritativnije publikacije koje možete pronaći na kraju svake stranice ili u omiljenoj pretraživačkoj mašini.
Dakle, modul scipy.optimize uključuje implementaciju sledećih procedura:
- Uslovne i bezuslovne minimizacije skalarnih funkcija više varijabli (minim) koristeći razne algoritme (Nelder-Mead simplex, BFGS, Newtonovi povezani gradijenti, in )
- Globalne optimizacije (na primer: , )
- Minimizaciju ostataka (least_squares) i algoritme za prilagođavanje krivulja nelinearnim MNK (curve_fit)
- Minimizaciju skalarnih funkcija jedne varijable (minim_scalar) i pretragu korena (root_scalar)
- Višedimenzionalni rešavači sistema jednačina (root) koristeći razne algoritme (hipotetički Powell, ili velike metode kao što su ).
U ovom članku razmotrićemo samo prvu stavku sa ovog spiska.
Bezuslovna minimizacija skalarne funkcije više varijabli
Funkcija minim iz paketa scipy.optimize pruža opšti interfejs za rešavanje problema uslovne i bezuslovne minimizacije skalarnih funkcija više varijabli. Da bismo demonstrirali njeno funkcionisanje, potrebna nam je odgovarajuća funkcija više varijabli koju ćemo različito minimizovati.
Za te svrhe savršeno je pogodna Rosenbrockova funkcija od N varijabli koja ima oblik:
Čeprav je funkcija Rosenbrock ter njene Jakobijeve in Hessove matrike (prve in druge derivacije) že definirane v paketu scipy.optimize, jo bomo določili sami.
import numpy as np
def rosen(x):
"""Funkcija Rosenbrock"""
return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)Za boljšo predstavo bomo v 3D narisali vrednosti funkcije Rosenbrock glede na dve spremenljivki.
Koda za risanje
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# Nastavimo 3D graf
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Nastavimo kot gledanja
ax.view_init(45, 30)
# Ustvarimo podatke za 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]))
# Narišemo površino
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
Zavedamo se, da je minimum enak 0 pri 
, poglejmo primere, kako določiti najmanjšo vrednost funkcije Rosenbrock s pomočjo različnih postopkov scipy.optimize.
Metoda simplex Nelder-Mead (Nelder-Mead)
Naj gre za začetno točko x0 v 5-dimenzionalnem prostoru. Poiščimo najbližjo točko minimuma funkcije Rosenbrock s pomočjo algoritma (algoritem naveden kot vrednost 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)Optimizacija je bila uspešno zaključena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 339
Ocene funkcij: 571
[1. 1. 1. 1. 1.]Metoda simplex je najpreprostejši način, kako minimizirati jasno definirano in precej gladko funkcijo. Ne zahteva izračunavanja odvoda funkcije, dovolj je, da navedemo njene vrednosti. Metoda Nelder-Mead je dobra izbira za enostavne naloge minimizacije. Vendar pa, ker ne uporablja ocen gradienta, lahko iščenje minimuma traja dlje.
Metoda Powell
Drug algoritem optimizacije, ki uporablja le vrednosti funkcij, je Za njeno uporabo je potrebno nastaviti method = ‘powell’ v funkciji 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)Optimizacija je bila uspešno zaključena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 19
Ocene funkcij: 1622
[1. 1. 1. 1. 1.]Algoritem Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Za hitrejšo konvergenco k rešitvi, postopek uporablja gradient ciljne funkcije. Gradient je lahko dan v obliki funkcije ali izračunan s pomočjo prvih razlik. V vsakem primeru običajno metoda BFGS zahteva manj klicev funkcij kot simplex metoda.
Izračunajmo analitično derivacijo Rosenbrockove funkcije:
Ta izraz velja za derivacije vseh spremenljivk, razen prve in zadnje, ki sta določeni 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 za izračun gradienta se navede 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)Optimizacija je bila uspešno končana.
Trenutna vrednost funkcije: 0.000000
Iteracije: 25
Klici funkcij: 30
Klici gradienta: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Algoritem konjugiranih gradientov (Newton)
Algoritem je modificirana Newtonova metoda.
Newtonova metoda temelji na aproksimaciji funkcije v lokalnem območju s polinomom druge stopnje:
kjer 
je matrika drugih derivatov (Hessejeva matrica, hessian).
Če je hessian pozitivno določen, lahko lokalni minimum te funkcije najdemo tako, da ničelni gradient kvadratne forme enačimo na nič. Rezultat bo izraz:
Obratni hessian se izračuna s pomočjo metode konjugiranih gradientov. Spodaj je prikazan primer uporabe te metode za minimizacijo Rosenbrockove funkcije. Za uporabo metode Newton-CG je treba navesti funkcijo, ki izračuna hessian.
Hessian Rosenbrockove funkcije v analitični obliki je enak:
kjer 
in 
, določajo matriko 
.
Preostali nenulni elementi matrike so enaki:
Na primer, v petdimenzionalnem prostoru N = 5, Hessejeva matrica za Rosenbrockovo funkcijo ima pasovno obliko:
Koda, ki izračuna ta hessian skupaj s kodo za minimizacijo Rosenbrockove funkcije z metodo konjugiranih gradientov (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)Optimizacija je uspešno zaključena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 24
Evaluacije funkcije: 33
Evaluacije gradijenta: 56
Evaluacije hessian-a: 24
[1. 1. 1. 0.99999999 0.99999999]Primer sa definicijom funkcije proizvoda Hessiana i proizvoljnog vektora
U stvarnim zadacima, izračunavanje i čuvanje cele Hessian matrice može zahtevati značajne resurse vremena i memorije. Pri tome, zapravo nije potrebno zadavati samu Hessian matricu, jer za proces minimizacije treba samo vektor koji je jednak proizvodu Hessiana sa drugim proizvoljnim vektorom. Na taj način, sa računarskog stanovišta, daleko je povoljnije odmah definisati funkciju koja vraća rezultat proizvoda Hessiana sa proizvoljnim vektorom.
Razmotrićemo funkciju hess, koja prima vektor minimizacije kao prvi argument, a proizvoljni vektor kao drugi argument (uz druge argumente minimizovane funkcije). U našem slučaju, izračunavanje proizvoda Hessiana funkcije Rosenbrock sa proizvoljnim vektorom nije previše složeno. Ako p — je proizvoljni vektor, onda proizvod 
ima oblik:
Funkcija koja izračunava proizvod Hessiana i proizvoljnog vektora, prosledjuje se kao vrednost argumenta hessp funkciji 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})
Optimizacija je uspešno zaključena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 24
Evaluacije funkcije: 33
Evaluacije gradijenta: 56
Evaluacije hessian-a: 66Algoritam poverljive oblasti (trust region) spojene gradijente (Newton)
Loša određeničnost Hessian matrice i netačni pravci pretrage mogu dovesti do toga da algoritam spojene gradijente Newton postane neučinkovit. U takvim slučajevima, prednost se daje (trust-region) spojene gradijente Newton-u.
Primer sa definicijom Hessian matrice:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimizacija je uspešno završena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 20
Evaluacije funkcije: 21
Evaluacije gradijenta: 20
Evaluacije Hese: 19
[1. 1. 1. 1. 1.]Primer sa funkcijom proizvoda Hese i proizvoljnog vektora:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimizacija je uspešno završena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 20
Evaluacije funkcije: 21
Evaluacije gradijenta: 20
Evaluacije Hese: 0
[1. 1. 1. 1. 1.]Metode tipa Krylov
Slično metodi trust-ncg, metode tipa Krylov su dobro prilagođene za rešavanje velikih problema, jer koriste samo matrično-vektorske proizvode. Njihova suština leži u rešavanju problema unutar poverljive oblasti, ograničene skraćenim Krylovim podprostorom. Za nedeterminističke probleme bolje je koristiti ovu metodu, jer zahteva manje nelinearnih iteracija zahvaljujući manjem broju matrično-vektorskih proizvoda po podproblem, u poređenju sa metodom trust-ncg. Pored toga, rešenje kvadratnog podproblema je preciznije nego kod metode trust-ncg.
Primer sa definicijom Hessian matrice:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Optimizacija je uspešno završena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 19
Evaluacije funkcije: 20
Evaluacije gradijenta: 20
Evaluacije Hese: 18
print(res.x)
[1. 1. 1. 1. 1.]
Primer sa funkcijom proizvoda Hese i proizvoljnog vektora:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Optimizacija je uspešno završena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 19
Evaluacije funkcije: 20
Evaluacije gradijenta: 20
Evaluacije Hese: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algoritam približnog rešenja u poverljivoj oblasti
Sve metode (Newton-CG, trust-ncg i trust-krylov) su veoma pogodna za rešavanje velikih problema (sa hiljadama varijabli). To je zato što osnovni algoritam konjugovanih gradijenata podrazumeva približno nalaženje inverzne Heseove matrice. Rešenje se pronalazi iterativno, bez eksplicitnog razlaganja Heseove matrice. Pošto je potrebno definisati samo funkciju za proizvod Heseove matrice i proizvoljnog vektora, ovaj algoritam je posebno dobar za rad sa retkim (traka dijagonalnim) matricama. To obezbeđuje niske troškove memorije i značajnu uštedu vremena.
V nalogah srednje velikosti stroški shranjevanja in faktoriranja Hessiane nimajo odločilnega pomena. To pomeni, da je mogoče doseči rešitev z manjšim številom iteracij, tako da se podnaloge zaupne regije rešijo skoraj natančno. Za to se nekatera nelinearna enačba iterativno rešujejo za vsak kvadratni podproblem. Ta rešitev običajno zahteva 3 ali 4 razcepitve Choleskega Hessiane matrike. Kot rezultat, metoda konvergira z manjšim številom iteracij in zahteva manj izračunov ciljne funkcije kot druge izvedene metode zaupanja. Ta algoritem predvideva le določitev popolne Hessiane matrike in ne podpira možnosti uporabe funkcije produkta Hessiane 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.xOptimizacija je uspešno zaključena.
Trenutna vrednost funkcije: 0.000000
Iteracije: 13
Ocene funkcij: 14
Ocene gradienta: 13
Ocene Hessiane: 14
array([1., 1., 1., 1., 1.])Na tem se, mislim, zaustavimo. V naslednjem članku bom poskusil predstaviti najbolj zanimivo o pogojni minimizaciji, aplikaciji minimizacije pri reševanju problemov aproksimacije, minimizaciji funkcije ene spremenljivke, poljubnih minimizatorjih in iskanju korenov sistema enačb s pomočjo paketa scipy.optimize.
Izvor:
Izvor: habr.com
