SciPy (ausgesprochen wie „sai pai“) ist ein Paket für angewandte mathematische Verfahren, das auf der Erweiterung von Numpy Python basiert. Mit SciPy wird eine interaktive Python-Sitzung zu einer vollständigen Umgebung für die Datenverarbeitung und Prototypenerstellung komplexer Systeme, ähnlich wie MATLAB, IDL, Octave, R-Lab und SciLab. Heute möchte ich kurz darlegen, wie einige bekannte Optimierungsalgorithmen im Paket scipy.optimize angewendet werden sollten. Eine detaillierte und aktuelle Dokumentation zur Verwendung der Funktionen kann jederzeit mit dem Befehl help() oder mit Shift+Tab abgerufen werden.
Einführung
Um mich selbst und die Leser von der Suche nach und dem Lesen der ursprünglichen Quellen zu befreien, werden die Verweise auf die Methodenbeschreibungen hauptsächlich auf Wikipedia sein. In der Regel reicht diese Information aus, um die Methoden in groben Zügen und die Bedingungen ihrer Anwendung zu verstehen. Um das Wesen der mathematischen Methoden zu begreifen, gehen wir über die Links zu autoritäreren Veröffentlichungen, die am Ende jedes Artikels oder in Ihrer bevorzugten Suchmaschine zu finden sind.
Das Modul scipy.optimize umfasst die Implementierung der folgenden Verfahren:
- Bedingte und unbelegte Minimierung von skalaren Funktionen mehrerer Variablen (minim) mithilfe verschiedener Algorithmen (Nelder-Mead-Simplex, BFGS, Newtons konjugierte Gradienten, und )
- Globale Optimierung (zum Beispiel: , )
- Minimierung von Residuen (least_squares) und Algorithmen zur Anpassung von Kurven über nichtlineare MGNK (curve_fit)
- Minimierung von skalaren Funktionen einer Variablen (minim_scalar) und Wurzelverwendung (root_scalar)
- Mehrdimensionale Lösungsverfahren für Gleichungssysteme (root) unter Verwendung verschiedener Algorithmen (hybrider Powell, oder großmaßstäbliche Methoden wie ).
In diesem Artikel werden wir nur den ersten Punkt aus dieser Liste betrachten.
Unbedingte Minimierung einer skalaren Funktion mehrerer Variablen
Die Funktion minim aus dem Paket scipy.optimize bietet eine allgemeine Schnittstelle zur Lösung von Aufgaben zur bedingten und unbedingten Minimierung skalaren Funktionen mehrerer Variablen. Um ihre Funktionsweise zu demonstrieren, benötigen wir eine geeignete Funktion mehrerer Variablen, die wir unterschiedlich minimieren werden.
Für diese Zwecke eignet sich die Rosenbrock-Funktion von N Variablen, die wie folgt aussieht:
Obwohl die Rosenbrock-Funktion und ihre Jacobi- und Hessischen Matrizen (erste und zweite Ableitung) bereits im Paket scipy.optimize definiert sind, werden wir sie selbst definieren.
import numpy as np
def rosen(x):
"""Die Rosenbrock-Funktion"""
return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)Um die Werte der Rosenbrock-Funktion mit zwei Variablen zu veranschaulichen, zeichnen wir diese in 3D auf.
Code zur Darstellung
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# 3D-Diagramm einrichten
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Sichtwinkel festlegen
ax.view_init(45, 30)
# Daten für das Diagramm erstellen
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]))
# Oberfläche zeichnen
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
Da wir im Voraus wissen, dass das Minimum bei 0 liegt, 
, betrachten wir Beispiele, wie man den minimalen Wert der Rosenbrock-Funktion mit verschiedenen scipy.optimize-Verfahren bestimmen kann.
Nelder-Mead-Simplex-Methode
Angenommen, wir haben einen Startpunkt x0 im 5-dimensionalen Raum. Finden wir den nächstgelegenen Punkt zum Minimum der Rosenbrock-Funktion mit dem Algorithmus (der Algorithmus wird als Wert des Parameters method angegeben):
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)Die Optimierung wurde erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 339
Funktionsauswertungen: 571
[1. 1. 1. 1. 1.]Die Simplex-Methode ist der einfachste Weg, um eine explizit definierte und relativ glatte Funktion zu minimieren. Sie erfordert keine Berechnung der Ableitungen der Funktion, es genügt, nur ihre Werte festzulegen. Die Nelder-Mead-Methode ist eine gute Wahl für einfache Minimierungsprobleme. Da sie jedoch keine Gradientenbewertungen verwendet, kann es länger dauern, ein Minimum zu finden.
Powell-Methode
Ein weiterer Optimierungsalgorithmus, der nur Funktionswerte berechnet, ist die . Um sie zu verwenden, muss method = 'powell' in der Funktion minim gesetzt werden.
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)Die Optimierung wurde erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 19
Funktionsauswertungen: 1622
[1. 1. 1. 1. 1.]Broyden-Fletcher-Goldfarb-Shanno-Algorithmus (BFGS)
Um eine schnellere Konvergenz zur Lösung zu erreichen, verwendet das Verfahren verwendet den Gradienten der Zielfunktion. Der Gradient kann als Funktion angegeben oder mit Hilfe von Differenzen erster Ordnung berechnet werden. In jedem Fall erfordert das BFGS-Verfahren normalerweise weniger Funktionsaufrufe als das Simplexverfahren.
Lassen Sie uns die Ableitung der Rosenbrock-Funktion analytisch bestimmen:
Dieser Ausdruck gilt für die Ableitungen aller Variablen, außer der ersten und letzten, die wie folgt definiert sind:
Sehen wir uns die Python-Funktion an, die diesen Gradient berechnet:
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 derDie Funktion zur Berechnung des Gradienten wird als Wert des Parameters jac der Funktion minim angegeben, wie unten gezeigt.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Optimierung erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 25
Funktionsbewertungen: 30
Gradientbewertungen: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Der Algorithmus der konjugierten Gradienten (Newton)
Algorithmus ist eine modifizierte Methode von Newton.
Die Newton-Methode basiert auf der Annäherung der Funktion in einem lokalen Bereich durch ein Polynom zweiten Grades:
wo 
ist die Matrix der zweiten Ableitungen (Hessische Matrix, Hesse).
Wenn die Hesse positiv definit ist, kann das lokale Minimum dieser Funktion gefunden werden, indem der Gradienten null für die quadratische Form gleich null gesetzt wird. Das ergibt den Ausdruck:
Die umgekehrte Hesse wird mit dem konjugierten Gradientverfahren berechnet. Ein Beispiel für die Anwendung dieser Methode zur Minimierung der Rosenbrock-Funktion ist unten aufgeführt. Um die Newton-CG-Methode zu verwenden, muss eine Funktion angegeben werden, die die Hesse berechnet.
Die Hesse der Rosenbrock-Funktion ist analytisch gleich:
wo 
und 
, definiert die Matrix 
.
Die anderen von null verschiedenen Elemente der Matrix sind gleich:
Zum Beispiel hat im fünfdimensionalen Raum N = 5 die Hessische Matrix für die Rosenbrock-Funktion eine bandförmige Form:
Code, der diese Hesse berechnet, zusammen mit dem Code zur Minimierung der Rosenbrock-Funktion mit Hilfe des konjugierten Gradientverfahrens (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)Die Optimierung wurde erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 24
Funktionsauswertungen: 33
Gradientenbewertungen: 56
Hesse-Bewertungen: 24
[1. 1. 1. 0.99999999 0.99999999]Beispiel zur Definition der Hesse-Matrix und eines beliebigen Vektors
In praktischen Anwendungen kann die Berechnung und Speicherung der gesamten Hesse-Matrix erhebliche Ressourcen an Zeit und Speicherplatz erfordern. Dabei ist es tatsächlich nicht notwendig, die Hesse-Matrix selbst anzugeben, da für das Minimierungsverfahren nur ein Vektor benötigt wird, der dem Produkt der Hesse-Matrix mit einem anderen beliebigen Vektor entspricht. Aus rechentechnischer Sicht ist es daher viel vorteilhafter, sofort eine Funktion zu definieren, die das Ergebnis des Produkts der Hesse-Matrix mit einem beliebigen Vektor zurückgibt.
Betrachten wir die Funktion hess, die den Minimierungsvektor als ersten Argument und den beliebigen Vektor als zweiten Argument (neben anderen Argumenten der zu minimierenden Funktion) akzeptiert. In unserem Fall ist es nicht sehr schwierig, das Produkt der Hesse-Matrix der Rosenbrock-Funktion mit einem beliebigen Vektor zu berechnen. Wenn p – ein beliebiger Vektor ist, dann hat das Produkt 
die folgende Form:
Die Funktion, die das Produkt der Hesse-Matrix und eines beliebigen Vektors berechnet, wird als Wert des Arguments hessp der Funktion minimize übergeben:
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})
Die Optimierung wurde erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 24
Funktionsauswertungen: 33
Gradientenbewertungen: 56
Hesse-Bewertungen: 66Der Trust-Region-Algorithmus der konjugierten Gradienten (Newton-Methode)
Eine schlechtere Konditionierung der Hesse-Matrix und fehlerhafte Suchrichtungen können dazu führen, dass der Newtonsche konjugierte Gradient Algorithmus ineffektiv ist. In solchen Fällen wird der (trust-region) der konjugierten Gradienten von Newton bevorzugt.
Beispiel zur Definition der Hesse-Matrix:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Die Optimierung wurde erfolgreich abgeschlossen.
Aktueller Funktionswert: 0.000000
Iterationen: 20
Funktionsauswertungen: 21
Gradientenauswertungen: 20
Hessianauswertungen: 19
[1. 1. 1. 1. 1.]Beispiel mit der Funktion des Produktes von Hessian und einem beliebigen Vektor:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Die Optimierung wurde erfolgreich abgeschlossen.
Aktueller Funktionswert: 0.000000
Iterationen: 20
Funktionsauswertungen: 21
Gradientenauswertungen: 20
Hessianauswertungen: 0
print(res.x)
[1. 1. 1. 1. 1.]Krylov-Methoden
Ähnlich wie die Methode trust-ncg eignen sich Krylov-Methoden gut zur Lösung großskaliger Aufgaben, da sie nur Matrix-Vektor-Produkte verwenden. Ihr Wesenszug besteht darin, das Problem in einem Vertrauensbereich zu lösen, der durch den abgeschnittenen Krylov-Unterraum begrenzt ist. Für unbestimmte Probleme ist diese Methode besser geeignet, da sie weniger nichtlineare Iterationen benötigt, bedingt durch die geringere Anzahl an Matrix-Vektor-Produkten pro Teilaufgabe im Vergleich zur Methode trust-ncg. Darüber hinaus wird die Lösung der quadratischen Teilaufgabe genauer gefunden als mit der Methode trust-ncg.
Beispiel zur Definition der Hesse-Matrix:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Die Optimierung wurde erfolgreich abgeschlossen.
Aktueller Funktionswert: 0.000000
Iterationen: 19
Funktionsauswertungen: 20
Gradientenauswertungen: 20
Hessianauswertungen: 18
print(res.x)
[1. 1. 1. 1. 1.]
Beispiel mit der Funktion des Produktes von Hessian und einem beliebigen Vektor:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Die Optimierung wurde erfolgreich abgeschlossen.
Aktueller Funktionswert: 0.000000
Iterationen: 19
Funktionsauswertungen: 20
Gradientenauswertungen: 20
Hessianauswertungen: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algorithmus für die approximative Lösung im Vertrauensbereich
Alle Methoden (Newton-CG, trust-ncg und trust-krylov) eignen sich gut für die Lösung großskaliger Probleme (mit Tausenden von Variablen). Dies liegt daran, dass dem zugrunde liegenden Algorithmus der conjugierten Gradienten die approximative Bestimmung der inversen Hessian-Matrix zugrunde liegt. Die Lösung erfolgt iterativ, ohne dass die Hessian explizit zerlegt wird. Da nur die Funktion für das Produkt von Hessian und einem beliebigen Vektor bestimmt werden muss, ist dieser Algorithmus besonders gut für die Arbeit mit spärlichen (Banden-diagonalen) Matrizen geeignet. Dies sorgt für niedrige Speicherkosten und erhebliche Zeitersparnis.
Bei mittleren Aufgaben sind die Kosten für die Speicherung und Faktorisierung der Hesse-Matrix nicht entscheidend. Das bedeutet, dass eine Lösung mit weniger Iterationen erzielt werden kann, indem die Teilprobleme des Vertrauensbereichs fast genau gelöst werden. Dazu werden einige nichtlineare Gleichungen iterativ für jedes quadratische Teilproblem gelöst. Solch eine Lösung erfordert normalerweise 3 oder 4 Cholesky-Zerlegungen der Hesse-Matrix. In der Folge konvergiert die Methode mit weniger Iterationen und benötigt weniger Berechnungen der Ziel-Funktion als andere implementierte Vertrauensbereichsverfahren. Dieser Algorithmus setzt lediglich die Bestimmung der vollständigen Hesse-Matrix voraus und unterstützt nicht die Möglichkeit, die Funktion des Hesse-Produkts und eines beliebigen Vektors zu verwenden.
Beispiel zur Minimierung der Rosenbrock-Funktion:
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimierung erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 13
Funktionsauswertungen: 14
Gradientenauswertungen: 13
Hessianauswertungen: 14
array([1., 1., 1., 1., 1.])Damit möchte ich wohl schließen. Im nächsten Artikel werde ich versuchen, das Interessanteste über bedingte Minimierung, den Einsatz von Minimierung zur Lösung von Approximationsaufgaben, die Minimierung einer Funktion einer Variablen, beliebige Minimierer und die Wurzelsuche von Gleichungssystemen mithilfe des Pakets scipy.optimize zu erzählen.
Quelle:
Quelle: habr.com
