SciPy (pronunțat sai pie) este un pachet de aplicații matematice bazat pe extensia Numpy Python. Cu SciPy, sesiunea dvs. interactivă Python devine același mediu complet de știință a datelor și același mediu complex de prototipare a sistemului ca MATLAB, IDL, Octave, R-Lab și SciLab. Astăzi vreau să vă vorbesc pe scurt despre cum să folosiți niște algoritmi de optimizare cunoscuți în pachetul scipy.optimize. Ajutor mai detaliat și actualizat cu privire la utilizarea funcțiilor poate fi întotdeauna obținut folosind comanda help() sau folosind Shift+Tab.
Introducere
Pentru a vă salva pe dumneavoastră și pe cititori de căutarea și citirea surselor primare, link-urile către descrierile metodelor vor fi în principal pe Wikipedia. De regulă, aceste informații sunt suficiente pentru a înțelege metodele în termeni generali și condițiile de aplicare a acestora. Pentru a înțelege esența metodelor matematice, urmărește linkurile către publicații mai autorizate, care pot fi găsite la sfârșitul fiecărui articol sau în motorul tău de căutare preferat.
Deci, modulul scipy.optimize include implementarea următoarelor proceduri:
- Minimizarea condiționată și necondiționată a funcțiilor scalare ale mai multor variabile (minim) folosind diverși algoritmi (Nelder-Mead simplex, BFGS, gradienți conjugați Newton, и )
- Optimizare globală (de exemplu: , )
- Minimizarea reziduurilor (least_squares) și algoritmi de potrivire a curbelor care utilizează cele mai mici pătrate neliniare (curve_fit)
- Minimizarea funcțiilor scalare ale unei variabile (minim_scalar) și căutarea rădăcinilor (root_scalar)
- Rezolvatori multidimensionali de sisteme de ecuații (rădăcină) folosind diverși algoritmi (hibrid Powell, sau metode pe scară largă precum ).
În acest articol vom lua în considerare doar primul articol din întreaga listă.
Minimizarea necondiționată a unei funcții scalare a mai multor variabile
Funcția minim din pachetul scipy.optimize oferă o interfață generală pentru rezolvarea problemelor de minimizare condiționată și necondiționată a funcțiilor scalare ale mai multor variabile. Pentru a demonstra cum funcționează, vom avea nevoie de o funcție adecvată a mai multor variabile, pe care le vom minimiza în diferite moduri.
În aceste scopuri, funcția Rosenbrock a N variabile este perfectă, care are forma:
În ciuda faptului că funcția Rosenbrock și matricele sale Jacobi și Hessian (prima și, respectiv, a doua derivată) sunt deja definite în pachetul scipy.optimize, o vom defini noi înșine.
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)Pentru claritate, să desenăm în 3D valorile funcției Rosenbrock a două variabile.
Cod de desen
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()
Știind dinainte că minimul este 0 la , să ne uităm la exemple despre cum să determinați valoarea minimă a funcției Rosenbrock folosind diverse proceduri scipy.optimize.
Metoda Nelder-Mead simplex
Să existe un punct inițial x0 în spațiul cu 5 dimensiuni. Să găsim punctul minim al funcției Rosenbrock cel mai apropiat de acesta folosind algoritmul (algoritmul este specificat ca valoare a parametrului metodei):
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.]Metoda simplex este cea mai simplă modalitate de a minimiza o funcție definită în mod explicit și destul de netedă. Nu necesită calcularea derivatelor unei funcții; este suficient să specificați doar valorile acesteia. Metoda Nelder-Mead este o alegere bună pentru probleme simple de minimizare. Cu toate acestea, deoarece nu utilizează estimări de gradient, poate dura mai mult pentru a găsi minimul.
metoda Powell
Un alt algoritm de optimizare în care sunt calculate doar valorile funcției este . Pentru a-l folosi, trebuie să setați metoda = 'powell' în funcția 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.]Algoritmul Broyden-Fletcher-Goldfarb-Shanno (BFGS).
Pentru a obține o convergență mai rapidă către o soluție, procedura folosește gradientul funcției obiectiv. Gradientul poate fi specificat ca funcție sau calculat folosind diferențe de ordinul întâi. În orice caz, metoda BFGS necesită de obicei mai puține apeluri de funcție decât metoda simplex.
Să găsim derivata funcției Rosenbrock în formă analitică:
Această expresie este valabilă pentru derivatele tuturor variabilelor, cu excepția primei și ultimei, care sunt definite ca:
Să ne uităm la funcția Python care calculează acest 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 derFuncția de calcul a gradientului este specificată ca valoare a parametrului jac al funcției minim, așa cum se arată mai jos.
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]Algoritm de gradient conjugat (Newton)
Algoritmul este o metodă modificată a lui Newton.
Metoda lui Newton se bazează pe aproximarea unei funcții dintr-o zonă locală printr-un polinom de gradul doi:
unde este matricea derivatelor secunde (matrice Hessian, Hessian).
Dacă Hessianul este definit pozitiv, atunci minimul local al acestei funcții poate fi găsit prin echivalarea gradientului zero al formei pătratice cu zero. Rezultatul va fi expresia:
Hessianul invers este calculat folosind metoda gradientului conjugat. Un exemplu de utilizare a acestei metode pentru a minimiza funcția Rosenbrock este dat mai jos. Pentru a utiliza metoda Newton-CG, trebuie să specificați o funcție care calculează Hessianul.
Hessianul funcției Rosenbrock în formă analitică este egal cu:
unde и , definiți matricea .
Elementele rămase nenule ale matricei sunt egale cu:
De exemplu, într-un spațiu cu cinci dimensiuni N = 5, matricea Hessiană pentru funcția Rosenbrock are forma unei benzi:
Cod care calculează acest Hessian împreună cu codul pentru minimizarea funcției Rosenbrock folosind metoda gradientului conjugat (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]Un exemplu cu definiția funcției de produs a Hessianului și a unui vector arbitrar
În problemele din lumea reală, calcularea și stocarea întregii matrice Hessian pot necesita resurse semnificative de timp și memorie. În acest caz, de fapt, nu este nevoie să specificați matricea Hessian în sine, deoarece procedura de minimizare necesită doar un vector egal cu produsul lui Hessian cu un alt vector arbitrar. Astfel, din punct de vedere computațional, este mult de preferat să definiți imediat o funcție care returnează rezultatul produsului Hessianului cu un vector arbitrar.
Luați în considerare funcția hess, care ia vectorul de minimizare ca prim argument și un vector arbitrar ca al doilea argument (împreună cu alte argumente ale funcției de minimizat). În cazul nostru, calcularea produsului Hessian al funcției Rosenbrock cu un vector arbitrar nu este foarte dificilă. Dacă p este un vector arbitrar, apoi produsul se pare ca:
Funcția care calculează produsul dintre Hessian și un vector arbitrar este transmisă ca valoare a argumentului hessp la funcția de minimizare:
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: 66Algoritmul conjugat pentru regiunea de încredere a gradientului (Newton)
Condiționarea proastă a matricei Hessian și direcțiile de căutare incorecte pot face ca algoritmul de gradient conjugat al lui Newton să fie ineficient. În astfel de cazuri, se acordă preferință (încredere-regiune) gradienți Newton conjugați.
Exemplu cu definiția matricei Hessian:
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.]Exemplu cu funcția de produs a Hessianului și a unui vector arbitrar:
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 de tip Krylov
La fel ca metoda trust-ncg, metodele de tip Krylov sunt potrivite pentru rezolvarea problemelor la scară mare, deoarece folosesc numai produse matrice-vectorale. Esența lor este de a rezolva o problemă într-o regiune de încredere limitată de un subspațiu Krylov trunchiat. Pentru problemele incerte, este mai bine să folosiți această metodă, deoarece folosește un număr mai mic de iterații neliniare datorită numărului mai mic de produse matrice-vector per subproblemă, comparativ cu metoda trust-ncg. În plus, soluția subproblemei pătratice este găsită cu mai multă acuratețe decât folosind metoda trust-ncg.
Exemplu cu definiția matricei Hessian:
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.]
Exemplu cu funcția de produs a Hessianului și a unui vector arbitrar:
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.]
Algoritm pentru soluția aproximativă în regiunea de încredere
Toate metodele (Newton-CG, trust-ncg și trust-krylov) sunt potrivite pentru rezolvarea problemelor la scară largă (cu mii de variabile). Acest lucru se datorează faptului că algoritmul de gradient conjugat de bază implică o determinare aproximativă a matricei hessiene inverse. Soluția se găsește iterativ, fără extinderea explicită a hessianului. Deoarece trebuie doar să definiți o funcție pentru produsul dintre un vector Hessian și un vector arbitrar, acest algoritm este deosebit de bun pentru a lucra cu matrici rare (diagonale de bandă). Acest lucru asigură costuri reduse de memorie și economii semnificative de timp.
Pentru problemele de dimensiuni medii, costul de stocare și factoring Hessian nu este critic. Aceasta înseamnă că este posibil să se obțină o soluție în mai puține iterații, rezolvând aproape exact subproblemele regiunii de încredere. Pentru a face acest lucru, unele ecuații neliniare sunt rezolvate iterativ pentru fiecare subproblemă pătratică. O astfel de soluție necesită de obicei 3 sau 4 descompuneri Cholesky ale matricei Hessian. Ca rezultat, metoda converge în mai puține iterații și necesită mai puține calcule ale funcției obiective decât alte metode implementate de regiune de încredere. Acest algoritm implică doar determinarea matricei Hessian complete și nu acceptă capacitatea de a utiliza funcția de produs a Hessianului și un vector arbitrar.
Exemplu cu minimizarea funcției 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.])Probabil ne vom opri aici. În articolul următor voi încerca să spun cele mai interesante lucruri despre minimizarea condiționată, aplicarea minimizării în rezolvarea problemelor de aproximare, minimizarea unei funcții a unei variabile, minimizatorii arbitrari și găsirea rădăcinilor unui sistem de ecuații folosind scipy.optimize. pachet.
Sursa:
Sursa: www.habr.com