SciPy (ausgesprochen „Sai Pie“) ist ein mathematisches Anwendungspaket, das auf der Numpy-Python-Erweiterung basiert. Mit SciPy wird Ihre interaktive Python-Sitzung zur gleichen vollständigen Datenwissenschafts- und komplexen System-Prototyping-Umgebung wie MATLAB, IDL, Octave, R-Lab und SciLab. Heute möchte ich kurz darüber sprechen, wie einige bekannte Optimierungsalgorithmen im Paket scipy.optimize verwendet werden. Ausführlichere und aktuellere Hilfe zur Verwendung von Funktionen erhalten Sie jederzeit mit dem Befehl help() oder mit Umschalt+Tab.
Einführung
Um sich und den Lesern das Suchen und Lesen von Primärquellen zu ersparen, finden sich Links zu Methodenbeschreibungen überwiegend auf Wikipedia. In der Regel reichen diese Informationen aus, um die Methoden im Allgemeinen und die Bedingungen für ihre Anwendung zu verstehen. Um die Essenz mathematischer Methoden zu verstehen, folgen Sie den Links zu maßgeblicheren Veröffentlichungen, die Sie am Ende jedes Artikels oder in Ihrer bevorzugten Suchmaschine finden.
Das Modul scipy.optimize umfasst also die Implementierung der folgenden Verfahren:
- Bedingte und unbedingte Minimierung von Skalarfunktionen mehrerer Variablen (Minim) unter Verwendung verschiedener Algorithmen (Nelder-Mead-Simplex, BFGS, Newton-konjugierte Gradienten, и )
- Globale Optimierung (zum Beispiel: , )
- Minimierung von Rückständen (least_squares) und Kurvenanpassungsalgorithmen unter Verwendung nichtlinearer kleinster Quadrate (curve_fit)
- Skalarfunktionen einer Variablen minimieren (minim_scalar) und nach Wurzeln suchen (root_scalar)
- Mehrdimensionale Löser von Gleichungssystemen (Wurzel) unter Verwendung verschiedener Algorithmen (hybrider Powell, oder groß angelegte Methoden wie z ).
In diesem Artikel betrachten wir nur den ersten Punkt dieser gesamten Liste.
Unbedingte Minimierung einer Skalarfunktion mehrerer Variablen
Die Minim-Funktion aus dem Paket scipy.optimize stellt eine allgemeine Schnittstelle zur Lösung bedingter und unbedingter Minimierungsprobleme von Skalarfunktionen mehrerer Variablen bereit. Um zu demonstrieren, wie es funktioniert, benötigen wir eine geeignete Funktion mehrerer Variablen, die wir auf unterschiedliche Weise minimieren.
Für diese Zwecke ist die Rosenbrock-Funktion von N Variablen perfekt, die die Form hat:
Obwohl die Rosenbrock-Funktion und ihre Jacobi- und Hessischen Matrizen (die erste bzw. zweite Ableitung) bereits im Paket scipy.optimize definiert sind, werden wir sie selbst definieren.
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)Zur Verdeutlichung zeichnen wir die Werte der Rosenbrock-Funktion zweier Variablen in 3D.
Zeichnungscode
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()
Im Voraus wissen, dass das Minimum bei 0 liegt Sehen wir uns Beispiele an, wie der Mindestwert der Rosenbrock-Funktion mithilfe verschiedener scipy.optimize-Prozeduren ermittelt werden kann.
Nelder-Mead-Simplex-Methode
Es gebe einen Anfangspunkt x0 im 5-dimensionalen Raum. Finden wir mithilfe des Algorithmus den nächstgelegenen Minimalpunkt der Rosenbrock-Funktion (Der Algorithmus wird als Wert des Methodenparameters 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 339
Function evaluations: 571
[1. 1. 1. 1. 1.]Die Simplex-Methode ist die einfachste Möglichkeit, eine explizit definierte und ziemlich glatte Funktion zu minimieren. Es ist nicht erforderlich, die Ableitungen einer Funktion zu berechnen; es reicht aus, nur ihre Werte anzugeben. Für einfache Minimierungsprobleme ist die Nelder-Mead-Methode eine gute Wahl. Da jedoch keine Gradientenschätzungen verwendet werden, kann es länger dauern, das Minimum zu finden.
Powell-Methode
Ein weiterer Optimierungsalgorithmus, bei dem nur die Funktionswerte berechnet werden, ist . Um es zu verwenden, müssen Sie method = 'powell' in der Minim-Funktion festlegen.
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)-Algorithmus
Um eine schnellere Konvergenz zu einer Lösung zu erreichen, ist das Verfahren verwendet den Gradienten der Zielfunktion. Der Gradient kann als Funktion angegeben oder anhand von Differenzen erster Ordnung berechnet werden. In jedem Fall erfordert die BFGS-Methode typischerweise weniger Funktionsaufrufe als die Simplex-Methode.
Finden wir die Ableitung der Rosenbrock-Funktion in analytischer Form:
Dieser Ausdruck gilt für die Ableitungen aller Variablen außer der ersten und der letzten, die wie folgt definiert sind:
Schauen 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 jac-Parameters der Minimfunktion angegeben, wie unten gezeigt.
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]Konjugierter Gradientenalgorithmus (Newton)
Algorithm ist eine modifizierte Newton-Methode.
Newtons Methode basiert auf der Approximation einer Funktion in einem lokalen Bereich durch ein Polynom zweiten Grades:
wo ist die Matrix der zweiten Ableitungen (Hessische Matrix, Hessisch).
Wenn die Hesse-Funktion positiv definit ist, kann das lokale Minimum dieser Funktion ermittelt werden, indem der Nullgradient der quadratischen Form mit Null gleichgesetzt wird. Das Ergebnis wird der Ausdruck sein:
Der inverse Hessesche Wert wird mithilfe der konjugierten Gradientenmethode berechnet. Nachfolgend finden Sie ein Beispiel für die Verwendung dieser Methode zur Minimierung der Rosenbrock-Funktion. Um die Newton-CG-Methode verwenden zu können, müssen Sie eine Funktion angeben, die den Hesse-Wert berechnet.
Die Hesse-Funktion der Rosenbrock-Funktion in analytischer Form ist gleich:
wo и , definieren Sie die Matrix .
Die verbleibenden Nicht-Null-Elemente der Matrix sind gleich:
In einem fünfdimensionalen Raum N = 5 hat die Hesse-Matrix für die Rosenbrock-Funktion beispielsweise die Form eines Bandes:
Code, der diese Hesse-Funktion zusammen mit Code zur Minimierung der Rosenbrock-Funktion mithilfe der konjugierten Gradientenmethode (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)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]Ein Beispiel mit der Definition der Produktfunktion des Hesseschen und eines beliebigen Vektors
Bei realen Problemen kann die Berechnung und Speicherung der gesamten Hesse-Matrix erhebliche Zeit- und Speicherressourcen erfordern. In diesem Fall ist es eigentlich nicht erforderlich, die Hesse-Matrix selbst anzugeben, weil Das Minimierungsverfahren erfordert nur einen Vektor, der dem Produkt des Hesseschen mit einem anderen beliebigen Vektor entspricht. Aus rechnerischer Sicht ist es daher viel vorzuziehen, sofort eine Funktion zu definieren, die das Ergebnis des Produkts der Hesse-Funktion mit einem beliebigen Vektor zurückgibt.
Betrachten Sie die Hess-Funktion, die den Minimierungsvektor als erstes Argument und einen beliebigen Vektor als zweites Argument (zusammen mit anderen Argumenten der zu minimierenden Funktion) verwendet. In unserem Fall ist die Berechnung des Produkts der Hesse-Funktion der Rosenbrock-Funktion mit einem beliebigen Vektor nicht sehr schwierig. Wenn p ein beliebiger Vektor ist, dann das Produkt hat die Form:
Die Funktion, die das Produkt aus dem Hesse-Wert und einem beliebigen Vektor berechnet, wird als Wert des Hessp-Arguments an die Minimierungsfunktion ü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})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 66Konjugierter Gradienten-Trust-Region-Algorithmus (Newton)
Eine schlechte Konditionierung der Hesse-Matrix und falsche Suchrichtungen können dazu führen, dass der konjugierte Gradientenalgorithmus von Newton unwirksam wird. In solchen Fällen wird bevorzugt (Vertrauensbereich) konjugierte Newton-Gradienten.
Beispiel mit der 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)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 19
[1. 1. 1. 1. 1.]Beispiel mit der Produktfunktion des Hesseschen 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)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 vom Krylov-Typ
Wie die Trust-NCG-Methode eignen sich auch Krylov-Methoden gut zur Lösung großräumiger Probleme, da sie nur Matrix-Vektor-Produkte verwenden. Ihr Kern besteht darin, ein Problem in einem Vertrauensbereich zu lösen, der durch einen verkürzten Krylov-Unterraum begrenzt ist. Für unsichere Probleme ist es besser, diese Methode zu verwenden, da sie im Vergleich zur Trust-NCG-Methode aufgrund der geringeren Anzahl von Matrix-Vektor-Produkten pro Teilproblem eine geringere Anzahl nichtlinearer Iterationen verwendet. Darüber hinaus wird die Lösung des quadratischen Teilproblems genauer gefunden als mit der Trust-NCG-Methode.
Beispiel mit der Definition der Hesse-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.]
Beispiel mit der Produktfunktion des Hesseschen und eines beliebigen Vektors:
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.]
Algorithmus zur Näherungslösung im Vertrauensbereich
Alle Methoden (Newton-CG, Trust-NCG und Trust-Krylov) eignen sich gut zur Lösung großräumiger Probleme (mit Tausenden von Variablen). Dies liegt daran, dass der zugrunde liegende konjugierte Gradientenalgorithmus eine näherungsweise Bestimmung der inversen Hesse-Matrix impliziert. Die Lösung wird iterativ gefunden, ohne explizite Erweiterung des Hessischen. Da Sie nur eine Funktion für das Produkt einer Hesse-Funktion und eines beliebigen Vektors definieren müssen, eignet sich dieser Algorithmus besonders gut für die Arbeit mit dünn besetzten Matrizen (Banddiagonalen). Dies sorgt für geringe Speicherkosten und eine erhebliche Zeitersparnis.
Bei mittelgroßen Problemen sind die Kosten für die Speicherung und Faktorisierung des Hessischen nicht kritisch. Dies bedeutet, dass es möglich ist, in weniger Iterationen eine Lösung zu erhalten, die die Teilprobleme der Vertrauensregion nahezu exakt löst. Dazu werden für jedes quadratische Teilproblem einige nichtlineare Gleichungen iterativ gelöst. Eine solche Lösung erfordert normalerweise 3 oder 4 Cholesky-Zerlegungen der Hesse-Matrix. Dadurch konvergiert die Methode in weniger Iterationen und erfordert weniger Zielfunktionsberechnungen als andere implementierte Konfidenzbereichsmethoden. Dieser Algorithmus umfasst nur die Bestimmung der vollständigen Hesse-Matrix und unterstützt nicht die Möglichkeit, die Produktfunktion der Hesse-Matrix und eines beliebigen Vektors zu verwenden.
Beispiel mit Minimierung der Rosenbrock-Funktion:
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.])Wir werden dort wahrscheinlich aufhören. Im nächsten Artikel werde ich versuchen, das Interessanteste über die bedingte Minimierung, die Anwendung der Minimierung bei der Lösung von Approximationsproblemen, die Minimierung einer Funktion einer Variablen, beliebige Minimierer und das Finden der Wurzeln eines Gleichungssystems mithilfe von scipy.optimize zu erzählen Paket.
Source:
Source: habr.com