SciPy (wymawiane sai pie) to pakiet aplikacji matematycznych oparty na rozszerzeniu Numpy Python. Dzięki SciPy Twoja interaktywna sesja Pythona staje się tym samym kompletnym środowiskiem do nauki o danych i złożonym prototypowaniem systemów, co MATLAB, IDL, Octave, R-Lab i SciLab. Dzisiaj chcę krótko porozmawiać o tym, jak wykorzystać niektóre dobrze znane algorytmy optymalizacyjne w pakiecie scipy.optimize. Bardziej szczegółową i aktualną pomoc dotyczącą korzystania z funkcji można zawsze uzyskać za pomocą polecenia help() lub kombinacji Shift+Tab.
Wprowadzenie
Aby oszczędzić sobie i czytelnikom poszukiwania i czytania źródeł pierwotnych, linki do opisów metod będą znajdować się głównie w Wikipedii. Z reguły informacje te są wystarczające, aby ogólnie zrozumieć metody i warunki ich stosowania. Aby zrozumieć istotę metod matematycznych, skorzystaj z linków do bardziej miarodajnych publikacji, które znajdziesz na końcu każdego artykułu lub w Twojej ulubionej wyszukiwarce.
Zatem moduł scipy.optimize obejmuje implementację następujących procedur:
- Minimalizacja warunkowa i bezwarunkowa funkcji skalarnych kilku zmiennych (minim) przy użyciu różnych algorytmów (simpleks Neldera-Meada, BFGS, gradienty sprzężone Newtona, и )
- Globalna optymalizacja (na przykład: , )
- Minimalizacja pozostałości (najmniejsze kwadraty) i algorytmy dopasowywania krzywych wykorzystujące nieliniową metodę najmniejszych kwadratów (curve_fit)
- Minimalizowanie funkcji skalarnych jednej zmiennej (minim_scalar) i wyszukiwanie pierwiastków (root_scalar)
- Wielowymiarowe rozwiązania układu równań (pierwiastek) wykorzystujące różne algorytmy (hybrydowe Powella, lub metody na dużą skalę, takie jak ).
W tym artykule rozważymy tylko pierwszą pozycję z całej listy.
Bezwarunkowa minimalizacja funkcji skalarnej kilku zmiennych
Funkcja minim z pakietu scipy.optimize zapewnia ogólny interfejs do rozwiązywania problemów minimalizacji warunkowej i bezwarunkowej funkcji skalarnych kilku zmiennych. Aby zademonstrować jak to działa, będziemy potrzebować odpowiedniej funkcji kilku zmiennych, które będziemy minimalizować na różne sposoby.
Do tych celów doskonale sprawdza się funkcja Rosenbrocka N zmiennych, która ma postać:
Pomimo tego, że funkcja Rosenbrocka oraz jej macierze Jacobiego i Hessego (odpowiednio pierwsza i druga pochodna) są już zdefiniowane w pakiecie scipy.optimize, zdefiniujemy ją 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)Dla przejrzystości narysujmy w 3D wartości funkcji Rosenbrocka dwóch zmiennych.
Kod rysunkowy
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()
Wiedząc z góry, że minimum wynosi 0 w , spójrzmy na przykłady wyznaczania minimalnej wartości funkcji Rosenbrocka przy użyciu różnych procedur scipy.optimize.
Metoda simpleks Neldera-Meada
Niech istnieje punkt początkowy x0 w przestrzeni 5-wymiarowej. Znajdźmy najbliższy jej minimalny punkt funkcji Rosenbrocka, korzystając z algorytmu (algorytm podawany jest jako wartość parametru metody):
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 to najprostszy sposób minimalizacji jawnie zdefiniowanej i dość gładkiej funkcji. Nie wymaga obliczania pochodnych funkcji, wystarczy podać jej wartości. Metoda Neldera-Meada jest dobrym wyborem w przypadku prostych problemów minimalizacji. Ponieważ jednak nie wykorzystuje szacunków gradientu, znalezienie minimum może zająć więcej czasu.
Metoda Powella
Innym algorytmem optymalizacji, w którym obliczane są tylko wartości funkcji, jest . Aby z niego skorzystać, musisz ustawić metodę = 'powell' w funkcji 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.]Algorytm Broydena-Fletchera-Goldfarba-Shanno (BFGS).
Aby uzyskać szybszą zbieżność do rozwiązania, należy zastosować procedurę wykorzystuje gradient funkcji celu. Gradient można określić jako funkcję lub obliczyć przy użyciu różnic pierwszego rzędu. W każdym razie metoda BFGS zazwyczaj wymaga mniejszej liczby wywołań funkcji niż metoda simpleks.
Znajdźmy pochodną funkcji Rosenbrocka w postaci analitycznej:
To wyrażenie obowiązuje dla pochodnych wszystkich zmiennych z wyjątkiem pierwszej i ostatniej, które są zdefiniowane jako:
Spójrzmy na funkcję Pythona, która oblicza ten 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 derFunkcja obliczania gradientu jest określona jako wartość parametru jac funkcji minim, jak pokazano poniżej.
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]Algorytm gradientu sprzężonego (Newton)
Algorytm jest zmodyfikowaną metodą Newtona.
Metoda Newtona polega na aproksymacji funkcji w obszarze lokalnym wielomianem drugiego stopnia:
gdzie jest macierzą drugich pochodnych (macierz Hesja, Hesja).
Jeśli Hesjan jest dodatnio określony, to minimum lokalne tej funkcji można znaleźć, przyrównując zerowy gradient postaci kwadratowej do zera. Wynikiem będzie wyrażenie:
Odwrotność Hesja oblicza się metodą gradientu sprzężonego. Poniżej podano przykład zastosowania tej metody do minimalizacji funkcji Rosenbrocka. Aby zastosować metodę Newtona-CG, należy określić funkcję obliczającą wartość Hesja.
Hesjan funkcji Rosenbrocka w formie analitycznej jest równy:
gdzie и , zdefiniuj macierz .
Pozostałe niezerowe elementy macierzy są równe:
Przykładowo w przestrzeni pięciowymiarowej N=5 macierz Hessego dla funkcji Rosenbrocka ma postać pasma:
Kod obliczający ten Hesjan wraz z kodem minimalizującym funkcję Rosenbrocka przy użyciu metody gradientu sprzężonego (Newtona):
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]Przykład z definicją funkcji iloczynu Hesja i dowolnego wektora
W przypadku problemów występujących w świecie rzeczywistym obliczenie i przechowywanie całej macierzy Hessego może wymagać znacznych zasobów czasu i pamięci. W tym przypadku właściwie nie ma potrzeby określania samej macierzy Hessego, ponieważ procedura minimalizacji wymaga jedynie wektora równego iloczynowi Hesja z innym dowolnym wektorem. Zatem z obliczeniowego punktu widzenia znacznie lepiej jest natychmiast zdefiniować funkcję, która zwraca wynik iloczynu Hessego za pomocą dowolnego wektora.
Rozważmy funkcję hessa, która jako pierwszy argument przyjmuje wektor minimalizacji, a jako drugi argument dowolny wektor (wraz z innymi argumentami funkcji, która ma zostać zminimalizowana). W naszym przypadku obliczenie iloczynu Hesja funkcji Rosenbrocka za pomocą dowolnego wektora nie jest bardzo trudne. Jeśli p jest dowolnym wektorem, to iloczyn ma postać:
Funkcja obliczająca iloczyn hesjanu i dowolnego wektora jest przekazywana jako wartość argumentu hessp do funkcji minimalizacji:
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: 66Algorytm regionu zaufania z gradientem sprzężonym (Newton)
Złe uwarunkowanie macierzy Hessego i nieprawidłowe kierunki poszukiwań mogą spowodować, że algorytm gradientu sprzężonego Newtona będzie nieskuteczny. W takich przypadkach preferowane jest (region zaufania) sprzężenie gradientów Newtona.
Przykład z definicją macierzy Hessego:
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.]Przykład z funkcją iloczynu Hesja i dowolnego wektora:
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.]Metody typu Kryłowa
Podobnie jak metoda trust-ncg, metody typu Kryłowa dobrze nadają się do rozwiązywania problemów na dużą skalę, ponieważ wykorzystują wyłącznie iloczyny macierzowo-wektorowe. Ich istotą jest rozwiązanie problemu w obszarze ufności ograniczonym przez obciętą podprzestrzeń Kryłowa. W przypadku problemów niepewnych lepiej jest zastosować tę metodę, gdyż wykorzystuje ona mniejszą liczbę nieliniowych iteracji ze względu na mniejszą liczbę iloczynów wektorów macierzowych na podproblem w porównaniu z metodą trust-ncg. Ponadto rozwiązanie podproblemu kwadratowego znajduje się dokładniej niż przy użyciu metody trust-ncg.
Przykład z definicją macierzy Hessego:
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.]
Przykład z funkcją iloczynu Hesja i dowolnego wektora:
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.]
Algorytm rozwiązania przybliżonego w obszarze ufności
Wszystkie metody (Newton-CG, trust-ncg i trust-krylov) dobrze nadają się do rozwiązywania problemów na dużą skalę (z tysiącami zmiennych). Wynika to z faktu, że leżący u podstaw algorytm gradientu sprzężonego implikuje przybliżone określenie odwrotnej macierzy Hessego. Rozwiązanie znajduje się iteracyjnie, bez wyraźnego rozszerzania Hesjanu. Ponieważ wystarczy zdefiniować funkcję dla iloczynu Hessego i dowolnego wektora, algorytm ten jest szczególnie dobry do pracy z macierzami rzadkimi (przekątna pasma). Zapewnia to niskie koszty pamięci i znaczną oszczędność czasu.
W przypadku problemów średniej wielkości koszt przechowywania i faktoringu Hesji nie jest krytyczny. Oznacza to, że możliwe jest uzyskanie rozwiązania w mniejszej liczbie iteracji, rozwiązując niemal dokładnie podproblemy obszaru zaufania. Aby to zrobić, niektóre równania nieliniowe są rozwiązywane iteracyjnie dla każdego podproblemu kwadratowego. Takie rozwiązanie wymaga zwykle 3 lub 4 rozkładów Choleskiego macierzy Hessego. W rezultacie metoda osiąga zbieżność w mniejszej liczbie iteracji i wymaga mniej obliczeń funkcji celu niż inne zaimplementowane metody obszaru ufności. Algorytm ten polega jedynie na wyznaczeniu pełnej macierzy Hessego i nie umożliwia wykorzystania funkcji iloczynu Hessego i dowolnego wektora.
Przykład z minimalizacją funkcji Rosenbrocka:
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.])Prawdopodobnie na tym się zatrzymamy. W kolejnym artykule postaram się opowiedzieć najciekawsze rzeczy o minimalizacji warunkowej, zastosowaniu minimalizacji w rozwiązywaniu problemów aproksymacyjnych, minimalizacji funkcji jednej zmiennej, arbitralnych minimalizatorach oraz znajdowaniu pierwiastków układu równań za pomocą funkcji scipy.optimize pakiet.
Źródło:
Źródło: www.habr.com