SciPy (uitgesproken als sai pie) is een wiskundig applicatiepakket gebaseerd op de Numpy Python-extensie. Met SciPy wordt uw interactieve Python-sessie dezelfde complete data science- en complexe systeemprototypingomgeving als MATLAB, IDL, Octave, R-Lab en SciLab. Vandaag wil ik het kort hebben over het gebruik van enkele bekende optimalisatie-algoritmen in het scipy.optimize-pakket. Meer gedetailleerde en actuele hulp bij het gebruik van functies kan altijd worden verkregen met behulp van de opdracht help() of met Shift+Tab.
Introductie
Om uzelf en uw lezers te behoeden voor het zoeken en lezen van primaire bronnen, zullen links naar beschrijvingen van methoden voornamelijk op Wikipedia staan. In de regel is deze informatie voldoende om de methoden in algemene termen en de voorwaarden voor hun toepassing te begrijpen. Om de essentie van wiskundige methoden te begrijpen, volgt u de links naar meer gezaghebbende publicaties, die u aan het einde van elk artikel of in uw favoriete zoekmachine kunt vinden.
De module scipy.optimize omvat dus de implementatie van de volgende procedures:
- Voorwaardelijke en onvoorwaardelijke minimalisatie van scalaire functies van verschillende variabelen (minim) met behulp van verschillende algoritmen (Nelder-Mead simplex, BFGS, Newton conjugaatgradiënten, и )
- Globale optimalisatie (bijvoorbeeld: , )
- Minimaliseren van reststromen (least_squares) en curve-fitting-algoritmen met behulp van niet-lineaire kleinste kwadraten (curve_fit)
- Minimaliseren van scalaire functies van één variabele (minim_scalar) en zoeken naar wortels (root_scalar)
- Multidimensionale oplossers van een systeem van vergelijkingen (wortel) met behulp van verschillende algoritmen (hybride Powell, of grootschalige methoden zoals ).
In dit artikel beschouwen we alleen het eerste item uit deze hele lijst.
Onvoorwaardelijke minimalisatie van een scalaire functie van verschillende variabelen
De minim-functie uit het scipy.optimize-pakket biedt een algemene interface voor het oplossen van voorwaardelijke en onvoorwaardelijke minimalisatieproblemen van scalaire functies van verschillende variabelen. Om te demonstreren hoe het werkt, hebben we een geschikte functie van verschillende variabelen nodig, die we op verschillende manieren zullen minimaliseren.
Voor deze doeleinden is de Rosenbrock-functie van N variabelen perfect, die de vorm heeft:
Ondanks het feit dat de Rosenbrock-functie en zijn Jacobi- en Hessische matrices (respectievelijk de eerste en tweede afgeleide) al zijn gedefinieerd in het scipy.optimize-pakket, zullen we deze zelf definiëren.
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)Laten we voor de duidelijkheid de waarden van de Rosenbrock-functie van twee variabelen in 3D tekenen.
Tekencode
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()
Van tevoren weten dat het minimum 0 is Laten we eens kijken naar voorbeelden van hoe we de minimumwaarde van de Rosenbrock-functie kunnen bepalen met behulp van verschillende scipy.optimize-procedures.
Nelder-Mead simplex-methode
Stel dat er een beginpunt x0 is in de vijfdimensionale ruimte. Laten we met behulp van het algoritme het minimumpunt van de Rosenbrock-functie vinden dat er het dichtst bij ligt (het algoritme wordt gespecificeerd als de waarde van de methodeparameter):
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.]De simplexmethode is de eenvoudigste manier om een expliciet gedefinieerde en redelijk vloeiende functie te minimaliseren. Het is niet nodig om de afgeleiden van een functie te berekenen; het is voldoende om alleen de waarden ervan te specificeren. De Nelder-Mead-methode is een goede keuze voor eenvoudige minimalisatieproblemen. Omdat er echter geen gradiëntschattingen worden gebruikt, kan het langer duren om het minimum te vinden.
Powell-methode
Een ander optimalisatie-algoritme waarbij alleen de functiewaarden worden berekend is . Om het te gebruiken, moet je method = 'powell' instellen in de minim-functie.
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
Om een snellere convergentie naar een oplossing te verkrijgen, is de procedure gebruikt de gradiënt van de objectieve functie. De gradiënt kan worden gespecificeerd als een functie of worden berekend met behulp van eerste-ordeverschillen. Hoe dan ook vereist de BFGS-methode doorgaans minder functieaanroepen dan de simplex-methode.
Laten we de afgeleide van de Rosenbrock-functie in analytische vorm vinden:
Deze uitdrukking is geldig voor de afgeleiden van alle variabelen behalve de eerste en de laatste, die zijn gedefinieerd als:
Laten we eens kijken naar de Python-functie die deze gradiënt berekent:
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 derDe gradiëntberekeningsfunctie wordt gespecificeerd als de waarde van de jac-parameter van de minim-functie, zoals hieronder weergegeven.
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]Conjugaatgradiëntalgoritme (Newton)
Het algoritme is een aangepaste methode van Newton.
De methode van Newton is gebaseerd op het benaderen van een functie in een lokaal gebied met een polynoom van de tweede graad:
waar is de matrix van tweede afgeleiden (Hessische matrix, Hessische).
Als de Hessiaan positief definitief is, kan het lokale minimum van deze functie worden gevonden door de nulgradiënt van de kwadratische vorm gelijk te stellen aan nul. Het resultaat is de uitdrukking:
De inverse Hessiaan wordt berekend met behulp van de geconjugeerde gradiëntmethode. Hieronder wordt een voorbeeld gegeven van het gebruik van deze methode om de Rosenbrock-functie te minimaliseren. Als u de Newton-CG-methode wilt gebruiken, moet u een functie opgeven die de Hessiaan berekent.
De Hessiaan van de Rosenbrock-functie in analytische vorm is gelijk aan:
waar и , definieer de matrix .
De overige niet-nul elementen van de matrix zijn gelijk aan:
In een vijfdimensionale ruimte N = 5 heeft de Hessische matrix voor de Rosenbrock-functie bijvoorbeeld de vorm van een band:
Code die deze Hessiaan berekent, samen met code voor het minimaliseren van de Rosenbrock-functie met behulp van de conjugaatgradiëntmethode (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]Een voorbeeld met de definitie van de productfunctie van de Hessische en een willekeurige vector
Bij problemen in de echte wereld kan het berekenen en opslaan van de gehele Hessische matrix aanzienlijke tijd- en geheugenbronnen vergen. In dit geval is het eigenlijk niet nodig om de Hessische matrix zelf te specificeren, omdat de minimalisatieprocedure vereist alleen een vector die gelijk is aan het product van de Hessiaan met een andere willekeurige vector. Vanuit computationeel oogpunt verdient het dus veel de voorkeur om onmiddellijk een functie te definiëren die het resultaat retourneert van het product van de Hessiaan met een willekeurige vector.
Beschouw de hess-functie, die de minimalisatievector als eerste argument neemt, en een willekeurige vector als tweede argument (samen met andere argumenten van de functie die moet worden geminimaliseerd). In ons geval is het berekenen van het product van de Hessiaan van de Rosenbrock-functie met een willekeurige vector niet erg moeilijk. Als p een willekeurige vector is, dan is het product heeft de vorm:
De functie die het product van de Hessische en een willekeurige vector berekent, wordt doorgegeven als de waarde van het hessp-argument aan de minimalisatiefunctie:
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: 66Algoritme voor geconjugeerde gradiëntvertrouwensregio (Newton)
Slechte conditionering van de Hessische matrix en onjuiste zoekrichtingen kunnen ervoor zorgen dat het conjugaatgradiëntalgoritme van Newton ineffectief is. In dergelijke gevallen wordt de voorkeur gegeven (vertrouwensregio) conjugeerde Newton-gradiënten.
Voorbeeld met de definitie van de Hessische 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.]Voorbeeld met de productfunctie van de Hessische en een willekeurige vector:
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.]Methoden van het Krylov-type
Net als de trust-ncg-methode zijn methoden van het Krylov-type zeer geschikt voor het oplossen van grootschalige problemen, omdat ze alleen matrix-vectorproducten gebruiken. Hun essentie is het oplossen van een probleem in een vertrouwensgebied dat wordt begrensd door een afgeknotte Krylov-subruimte. Voor onzekere problemen is het beter om deze methode te gebruiken, omdat deze een kleiner aantal niet-lineaire iteraties gebruikt vanwege het kleinere aantal matrix-vectorproducten per subprobleem, vergeleken met de trust-ncg-methode. Bovendien wordt de oplossing voor het kwadratische deelprobleem nauwkeuriger gevonden dan met behulp van de trust-ncg-methode.
Voorbeeld met de definitie van de Hessische 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.]
Voorbeeld met de productfunctie van de Hessische en een willekeurige vector:
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 voor geschatte oplossing in het vertrouwensgebied
Alle methoden (Newton-CG, trust-ncg en trust-krylov) zijn zeer geschikt voor het oplossen van grootschalige problemen (met duizenden variabelen). Dit is te wijten aan het feit dat het onderliggende geconjugeerde gradiëntalgoritme een geschatte bepaling van de inverse Hessische matrix impliceert. De oplossing wordt iteratief gevonden, zonder expliciete uitbreiding van de Hessiaan. Omdat u alleen een functie hoeft te definiëren voor het product van een Hessische en een willekeurige vector, is dit algoritme vooral goed voor het werken met dunne (banddiagonaal) matrices. Dit zorgt voor lage geheugenkosten en aanzienlijke tijdsbesparing.
Voor middelgrote problemen zijn de kosten voor het opslaan en factoriseren van de Hessiaan niet van cruciaal belang. Dit betekent dat het mogelijk is om in minder iteraties een oplossing te verkrijgen, waarbij de deelproblemen van de vertrouwensregio vrijwel exact worden opgelost. Om dit te doen, worden enkele niet-lineaire vergelijkingen iteratief opgelost voor elk kwadratisch deelprobleem. Een dergelijke oplossing vereist gewoonlijk 3 of 4 Cholesky-ontledingen van de Hessische matrix. Als gevolg hiervan convergeert de methode in minder iteraties en vereist deze minder objectieve functieberekeningen dan andere geïmplementeerde methoden voor betrouwbaarheidsregio's. Dit algoritme impliceert alleen de bepaling van de volledige Hessische matrix en ondersteunt niet de mogelijkheid om de productfunctie van de Hessische en een willekeurige vector te gebruiken.
Voorbeeld met minimalisatie van de Rosenbrock-functie:
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.])Waarschijnlijk stoppen we daar. In het volgende artikel zal ik proberen de meest interessante dingen te vertellen over voorwaardelijke minimalisatie, de toepassing van minimalisatie bij het oplossen van benaderingsproblemen, het minimaliseren van een functie van één variabele, willekeurige minimalisatoren en het vinden van de wortels van een stelsel van vergelijkingen met behulp van de scipy.optimize pakket.
Bron:
Bron: www.habr.com