SciPy (izgovara se sai pie) je paket matematičkih aplikacija baziran na proširenju Numpy Python. Uz SciPy, vaša interaktivna Python sesija postaje ista potpuna nauka o podacima i okruženje za izradu prototipa kompleksnog sistema kao MATLAB, IDL, Octave, R-Lab i SciLab. Danas želim ukratko govoriti o tome kako koristiti neke dobro poznate algoritme optimizacije u paketu scipy.optimize. Detaljnija i ažurirana pomoć o korištenju funkcija uvijek se može dobiti pomoću naredbe help() ili pomoću Shift+Tab.
Uvod
Kako biste sebe i čitatelje spasili od pretraživanja i čitanja primarnih izvora, linkovi na opise metoda bit će uglavnom na Wikipediji. U pravilu, ove informacije su dovoljne za razumijevanje metoda općenito i uslova za njihovu primjenu. Da biste razumjeli suštinu matematičkih metoda, slijedite veze do mjerodavnijih publikacija, koje se nalaze na kraju svakog članka ili u vašoj omiljenoj tražilici.
Dakle, modul scipy.optimize uključuje implementaciju sljedećih procedura:
- Uslovna i bezuslovna minimizacija skalarnih funkcija nekoliko varijabli (minimum) korišćenjem različitih algoritama (Nelder-Mead simpleks, BFGS, Newton konjugirani gradijenti, и )
- Globalna optimizacija (na primjer: , )
- Minimiziranje ostataka (najmanji_kvadrati) i algoritmi prilagođavanja krivulje korištenjem nelinearnih najmanjih kvadrata (curve_fit)
- Minimiziranje skalarnih funkcija jedne varijable (minim_skalar) i traženje korijena (root_scalar)
- Multidimenzionalni rešavači sistema jednačina (koren) korišćenjem različitih algoritama (hibridni Powell, ili metode velikih razmjera kao npr ).
U ovom članku ćemo razmotriti samo prvu stavku sa cijele liste.
Bezuslovna minimizacija skalarne funkcije nekoliko varijabli
Minimalna funkcija iz paketa scipy.optimize pruža opšti interfejs za rešavanje problema uslovne i bezuslovne minimizacije skalarnih funkcija nekoliko varijabli. Da bismo demonstrirali kako to funkcionira, trebat će nam odgovarajuća funkcija nekoliko varijabli, koje ćemo minimizirati na različite načine.
Za ove svrhe, Rosenbrockova funkcija od N varijabli je savršena, koja ima oblik:
Unatoč činjenici da su Rosenbrockova funkcija i njene Jacobijeve i Hessijeve matrice (prvi i drugi derivati, respektivno) već definirane u paketu scipy.optimize, mi ćemo je sami definirati.
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)Radi jasnoće, nacrtajmo u 3D vrijednosti Rosenbrockove funkcije dvije varijable.
Kod za crtanje
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()
Znajući unaprijed da je minimum 0 at , pogledajmo primjere kako odrediti minimalnu vrijednost funkcije Rosenbrock koristeći različite scipy.optimize procedure.
Nelder-Mead simpleks metoda
Neka postoji početna tačka x0 u 5-dimenzionalnom prostoru. Nađimo minimalnu tačku Rosenbrockove funkcije koja joj je najbliža pomoću algoritma (algoritam je specificiran kao vrijednost 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.]Simpleks metoda je najjednostavniji način da se minimizira eksplicitno definirana i prilično glatka funkcija. Ne zahtijeva izračunavanje izvoda funkcije, dovoljno je navesti samo njene vrijednosti. Nelder-Meadova metoda je dobar izbor za jednostavne probleme minimizacije. Međutim, budući da ne koristi procjene gradijenta, može potrajati duže da se pronađe minimum.
Powellova metoda
Drugi algoritam optimizacije u kojem se izračunavaju samo vrijednosti funkcije je . Da biste ga koristili, trebate postaviti method = 'powell' u minimalnoj funkciji.
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.]Broyden-Fletcher-Goldfarb-Shanno (BFGS) algoritam
Da bi se postigla brža konvergencija do rješenja, postupak koristi gradijent funkcije cilja. Gradijent se može specificirati kao funkcija ili izračunati korištenjem razlika prvog reda. U svakom slučaju, BFGS metoda obično zahtijeva manje poziva funkcije nego simpleks metoda.
Nađimo izvod Rosenbrockove funkcije u analitičkom obliku:
Ovaj izraz vrijedi za derivate svih varijabli osim prve i posljednje, koje su definirane kao:
Pogledajmo Python funkciju koja izračunava ovaj gradijent:
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 proračuna gradijenta specificirana je kao vrijednost parametra jac minimalne funkcije, kao što je prikazano ispod.
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]Algoritam konjugiranog gradijenta (njutn)
Algoritam je modificirana Newtonova metoda.
Newtonova metoda temelji se na aproksimaciji funkcije u lokalnom području polinomom drugog stepena:
gdje je matrica drugih izvoda (Hessian matrica, Hessian).
Ako je Hessian pozitivno određen, tada se lokalni minimum ove funkcije može naći izjednačavanjem nulti gradijenta kvadratnog oblika sa nulom. Rezultat će biti izraz:
Inverzni Hessian se izračunava metodom konjugovanog gradijenta. Primjer korištenja ove metode za minimiziranje Rosenbrockove funkcije je dat u nastavku. Da biste koristili Newton-CG metodu, morate navesti funkciju koja izračunava Hessian.
Hessian Rosenbrockove funkcije u analitičkom obliku jednak je:
gdje и , definirati matricu .
Preostali elementi matrice koji nisu nula jednaki su:
Na primjer, u petodimenzionalnom prostoru N = 5, Hessiova matrica za Rosenbrockovu funkciju ima oblik trake:
Kod koji izračunava ovaj Hessian zajedno s kodom za minimiziranje Rosenbrockove funkcije korištenjem metode konjugiranog gradijenta (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]Primjer sa definicijom funkcije proizvoda Hessiana i proizvoljnog vektora
U problemima iz stvarnog svijeta, računanje i pohranjivanje cijele Hessian matrice može zahtijevati značajno vrijeme i memorijske resurse. U ovom slučaju zapravo nema potrebe specificirati samu Hessian matricu, jer procedura minimizacije zahtijeva samo vektor jednak proizvodu Hessiana sa drugim proizvoljnim vektorom. Dakle, sa računske tačke gledišta, mnogo je bolje odmah definisati funkciju koja vraća rezultat Hessianovog proizvoda sa proizvoljnim vektorom.
Uzmite u obzir hess funkciju, koja uzima vektor minimizacije kao prvi argument, a proizvoljni vektor kao drugi argument (zajedno sa drugim argumentima funkcije koju treba minimizirati). U našem slučaju, izračunavanje proizvoda Hessiana Rosenbrockove funkcije sa proizvoljnim vektorom nije teško. Ako p je proizvoljan vektor, onda je proizvod izgleda kao:
Funkcija koja izračunava proizvod Hessiana i proizvoljnog vektora prosljeđuje se kao vrijednost argumenta hessp 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: 66Algoritam regije povjerenja konjugiranog gradijenta (Newton)
Loše kondicioniranje Hessian matrice i netačni smjerovi traženja mogu uzrokovati da Njutnov algoritam konjugovanog gradijenta bude neefikasan. U takvim slučajevima prednost se daje (regija povjerenja) konjugiraju Newton gradijente.
Primjer 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 19
[1. 1. 1. 1. 1.]Primjer s funkcijom proizvoda Hessiana i proizvoljnim vektorom:
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 Krilov
Kao i metoda trust-ncg, metode tipa Krilov su vrlo pogodne za rješavanje problema velikih razmjera jer koriste samo proizvode matričnog vektora. Njihova suština je rješavanje problema u području povjerenja ograničenom skraćenim Krilovskim podprostorom. Za nesigurne probleme, bolje je koristiti ovu metodu, jer koristi manji broj nelinearnih iteracija zbog manjeg broja matrično-vektorskih proizvoda po podproblemu, u poređenju sa metodom trust-ncg. Osim toga, rješenje kvadratnog podproblema se pronalazi preciznije nego korištenjem trust-ncg metode.
Primjer sa definicijom Hessian matrice:
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.]
Primjer s funkcijom proizvoda Hessiana i proizvoljnim vektorom:
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.]
Algoritam za približno rješenje u području povjerenja
Sve metode (Newton-CG, trust-ncg i trust-krylov) su pogodne za rješavanje problema velikih razmjera (sa hiljadama varijabli). To je zbog činjenice da osnovni algoritam konjugiranog gradijenta podrazumijeva približno određivanje inverzne Hesove matrice. Rješenje se nalazi iterativno, bez eksplicitnog proširenja Hessiana. Budući da trebate samo definirati funkciju za proizvod Hessiana i proizvoljnog vektora, ovaj algoritam je posebno dobar za rad sa rijetkim (dijagonalnim) matricama. Ovo obezbeđuje niske troškove memorije i značajnu uštedu vremena.
Za probleme srednje veličine, troškovi skladištenja i faktoringa Hessiana nisu kritični. To znači da je moguće dobiti rješenje u manje iteracija, rješavajući podprobleme regije povjerenja gotovo tačno. Da bi se to postiglo, neke nelinearne jednadžbe se rješavaju iterativno za svaki kvadratni podproblem. Takvo rješenje obično zahtijeva 3 ili 4 Choleskyjeve dekompozicije Hesove matrice. Kao rezultat toga, metoda konvergira u manje iteracija i zahtijeva manje proračuna funkcije cilja nego druge implementirane metode regiona povjerenja. Ovaj algoritam uključuje samo određivanje kompletne Hessian matrice i ne podržava mogućnost korištenja funkcije proizvoda Hessiana i proizvoljnog vektora.
Primjer sa minimiziranjem 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.])Verovatno ćemo tu stati. U sljedećem članku pokušat ću ispričati najzanimljivije stvari o uslovnoj minimizaciji, primjeni minimizacije u rješavanju aproksimacijskih problema, minimiziranju funkcije jedne varijable, proizvoljnim minimizatorima i pronalaženju korijena sistema jednadžbi pomoću scipy.optimize paket.
izvor:
izvor: www.habr.com