SciPy (izgovara se sai pie) je matematički aplikacijski paket temeljen na ekstenziji Numpy Python. Uz SciPy, vaša interaktivna Python sesija postaje ista potpuna znanost o podacima i okruženje za izradu prototipa složenog sustava kao MATLAB, IDL, Octave, R-Lab i SciLab. Danas želim ukratko govoriti o tome kako koristiti neke dobro poznate optimizacijske algoritme u paketu scipy.optimize. Detaljniju i ažurniju pomoć o korištenju funkcija uvijek možete dobiti pomoću naredbe help() ili pomoću Shift+Tab.
Uvod
Kako biste sebe i čitatelje poštedjeli traženja i čitanja primarnih izvora, poveznice na opise metoda bit će uglavnom na Wikipediji. Ti su podaci u pravilu dovoljni za razumijevanje metoda u općim crtama i uvjeta njihove primjene. Da biste razumjeli bit matematičkih metoda, slijedite poveznice na autoritativnije publikacije koje možete pronaći na kraju svakog članka ili u svojoj omiljenoj tražilici.
Dakle, modul scipy.optimize uključuje implementaciju sljedećih postupaka:
- Uvjetna i bezuvjetna minimizacija skalarnih funkcija više varijabli (minim) pomoću različitih algoritama (Nelder-Mead simplex, BFGS, Newton konjugirani gradijenti, и )
- Globalna optimizacija (na primjer: , )
- Minimiziranje ostataka (least_squares) i algoritmi za prilagođavanje krivulje koji koriste nelinearne najmanje kvadrate (curve_fit)
- Minimiziranje skalarnih funkcija jedne varijable (minim_scalar) i traženje korijena (root_scalar)
- Višedimenzionalni rješavači sustava jednadžbi (korijen) pomoću različitih algoritama (hibridni Powell, ili metode velikih razmjera kao što su ).
U ovom ćemo članku razmotriti samo prvu stavku s cijelog popisa.
Bezuvjetna minimizacija skalarne funkcije više varijabli
Funkcija minim iz paketa scipy.optimize pruža opće sučelje za rješavanje problema uvjetne i bezuvjetne minimizacije skalarnih funkcija nekoliko varijabli. Da bismo demonstrirali kako to radi, trebat će nam odgovarajuća funkcija nekoliko varijabli, koje ćemo minimizirati na različite načine.
Za ove svrhe savršena je Rosenbrockova funkcija N varijabli koja ima oblik:
Unatoč činjenici da su Rosenbrockova funkcija i njezine Jacobijeve i Hessianove matrice (prva odnosno druga derivacija) već definirane u paketu scipy.optimize, definirat ćemo je 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)Radi jasnoće, nacrtajmo u 3D vrijednosti Rosenbrockove funkcije dviju varijabli.
Kod crtanja
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 Rosenbrockove funkcije pomoću različitih postupaka scipy.optimize.
Nelder-Meadova simpleks metoda
Neka postoji početna točka x0 u 5-dimenzionalnom prostoru. Pronađimo točku minimuma Rosenbrockove funkcije koja joj je najbliža pomoću algoritma (algoritam je naveden 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.]Simpleksna metoda je najjednostavniji način minimiziranja eksplicitno definirane i prilično glatke funkcije. Ne zahtijeva izračunavanje derivacija funkcije, dovoljno je navesti samo njezine vrijednosti. Nelder-Mead metoda je dobar izbor za jednostavne probleme minimizacije. Međutim, budući da ne koristi procjene gradijenata, može biti potrebno više vremena da se pronađe minimum.
Powellova metoda
Drugi optimizacijski algoritam u kojem se izračunavaju samo vrijednosti funkcije je . Da biste ga koristili, trebate postaviti method = 'powell' u 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)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 prema rješenju, postupak koristi gradijent funkcije cilja. Gradijent se može navesti kao funkcija ili izračunati pomoću razlika prvog reda. U svakom slučaju, BFGS metoda obično zahtijeva manje poziva funkcija nego simpleks metoda.
Nađimo izvod Rosenbrockove funkcije u analitičkom obliku:
Ovaj izraz vrijedi za derivacije svih varijabli osim prve i zadnje, 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 izračuna gradijenta određena je kao vrijednost parametra jac funkcije minim, kao što je prikazano u nastavku.
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 (Newton)
Algoritam je modificirana Newtonova metoda.
Newtonova metoda temelji se na aproksimaciji funkcije u lokalnom području polinomom drugog stupnja:
gdje je matrica drugih izvoda (Hessian matrica, Hessian).
Ako je Hessian pozitivno određen, tada se lokalni minimum ove funkcije može pronaći izjednačavanjem nultog gradijenta kvadratne forme s nulom. Rezultat će biti izraz:
Inverzni Hessian izračunava se metodom konjugiranog gradijenta. Primjer korištenja ove metode za minimiziranje Rosenbrockove funkcije dan je 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 и , definirajte matricu .
Preostali elementi matrice koji nisu nula jednaki su:
Na primjer, u petodimenzionalnom prostoru N = 5 Hessova matrica za Rosenbrockovu funkciju ima oblik trake:
Kod koji izračunava ovaj Hessian zajedno s kodom za minimiziranje Rosenbrockove funkcije pomoću 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 s definicijom umnoška funkcije Hessiana i proizvoljnog vektora
U problemima 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 Hessianovu matricu, jer postupak minimizacije zahtijeva samo vektor jednak umnošku Hessiana s drugim proizvoljnim vektorom. Prema tome, s računalne točke gledišta, mnogo je bolje odmah definirati funkciju koja vraća rezultat umnoška Hessiana s proizvoljnim vektorom.
Razmotrimo funkciju Hess, koja uzima vektor minimizacije kao prvi argument, a proizvoljni vektor kao drugi argument (zajedno s ostalim argumentima funkcije koju treba minimizirati). U našem slučaju, izračunavanje umnoška Hessiana Rosenbrockove funkcije s proizvoljnim vektorom nije jako teško. Ako p je proizvoljni vektor, zatim produkt izgleda kao:
Funkcija koja izračunava umnožak 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 povjerene regije konjugiranog gradijenta (Newton)
Loše uvjetovanje Hessove matrice i netočni smjerovi pretraživanja mogu uzrokovati da Newtonov algoritam konjugiranog gradijenta bude neučinkovit. U takvim slučajevima prednost se daje (trust-region) konjugirani Newtonovi gradijenti.
Primjer s definicijom Hessove 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 umnoškom funkcije Hessiana 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)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
Poput trust-ncg metode, metode Krylovljevog tipa dobro su prikladne za rješavanje velikih problema jer koriste samo matrično-vektorske produkte. Njihova bit je riješiti problem u području povjerenja ograničenom krnjim Krylovljevim podprostorom. Za neizvjesne probleme bolje je koristiti ovu metodu, budući da koristi manji broj nelinearnih iteracija zbog manjeg broja matrica-vektor proizvoda po podproblemu, u usporedbi s trust-ncg metodom. Osim toga, rješenje kvadratnog podproblema nalazi se točnije nego korištenjem trust-ncg metode.
Primjer s definicijom Hessove 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 umnoškom funkcije Hessiana i proizvoljnog vektora:
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 pouzdanosti
Sve metode (Newton-CG, trust-ncg i trust-krylov) dobro su prikladne za rješavanje velikih problema (s tisućama varijabli). To je zbog činjenice da temeljni algoritam konjugiranog gradijenta podrazumijeva približno određivanje inverzne Hessove matrice. Rješenje se nalazi iterativno, bez eksplicitnog proširenja Hessiana. Budući da samo trebate definirati funkciju za umnožak Hessovog i proizvoljnog vektora, ovaj je algoritam posebno dobar za rad s rijetkim matricama (pojasna dijagonala). To osigurava niske troškove memorije i značajnu uštedu vremena.
Za probleme srednje veličine, trošak pohranjivanja i faktoringa Hessana nije kritičan. To znači da je moguće dobiti rješenje u manje ponavljanja, rješavajući gotovo točno podprobleme regije povjerenja. Da bi se to postiglo, neke se nelinearne jednadžbe rješavaju iterativno za svaki kvadratni podproblem. Takvo rješenje obično zahtijeva 3 ili 4 Choleskyjeve dekompozicije Hessianove matrice. Kao rezultat toga, metoda konvergira u manje ponavljanja i zahtijeva manje proračuna objektivne funkcije nego druge implementirane metode područja pouzdanosti. Ovaj algoritam uključuje samo određivanje potpune Hessianove matrice i ne podržava mogućnost upotrebe funkcije produkta Hessian i proizvoljnog vektora.
Primjer s minimizacijom 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.])Tu ćemo vjerojatno stati. U sljedećem ću članku pokušati ispričati najzanimljivije stvari o uvjetnoj minimizaciji, primjeni minimizacije u rješavanju problema aproksimacije, minimiziranju funkcije jedne varijable, proizvoljnim minimizatorima i pronalaženju korijena sustava jednadžbi pomoću scipy.optimize paket.
Izvor:
Izvor: www.habr.com