SciPy (ausgesprochen wie "sai pai") ist ein Paket für angewandte mathematische Verfahren, das auf der Numpy-Erweiterung von Python basiert. Mit SciPy wird eine interaktive Python-Sitzung zu einer vollwertigen Umgebung für Datenverarbeitung und Prototypenerstellung komplexer Systeme, ähnlich wie MATLAB, IDL, Octave, R-Lab und SciLab. Heute möchte ich kurz erläutern, wie einige bekannte Optimierungsalgorithmen aus dem Paket scipy.optimize angewendet werden können. Detaillierte und aktuelle Informationen zur Nutzung der Funktionen erhalten Sie jederzeit mit dem Befehl help() oder durch Drücken von Shift+Tab.
Einführung
Um sowohl mir als auch den Lesern die Suche und das Lesen von Primärquellen zu ersparen, beziehen sich die Links zu den Methodenerklärungen hauptsächlich auf Wikipedia. In der Regel reicht diese Information aus, um die Methoden im Allgemeinen zu verstehen und die Bedingungen ihrer Anwendung zu erfahren. Für ein tiefgehendes Verständnis der mathematischen Methoden folgen wir den Links zu anerkannten Publikationen, die am Ende jedes Artikels oder über die gewählte Suchmaschine zu finden sind.
Das Modul scipy.optimize umfasst die Implementierung der folgenden Verfahren:
- Bedingte und unbedingte Minimierung von skalarer Funktionen mit mehreren Variablen (minim) unter Verwendung verschiedener Algorithmen (Nelder-Mead-Simplex, BFGS, Newtons Methode mit konjugierten Gradienten, und )
- Globale Optimierung (zum Beispiel: , )
- Minimierung der Residuen (least_squares) und Algorithmen zur Anpassung von Kurven mit nichtlinearem MSE (curve_fit)
- Minimierung skalarer Funktionen mit einer Variablen (minim_scalar) und Wurzelberechnung (root_scalar)
- Mehrdimensionale Lösungsverfahren für Gleichungssysteme (root) unter Einsatz verschiedener Algorithmen (hybride Powell-Methode, oder groß angelegte Methoden wie ).
In diesem Artikel betrachten wir nur den ersten Punkt aus dieser Liste.
Unbedingte Minimierung einer skalarer Funktion mit mehreren Variablen
Die Funktion minim aus dem Paket scipy.optimize bietet eine allgemeine Schnittstelle zur Lösung von Aufgaben der bedingten und unbedingten Minimierung von skalarer Funktionen mit mehreren Variablen. Um ihre Arbeitsweise zu demonstrieren, benötigen wir eine geeignete Funktion mit mehreren Variablen, die wir unterschiedlich minimieren werden.
Für diese Zwecke eignet sich hervorragend die Rosenbrock-Funktion von N Variablen, die wie folgt aussieht:
Obwohl die Rosenbrock-Funktion sowie ihre Jakobimatrix und Hessian (erste bzw. 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)Zur Veranschaulichung werden wir die Werte der Rosenbrock-Funktion von zwei Variablen in 3D zeichnen.
Der Code zur Visualisierung
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# 3D-Plot einrichten
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Betrachtungswinkel festlegen
ax.view_init(45, 30)
# Daten für den Plot 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 wissen, dass das Minimum bei 0 liegt, 
, betrachten wir Beispiele dafür, wie man den minimalen Wert der Rosenbrock-Funktion mit verschiedenen Verfahren von scipy.optimize bestimmen kann.
Nelder-Mead-Simplexverfahren
Angenommen, es gibt einen Startpunkt x0 im 5-dimensionalen Raum. Wir werden den Punkt finden, der am nächsten zum Minimum der Rosenbrock-Funktion liegt, mithilfe des Algorithmus. (Algorithmus als Parameterwert '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)Optimierung erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 339
Funktionsauswertungen: 571
[1. 1. 1. 1. 1.]Die Simplex-Methode ist der einfachste Weg, eine eindeutig definierte und relativ glatte Funktion zu minimieren. Sie erfordert keine Berechnung der Ableitungen der Funktion, es genügt, ihre Werte anzugeben. Der Nelder-Mead-Algorithmus ist eine gute Wahl für einfache Minimierungsprobleme. Da er jedoch keine Gradienteninformationen nutzt, kann die Suche nach dem Minimum länger dauern.
Powell-Methode
Ein weiterer Optimierungsalgorithmus, bei dem nur Funktionswerte berechnet werden, ist . Um sie zu verwenden, setzen Sie method = 'powell' in der Funktion 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)Optimierung 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 zu erreichen, wird das Verfahren verwendet den Gradienten der Zielfunktion. Der Gradient kann entweder als Funktion angegeben oder durch erste Differenzen berechnet werden. In jedem Fall benötigt das BFGS-Verfahren normalerweise weniger Funktionsaufrufe als das Simplex-Verfahren.
Lassen Sie uns die Ableitung der Rosenbrock-Funktion analytisch finden:
Dieser Ausdruck gilt für die Ableitungen aller Variablen, mit Ausnahme der ersten und letzten, die definiert sind als:
Sehen wir uns die Python-Funktion an, die diesen Gradienten 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 Gradientenberechnungsfunktion wird als Wert des Parameters jac der Funktion minim angegeben, wie unten dargestellt.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Optimierung wurde erfolgreich beendet.
Aktueller Funktionswert: 0,000000
Iterationen: 25
Funktionsauswertungen: 30
Gradientenauswertungen: 30
[1,00000004 1,0000001 1,00000021 1,00000044 1,00000092]Algorithmus der konjugierten Gradient (Newton)
Algorithmus ist eine modifizierte Methode von Newton.
Die Newton-Methode basiert auf der Approximierung einer Funktion in einem lokalen Bereich mit einem Polynom zweiten Grades:
wobei 
ist die Matrix der zweiten Ableitungen (Hessian-Matrix, Hessian).
Wenn der Hessian positiv definit ist, kann das lokale Minimum dieser Funktion gefunden werden, indem der Gradient der quadratischen Form auf null gesetzt wird. Es resultiert in folgendem Ausdruck:
Der inverse Hessian wird mit dem Verfahren der konjugierten Gradienten berechnet. Ein Beispiel für die Anwendung dieser Methode zur Minimierung der Rosenbrock-Funktion wird unten gezeigt. Um die Newton-CG-Methode zu verwenden, muss eine Funktion definiert werden, die den Hessian berechnet.
Der Hessian der Rosenbrock-Funktion im analytischen Format lautet:
wobei 
und 
, matrix definieren 
.
Die restlichen von Null verschiedenen Elemente der Matrix sind:
Zum Beispiel hat die Hessische Matrix für die Rosenbrock-Funktion im fünf-dimensionalen Raum (N = 5) eine bandartige Form:
Der Code, der diesen Hessian zusammen mit dem Code zur Minimierung der Rosenbrock-Funktion mittels des Verfahrens der konjugierten Gradienten (Newton) berechnet:
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)Optimierung erfolgreich abgeschlossen.
Aktueller Funktionswert: 0.000000
Iterationen: 24
Funktionsbewertungen: 33
Gradientbewertungen: 56
Hessianbewertungen: 24
[1. 1. 1. 0.99999999 0.99999999]Beispiel zur Bestimmung der Funktion des Hessianprodukts und eines beliebigen Vektors
In realen Anwendungen kann die Berechnung und Speicherung der gesamten Hessianmatrix erhebliche zeitliche und speichermäßige Ressourcen erfordern. Tatsächlich ist es nicht notwendig, die Hessianmatrix selbst festzulegen, da für das Minimierungsverfahren nur ein Vektor benötigt wird, der dem Produkt des Hessians mit einem anderen beliebigen Vektor entspricht. Daher ist es aus rechnerischer Sicht viel vorteilhafter, direkt eine Funktion zu definieren, die das Ergebnis des Produkts des Hessians mit einem beliebigen Vektor zurückgibt.
Betrachten wir die Funktion hess, die einen Minimierungsvektor als erstes Argument und einen beliebigen Vektor als zweites Argument (neben anderen Argumenten der zu minimierenden Funktion) akzeptiert. In unserem Fall ist es nicht sehr kompliziert, das Produkt des Hessematrix der Rosenbrock-Funktion mit einem beliebigen Vektor zu berechnen. Wenn p es sich um einen beliebigen Vektor handelt, dann hat das Produkt 
die folgende Form:
Die Funktion zur Berechnung des Produkts der Hessischen Matrix und eines beliebigen Vektors wird als Argumentwert 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
Gradientenauswertungen: 56
Hessian-Auswertungen: 66Der Algorithmus der Vertrauensregion (trust region) der konjugierten Gradienten (Newton)
Eine ungünstige Bedingung der Hessian-Matrix und inkorrekte Suchrichtungen können dazu führen, dass der Newtonsche konjugierte Gradientenalgorithmus ineffektiv ist. In solchen Fällen wird bevorzugt (Trust-Region) der konjugierten Gradienten von Newton.
Beispiel zur Bestimmung der Hessian-Matrix:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimierung erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 20
Funktionsauswertungen: 21
Gradientenbewertungen: 20
Hessianbewertungen: 19
[1. 1. 1. 1. 1.]Beispiel mit der Funktion des Produkts der Hessian-Matrix und eines beliebigen Vektors:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimierung erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 20
Funktionsauswertungen: 21
Gradientenbewertungen: 20
Hessianbewertungen: 0
[1. 1. 1. 1. 1.]Krylov-Methoden
Ähnlich wie beim trust-ncg-Verfahren eignen sich die Krylov-Methoden gut zur Lösung von großangelegten Aufgaben, da sie ausschließlich matrix-vektor-Produkte verwenden. Der Grundgedanke besteht darin, das Problem in einem vertrauenswürdigen Bereich zu lösen, der durch einen verkürzten Krylov-Unterraum begrenzt ist. Für nichtdefinierte Probleme ist es besser, diese Methode anzuwenden, da sie weniger nichtlineare Iterationen erfordert und deshalb weniger matrix-vektor-Produkte für jede Teilaufgabe benötigt als die Methode trust-ncg. Darüber hinaus wird die quadratische Teilaufgabe präziser gelöst als bei der Methode trust-ncg.
Beispiel zur Bestimmung der Hessian-Matrix:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Optimierung erfolgreich beendet.
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 Produkts der Hessian-Matrix und eines beliebigen Vektors:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Optimierung erfolgreich beendet.
Aktueller Funktionswert: 0.000000
Iterationen: 19
Funktionsauswertungen: 20
Gradientenauswertungen: 20
Hessianauswertungen: 0
print(res.x)
[1. 1. 1. 1. 1.]
Algorithmus zur approximativen Lösung im Vertrauensbereich
Alle Methoden (Newton-CG, trust-ncg und trust-krylov) eignen sich gut zur Lösung groß angelegter Aufgaben (mit Tausenden von Variablen). Dies liegt daran, dass der zugrunde liegende Algorithmus der konjugierten Gradienten eine approximative Berechnung der inversen Hesse-Matrix voraussetzt. Die Lösung wird iterativ gefunden, ohne dass eine explizite Zerlegung des Hessian durchgeführt wird. Da nur die Funktion für das Produkt aus dem Hessian und einem beliebigen Vektor bestimmt werden muss, ist dieser Algorithmus besonders geeignet für die Arbeit mit dünn besetzten (band-diagonalen) Matrizen. Dies führt zu geringen Speicheranforderungen und einer erheblichen Zeitersparnis.
Bei mittelgroßen Aufgaben spielen die Kosten für die Speicherung und die Faktorisierung der Hesse-Matrix keine entscheidende Rolle. Das bedeutet, dass eine Lösung in weniger Iterationen erzielt werden kann, indem die Teilaufgaben im Vertrauensbereich nahezu genau gelöst werden. Zu diesem Zweck werden einige nichtlineare Gleichungen iterativ für jede quadratische Teilaufgabe gelöst. Diese Lösung erfordert in der Regel 3 oder 4 Zerlegungen der Hesse-Matrix durch Cholesky. Infolgedessen konvergiert die Methode in weniger Iterationen und benötigt weniger Funktionsbewertungen als andere umgesetzte Methoden des Vertrauensbereichs. Dieser Algorithmus setzt lediglich die Bestimmung der vollständigen Hesse-Matrix voraus und unterstützt nicht die Möglichkeit, eine Funktion des Hesses 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
Funktionsbewertungen: 14
Gradientenbewertungen: 13
Hesse-Bewertungen: 14
array([1., 1., 1., 1., 1.])Das war's für jetzt. In meinem nächsten Artikel werde ich versuchen, die interessantesten Aspekte der bedingten Minimierung, der Anwendung der Minimierung zur Lösung von Approximationsproblemen, der Minimierung einer Funktion mit einer Variablen, der allgemeinen Minimierer und der Wurzelbestimmung in Gleichungssystemen mit dem Paket scipy.optimize zu erläutern.
Quelle:
Quelle: habr.com
