ProHoster > Blog > Nowinki internetowe > SciPy, optymalizacja

SciPy, optymalizacja

SciPy, optymalizacja

SciPy (wymawia się jak saj paj) to pakiet aplikacji matematycznych oparty na rozszerzeniu Numpy Python. Z SciPy interaktywny sesja Python przekształca się w pełnoprawne środowisko przetwarzania danych oraz prototypowania złożonych systemów, takie jak MATLAB, IDL, Octave, R-Lab i SciLab. Dziś chciałbym krótko omówić, jak zastosować niektóre znane algorytmy optymalizacji w pakiecie scipy.optimize. Szczegółowe i aktualne informacje na temat zastosowania funkcji można zawsze uzyskać za pomocą polecenia help() lub kombinacji Shift+Tab.

Wprowadzenie

Aby zaoszczędzić sobie i czytelnikom poszukiwania oraz czytania źródłowych materiałów, odniesienia do opisu metod będą głównie odsyłać do Wikipedii. Zazwyczaj te informacje są wystarczające do ogólnego zrozumienia metod oraz warunków ich zastosowania. Aby zrozumieć istotę metod matematycznych, przechodzimy do bardziej autorytatywnych publikacji, które można znaleźć na końcu każdego artykułu lub w ulubionej wyszukiwarce.

Moduł scipy.optimize obejmuje implementację następujących procedur:

  1. Warunkowej i bezwarunkowej minimalizacji funkcji skalarnej wielu zmiennych (minim) przy użyciu różnych algorytmów (simplex Neldera-Meada, BFGS, sprzężone gradienty Newtona, COBYLA i SLSQP)
  2. Globalnej optymalizacji (na przykład: basinhopping, diff_evolution)
  3. Minimalizacji reszt MNK (least_squares) oraz algorytmy dopasowania krzywych nieliniowym MNK (curve_fit)
  4. Minimalizacji funkcji skalarnej jednej zmiennej (minim_scalar) oraz znajdowania miejsc zerowych (root_scalar)
  5. Wielo wymiarowych rozwiązań systemu równań (root) z wykorzystaniem różnych algorytmów (hybrydowy Powell, Levenberga-Marquardta lub metod wielkoskalowych, takich jak Newtona-Kryłowa.).

W tym artykule omówimy tylko pierwszy punkt z tej listy.

Bezwarunkowa minimalizacja funkcji skalarnej wielu zmiennych

Funkcja minim z pakietu scipy.optimize zapewnia ogólny interfejs do rozwiązywania zadań warunkowej i bezwarunkowej minimalizacji funkcji skalarnej wielu zmiennych. Aby zilustrować jej działanie, potrzebujemy odpowiedniej funkcji wielu zmiennych, którą będziemy minimalizować na różne sposoby.

Do tego celu idealnie nadaje się funkcja Rosenbrocka od N zmiennych, która ma postać:

SciPy, optymalizacja

Mimo że funkcja Rosenbrocka oraz jej macierze Jacobiego i Hesjona (pierwszej i drugiej pochodnej odpowiednio) są już zdefiniowane w pakiecie scipy.optimize, zdefiniujmy ją samodzielnie.

import numpy as np

def rosen(x):
    """Funkcja Rosenbrocka"""
    return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

Dla lepszej wizualizacji narysujmy w 3D wartości funkcji Rosenbrocka w dwóch zmiennych.

Kod do rysowania

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter

# Ustawiamy wykres 3D
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')

# Ustawiamy kąt widzenia
ax.view_init(45, 30)

# Tworzymy dane do wykresu
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]))

# Rysujemy powierzchnię
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()

SciPy, optymalizacja

Z góry wiedząc, że minimum wynosi 0, przyjrzymy się przykładom tego, jak określić minimalną wartość funkcji Rosenbrocka za pomocą różnych procedur scipy.optimize. SciPy, optymalizacja, rozważmy przykłady tego, jak znaleźć minimum funkcji Rosenbrocka używając różnych procedur scipy.optimize.

Metoda sympleksowa Neldera-Meada

Niech będzie dana początkowa punkt x0 w pięcio wymiarowej przestrzeni. Znajdziemy najbliższy punkt minimum funkcji Rosenbrocka za pomocą algorytmu sympleksu Neldera-Meada (algorytm podany jako wartość parametru method):

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)

Optymalizacja zakończona pomyślnie.
         Bieżąca wartość funkcji: 0.000000
         Iteracje: 339
         Oceny funkcji: 571
[1. 1. 1. 1. 1.]

Metoda sympleksowa jest najprostszym sposobem minimalizacji wyraźnie określonej i zdrowo gładkiej funkcji. Nie wymaga obliczania pochodnych funkcji, wystarczy podać tylko jej wartości. Metoda Neldera-Meada jest dobrym wyborem dla prostych zadań minimalizacji. Jednak ponieważ nie używa oceny gradientów, może wymagać więcej czasu, aby znaleźć minimum.

Metoda Powella

Innym algorytmem optymalizacji, w którym obliczane są tylko wartości funkcji, jest metoda Powella. Aby z niej skorzystać, należy ustawić method = ‘powell’ w funkcji minimize.

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)

Optymalizacja zakończona pomyślnie.
         Bieżąca wartość funkcji: 0.000000
         Iteracje: 19
         Oceny funkcji: 1622
[1. 1. 1. 1. 1.]

Algorytm Broydena-Fletcher-Goldfarba-Shanno (BFGS)

Aby uzyskać szybszą zbieżność do rozwiązania, procedura BFGS używa gradientu funkcji celu. Gradient może być zdefiniowany jako funkcja lub obliczany za pomocą różnic pierwszego rzędu. W każdym razie metoda BFGS zazwyczaj wymaga mniej wywołań funkcji niż metoda sympleksowa.

Znajdźmy pochodną z funkcji Rosenbrocka w postaci analitycznej:

SciPy, optymalizacja

SciPy, optymalizacja

To wyrażenie jest prawdziwe dla pochodnych wszystkich zmiennych, z wyjątkiem pierwszej i ostatniej, które definiowane są jako:

SciPy, optymalizacja

SciPy, optymalizacja

Spójrzmy na funkcję Python, 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 der

Funkcja obliczania gradientu jest podawana jako wartość parametru jac w funkcji minim, jak pokazano poniżej.

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)

Optymalizacja zakończona pomyślnie.
         Bieżąca wartość funkcji: 0.000000
         Iteracje: 25
         Oceny funkcji: 30
         Oceny gradientu: 30
[1.00000004 1.0000001  1.00000021 1.00000044 1.00000092]

Algorytm sprzężonych gradientów (Newtona)

Algorytm sprzężonych gradientów Newtona jest zmodyfikowaną metodą Newtona.
Metoda Newtona opiera się na aproksymacji funkcji w lokalnym obszarze wielomianem drugiego stopnia:

SciPy, optymalizacja

gdzie SciPy, optymalizacja jest macierzą drugich pochodnych (macierz Hessego, hessian).
Jeśli hessian jest dodatnio określony, lokalne minimum tej funkcji można znaleźć, przyrównując zerowy gradient do zera w formie kwadratowej. W wyniku uzyskuje się wyrażenie:

SciPy, optymalizacja

Odwrotność hessiana jest obliczana za pomocą metody sprzężonych gradientów. Przykład zastosowania tej metody do minimalizacji funkcji Rosenbrocka przedstawiono poniżej. Aby użyć metody Newton-CG, należy określić funkcję, która oblicza hessian.
Hessian funkcji Rosenbrocka w postaci analitycznej wynosi:

SciPy, optymalizacja

SciPy, optymalizacja

gdzie SciPy, optymalizacja i SciPy, optymalizacja, określają macierz SciPy, optymalizacja.

Pozostałe niena zero elementy macierzy są równe:

SciPy, optymalizacja

SciPy, optymalizacja

SciPy, optymalizacja

SciPy, optymalizacja

Na przykład, w pięciowymiarowej przestrzeni N = 5, macierz Hessego dla funkcji Rosenbrocka ma postać taśmy:

SciPy, optymalizacja

Kod, który oblicza ten hessian wraz z kodem do minimalizacji funkcji Rosenbrocka za pomocą metody sprzężonych gradientów (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)

Optymalizacja zakończona pomyślnie.
         Aktualna wartość funkcji: 0.000000
         Iteracje: 24
         Oceny funkcji: 33
         Oceny gradientu: 56
         Oceny hesjanu: 24
[1.         1.         1.         0.99999999 0.99999999]

Przykład z definicją funkcji iloczynu hesjanu i dowolnego wektora

W rzeczywistych zadaniach obliczenie i przechowywanie całej macierzy Hessego może wymagać znacznych zasobów czasowych i pamięciowych. Nie ma jednak potrzeby definiowania samej macierzy Hessego, ponieważ do procedury minimalizacji potrzebny jest tylko wektor będący wynikiem iloczynu hesjanu z innym dowolnym wektorem. Z komputerowego punktu widzenia znacznie korzystniejsze jest zdefiniowanie funkcji, która od razu zwraca wynik iloczynu hesjanu z dowolnym wektorem.

Rozważmy funkcję hess, która przyjmuje wektor minimalizacji jako pierwszy argument, a dowolny wektor jako drugi argument (obok innych argumentów funkcji minimalizowanej). W naszym przypadku obliczenie iloczynu hesjanu funkcji Rosenbrocka z dowolnym wektorem nie jest szczególnie trudne. Jeżeli p — dowolny wektor, to iloczyn SciPy, optymalizacja ma postać:

SciPy, optymalizacja

Funkcja obliczająca iloczyn hesjanu i dowolnego wektora jest przekazywana jako wartość argumentu hessp funkcji minimize:

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})


Optymalizacja zakończona pomyślnie.
         Aktualna wartość funkcji: 0.000000
         Iteracje: 24
         Oceny funkcji: 33
         Oceny gradientu: 56
         Oceny hesjanu: 66

Algorytm zaufanej strefy (trust region) sprzężonych gradientów (Newtona)

Słaba kondycja macierzy Hessego i błędne kierunki poszukiwań mogą spowodować, że algorytm sprzężonych gradientów Newtona może być nieefektywny. W takich przypadkach preferuje się metodę zaufanej strefy (trust-region) sprzężonych 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)

Optymalizacja zakończona pomyślnie.
         Aktualna wartość funkcji: 0.000000
         Iteracje: 20
         Oceny funkcji: 21
         Oceny gradientu: 20
         Oceny hesjanu: 19
[1. 1. 1. 1. 1.]

Przykład z funkcją iloczynu hesjanu 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)

Optymalizacja zakończona pomyślnie.
         Aktualna wartość funkcji: 0.000000
         Iteracje: 20
         Oceny funkcji: 21
         Oceny gradientu: 20
         Oceny hesjanu: 0
[1. 1. 1. 1. 1.]

Metody Krylowa

Podobnie jak metoda trust-ncg, metody Krylowa są dobrze dostosowane do rozwiązywania dużych problemów, ponieważ wykorzystują jedynie iloczyny macierzowo-wektorowe. Ich istota polega na rozwiązywaniu problemu w przestrzeni zaufania, ograniczonej przyciętym podprzestrzenią Krylowa. W przypadku problemów nieokreślonych lepiej jest stosować tę metodę, ponieważ wykorzystuje mniej nieliniowych iteracji przy mniejszej liczbie iloczynów macierzowo-wektorowych na jeden podproblem, w porównaniu do metody trust-ncg. Ponadto, rozwiązanie kwadratowego podproblemu jest uzyskiwane z większą dokładnością niż w metodzie 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})

Optymalizacja zakończona pomyślnie.
         Aktualna wartość funkcji: 0.000000
         Iteracje: 19
         Oceny funkcji: 20
         Oceny gradientu: 20
         Oceny hesjanu: 18

print(res.x)

    [1. 1. 1. 1. 1.]

Przykład z funkcją iloczynu hesjanu i dowolnego wektora:

res = minimize(rosen, x0, method='trust-krylov',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})

Optymalizacja zakończona pomyślnie.
         Aktualna wartość funkcji: 0.000000
         Iteracje: 19
         Oceny funkcji: 20
         Oceny gradientu: 20
         Oceny hesjanu: 0

print(res.x)

    [1. 1. 1. 1. 1.]

Algorytm przybliżonego rozwiązania w obszarze zaufania

Wszystkie metody (Newton-CG, trust-ncg i trust-krylov) są dobrze dostosowane do rozwiązywania dużych problemów (z tysiącami zmiennych). Ma to związek z tym, że leżący u podstaw algorytm gradientów sprzężonych zakłada przybliżone znajdowanie odwrotnej macierzy hesjanu. Rozwiązanie znajduje się iteracyjnie, bez jawnego rozkładu hesjanu. Ponieważ wymaga się jedynie określenia funkcji dla iloczynu hesjanu i dowolnego wektora, ten algorytm jest szczególnie dobry w pracy z rzadkimi (pasmowymi) macierzami. Zapewnia to niskie koszty pamięci i znaczne oszczędności czasu.

W zadaniach średniej wielkości koszty przechowywania i faktoryzacji macierzy Hessego nie mają kluczowego znaczenia. Oznacza to, że można uzyskać rozwiązanie w mniejszej liczbie iteracji, rozwiązując podzadania obszaru zaufania prawie dokładnie. W tym celu niektóre nieliniowe równania są rozwiązywane iteracyjnie dla każdego kwadratowego podzadania. Takie rozwiązanie wymaga zazwyczaj 3 lub 4 dekompozycji choleskiego macierzy Hessego. W rezultacie metoda zbiega w mniejszej liczbie iteracji i wymaga mniej obliczeń funkcji celu niż inne wdrożone metody obszaru zaufania. Ten algorytm zakłada tylko określenie pełnej macierzy Hessego i nie obsługuje możliwości użycia funkcji iloczynu Hessego i dowolnego wektora.

Przykład minimalizacji funkcji Rosenbrocka:

res = minimize(rosen, x0, method='trust-exact',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})
res.x

Optymalizacja zakończona sukcesem.
         Obecna wartość funkcji: 0.000000
         Iteracje: 13
         Oceny funkcji: 14
         Oceny gradientu: 13
         Oceny Hessego: 14

array([1., 1., 1., 1., 1.])

Na tym, myślę, zakończymy. W następnym artykule postaram się opowiedzieć o najciekawszych aspektach minimalizacji warunkowej, zastosowaniu minimalizacji w rozwiązywaniu problemów aproksymacji, minimalizacji funkcji jednej zmiennej, dowolnych minimalizatorach oraz szukaniu pierwiastków systemu równań za pomocą pakietu scipy.optimize.

Źródło: https://docs.scipy.org/doc/scipy/reference/

Źródło: habr.com

Kup solidny hosting stron z ochroną przed DDoS, serwery VPS VDS 🔥 Kup solidny hosting stron z ochroną przed DDoS, serwery VPS VDS | ProHoster