SciPy (udtales sai pie) er en matematisk applikationspakke baseret på Numpy Python-udvidelsen. Med SciPy bliver din interaktive Python-session det samme komplette datavidenskabelige og komplekse systemprototypemiljø som MATLAB, IDL, Octave, R-Lab og SciLab. I dag vil jeg kort fortælle om, hvordan man bruger nogle velkendte optimeringsalgoritmer i scipy.optimize-pakken. Mere detaljeret og opdateret hjælp til brug af funktioner kan altid fås ved at bruge help()-kommandoen eller ved at bruge Shift+Tab.
Indledning
For at spare dig selv og læserne fra at søge og læse primære kilder, vil links til beskrivelser af metoder hovedsageligt være på Wikipedia. Som regel er disse oplysninger tilstrækkelige til at forstå metoderne i generelle vilkår og betingelserne for deres anvendelse. For at forstå essensen af matematiske metoder, følg linkene til mere autoritative publikationer, som kan findes i slutningen af hver artikel eller i din foretrukne søgemaskine.
Så scipy.optimize-modulet inkluderer implementeringen af følgende procedurer:
- Betinget og ubetinget minimering af skalarfunktioner af flere variable (minimum) ved hjælp af forskellige algoritmer (Nelder-Mead simplex, BFGS, Newton konjugerede gradienter, и )
- Global optimering (f.eks.: , )
- Minimering af rester (mindste kvadrater) og kurvetilpasningsalgoritmer ved hjælp af ikke-lineære mindste kvadrater (kurvetilpasning)
- Minimering af skalarfunktioner af en variabel (minim_scalar) og søgning efter rødder (root_scalar)
- Multidimensionelle løsere af ligningssystem (rod) ved hjælp af forskellige algoritmer (hybrid Powell, eller storskala metoder som f.eks ).
I denne artikel vil vi kun overveje det første element fra hele denne liste.
Ubetinget minimering af en skalarfunktion af flere variable
Minimfunktionen fra scipy.optimize-pakken giver en generel grænseflade til løsning af betingede og ubetingede minimeringsproblemer af skalarfunktioner af flere variable. For at demonstrere, hvordan det virker, skal vi bruge en passende funktion af flere variable, som vi vil minimere på forskellige måder.
Til disse formål er Rosenbrock-funktionen af N variable perfekt, som har formen:
På trods af at Rosenbrock-funktionen og dens Jacobi- og Hessian-matricer (henholdsvis den første og anden afledte) allerede er defineret i scipy.optimize-pakken, vil vi selv definere den.
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)For klarhedens skyld, lad os tegne i 3D værdierne af Rosenbrock-funktionen af to variable.
Tegningskode
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()
På forhånd at vide, at minimum er 0 kl , lad os se på eksempler på, hvordan man bestemmer minimumsværdien af Rosenbrock-funktionen ved hjælp af forskellige scipy.optimize-procedurer.
Nelder-Mead simplex metode
Lad der være et begyndelsespunkt x0 i 5-dimensionelt rum. Lad os finde minimumspunktet for Rosenbrock-funktionen tættest på den ved hjælp af algoritmen (algoritmen er angivet som værdien af metodeparameteren):
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.]Simplex-metoden er den enkleste måde at minimere en eksplicit defineret og ret glat funktion. Det kræver ikke at beregne afledte af en funktion; det er nok kun at angive dens værdier. Nelder-Mead metoden er et godt valg til simple minimeringsproblemer. Men da den ikke bruger gradientestimater, kan det tage længere tid at finde minimum.
Powell metode
En anden optimeringsalgoritme, hvor kun funktionsværdierne beregnes, er . For at bruge det skal du indstille metode = 'powell' i minimfunktionen.
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) algoritme
For at opnå hurtigere konvergens til en løsning, proceduren bruger gradienten af objektivfunktionen. Gradienten kan angives som en funktion eller beregnes ved hjælp af førsteordens forskelle. Under alle omstændigheder kræver BFGS-metoden typisk færre funktionskald end simplex-metoden.
Lad os finde den afledte af Rosenbrock-funktionen i analytisk form:
Dette udtryk er gyldigt for de afledte variabler undtagen den første og den sidste, som er defineret som:
Lad os se på Python-funktionen, der beregner denne 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 derGradientberegningsfunktionen er angivet som værdien af jac-parameteren for minimfunktionen, som vist nedenfor.
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]Konjugeret gradientalgoritme (Newton)
Algoritme er en modificeret Newtons metode.
Newtons metode er baseret på at tilnærme en funktion i et lokalområde ved et polynomium af anden grad:
где er matrixen af anden afledte (hessisk matrix, hessisk).
Hvis hessian er positiv bestemt, så kan det lokale minimum af denne funktion findes ved at sidestille nulgradienten af den kvadratiske form med nul. Resultatet bliver udtrykket:
Den omvendte hessian beregnes ved hjælp af den konjugerede gradientmetode. Et eksempel på brug af denne metode til at minimere Rosenbrock-funktionen er givet nedenfor. For at bruge Newton-CG metoden skal du angive en funktion, der beregner hessian.
Hessian af Rosenbrock-funktionen i analytisk form er lig med:
где и , definere matrixen .
De resterende ikke-nul elementer i matrixen er lig med:
For eksempel, i et femdimensionelt rum N = 5, har den hessiske matrix for Rosenbrock-funktionen form af et bånd:
Kode, der beregner denne hessian sammen med kode til minimering af Rosenbrock-funktionen ved hjælp af konjugeret gradient (Newton) metode:
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]Et eksempel med definition af produktfunktionen af hessian og en vilkårlig vektor
Ved problemer i den virkelige verden kan det kræve betydelige tids- og hukommelsesressourcer at beregne og lagre hele den hessiske matrix. I dette tilfælde er der faktisk ikke behov for at specificere selve den hessiske matrix, fordi minimeringsproceduren kræver kun en vektor svarende til produktet af hessian med en anden vilkårlig vektor. Fra et beregningsmæssigt synspunkt er det således meget at foretrække umiddelbart at definere en funktion, der returnerer resultatet af produktet af hessian med en vilkårlig vektor.
Overvej hess-funktionen, som tager minimeringsvektoren som det første argument, og en vilkårlig vektor som det andet argument (sammen med andre argumenter for funktionen, der skal minimeres). I vores tilfælde er det ikke særlig svært at beregne produktet af Hessian af Rosenbrock-funktionen med en vilkårlig vektor. Hvis p er en vilkårlig vektor, så produktet har formen:
Funktionen, der beregner produktet af hessian og en vilkårlig vektor, overføres som værdien af hessp-argumentet til minimeringsfunktionen:
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: 66Konjugeret gradient-tillidsregionsalgoritme (Newton)
Dårlig konditionering af den hessiske matrix og forkerte søgeretninger kan forårsage, at Newtons konjugerede gradientalgoritme er ineffektiv. I sådanne tilfælde gives fortrinsret til (trust-region) konjugerer Newton-gradienter.
Eksempel med definitionen af den hessiske matrix:
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.]Eksempel med produktfunktionen af hessian og en vilkårlig vektor:
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.]Krylov-type metoder
Ligesom trust-ncg-metoden er metoder af Krylov-typen velegnede til at løse store problemer, fordi de kun bruger matrix-vektorprodukter. Deres essens er at løse et problem i et tillidsområde begrænset af et afkortet Krylov-underrum. For usikre problemer er det bedre at bruge denne metode, da den bruger et mindre antal ikke-lineære iterationer på grund af det mindre antal matrix-vektorprodukter pr. underproblem sammenlignet med trust-ncg-metoden. Derudover findes løsningen på det kvadratiske delproblem mere præcist end at bruge trust-ncg metoden.
Eksempel med definitionen af den hessiske matrix:
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.]
Eksempel med produktfunktionen af hessian og en vilkårlig vektor:
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.]
Algoritme til omtrentlig løsning i konfidensområdet
Alle metoder (Newton-CG, trust-ncg og trust-krylov) er velegnede til at løse store problemer (med tusindvis af variabler). Dette skyldes det faktum, at den underliggende konjugerede gradientalgoritme indebærer en omtrentlig bestemmelse af den inverse hessiske matrix. Løsningen findes iterativt uden eksplicit udvidelse af hessisk. Da du kun skal definere en funktion for produktet af en hessian og en vilkårlig vektor, er denne algoritme især god til at arbejde med sparsomme (bånddiagonale) matricer. Dette giver lave hukommelsesomkostninger og betydelige tidsbesparelser.
For mellemstore problemer er omkostningerne ved opbevaring og faktorisering af hessian ikke kritiske. Det betyder, at det er muligt at opnå en løsning i færre iterationer, hvilket løser tillidsregionens delproblemer næsten nøjagtigt. For at gøre dette løses nogle ikke-lineære ligninger iterativt for hver andengradsdelopgave. En sådan løsning kræver normalt 3 eller 4 Cholesky-nedbrydninger af den hessiske matrix. Som et resultat konvergerer metoden i færre iterationer og kræver færre objektive funktionsberegninger end andre implementerede konfidensregionmetoder. Denne algoritme involverer kun bestemmelse af den komplette hessiske matrix og understøtter ikke evnen til at bruge produktfunktionen af hessisk og en vilkårlig vektor.
Eksempel med minimering af Rosenbrock-funktionen:
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.])Vi stopper nok der. I den næste artikel vil jeg forsøge at fortælle de mest interessante ting om betinget minimering, anvendelsen af minimering til at løse tilnærmelsesproblemer, minimering af en funktion af én variabel, vilkårlige minimizere og at finde rødderne til et ligningssystem ved hjælp af scipy.optimize pakke.
Kilde:
Kilde: www.habr.com