🙌 🚬 🏣 SciPy, Optimierung 🖕🏽 ❄️ 👄

SciPy (ausgesprochen sai pie) ist ein mathematisches Anwendungspaket, das auf der Erweiterung Numpy Python basiert. Mit SciPy wird aus einer interaktiven Python-Sitzung dieselbe komplexe Datenverarbeitungs- und Prototyping-Umgebung für komplexe Systeme 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. Eine detailliertere und aktuellere Hilfe zur Verwendung von Funktionen erhalten Sie immer mit dem Befehl help () oder mit Umschalt + Tab.

Einführung

Um sich und den Lesern das Durchsuchen und Lesen der Quelle zu ersparen, werden Links zu Methodenbeschreibungen hauptsächlich auf Wikipedia verweist. Diese Informationen reichen in der Regel aus, um die Methoden allgemein und die Bedingungen für ihre Anwendung zu verstehen. Um die Essenz mathematischer Methoden zu verstehen, folgen wir den Links zu maßgeblicheren Veröffentlichungen, die am Ende jedes Artikels oder in Ihrer bevorzugten Suchmaschine zu finden sind.

Das Modul scipy.optimize umfasst also die Implementierung der folgenden Verfahren:

Bedingte und bedingungslose Minimierung der Skalarfunktionen mehrerer Variablen (minim) unter Verwendung verschiedener Algorithmen (Nelder-Mead-Simplex, BFGS, konjugierte Newton-Gradienten, COBYLA und SLSQP )
Globale Optimierung (zB: basinhopping , diff_evolution )
Minimierung von Residuen der kleinsten Quadrate (Least_Squares) und Algorithmen zum Anpassen von Kurven an nichtlineare Least Squares (Curve_fit)
Minimierung der Skalarfunktionen einer Variablen (minim_scalar) und Auffinden von Wurzeln (root_scalar)
Mehrdimensionale Löser des Gleichungssystems (Wurzel) unter Verwendung verschiedener Algorithmen (Hybrid-Powell, Levenberg-Marquardt oder groß angelegte Methoden wie Newton-Krylov ).

In diesem Artikel wird nur das erste Element aus dieser gesamten Liste betrachtet.

Bedingungslose Minimierung der Skalarfunktion mehrerer Variablen

Die Minim-Funktion aus dem Paket scipy.optimize bietet eine gemeinsame Schnittstelle zur Lösung der Probleme der bedingten und bedingungslosen Minimierung von Skalarfunktionen mehrerer Variablen. Um seine Arbeit zu demonstrieren, 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:

$f \ left (\ mathbf {x} \ right) = \ sum_ {i = 1} ^ {N-1} [100 \ left (x_ {i + 1} -x_ {i} ^ {2} \ right) ^ {2} + \ left (1-x_ {i} \ right) ^ {2}]$

Obwohl die Rosenbrock-Funktion und ihre Jacobi- und Hessischen Matrizen (die erste bzw. zweite Ableitung) bereits im Paket scipy.optimize definiert sind, definieren wir sie selbst.

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 in 3D die Werte der Rosenbrock-Funktion zweier Variablen.

Code zum Rendern

 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 0 für ist $x_i = 1$ Betrachten Sie Beispiele für die Bestimmung des Mindestwerts der Rosenbrock-Funktion mithilfe verschiedener scipy.optimize-Verfahren.

Die Nelder-Mead-Simplex-Methode (Nelder-Mead)

Es sei ein Anfangspunkt x0 im 5-dimensionalen Raum. Ermitteln Sie den nächstgelegenen Mindestpunkt der Rosenbrock-Funktion mithilfe des Nelder-Mead-Simplex- Algorithmus (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 der einfachste Weg, um eine klar definierte und ziemlich reibungslose Funktion zu minimieren. Es ist nicht erforderlich, Ableitungen einer Funktion zu berechnen, sondern nur deren Werte anzugeben. Die Nelder-Mead-Methode ist eine gute Wahl für einfache Minimierungsprobleme. Da jedoch keine Gradientenschätzungen verwendet werden, kann es länger dauern, das Minimum zu finden.

Powell-Methode

Ein weiterer Optimierungsalgorithmus, bei dem nur Funktionswerte berechnet werden, ist die Powell-Methode . Um es zu verwenden, müssen Sie in der Minim-Funktion method = 'powell' setzen.

 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-Channo-Algorithmus (BFGS)

Um eine schnellere Konvergenz zur Lösung zu erhalten, verwendet das BFGS- Verfahren 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 normalerweise weniger Funktionsaufrufe als die Simplex-Methode.

Wir finden die Ableitung der Rosenbrock-Funktion in der analytischen Form:

$\ frac {\ partielle f} {\ partielle x_j} = \ summe \ grenzen_ {i = 1} ^ N 200 (x_i - x_ {i-1} ^ 2) (\ delta_ {i, j} - 2x_ {i -1, j}) - 2 (1 - x_ {i-1}) \ delta_ {i-1, j} =$

$= 200 (x_j - x_ {j-1} ^ 2) - 400x_j (x_ {j + 1} - x_j ^ 2) - 2 (1-x_j)$

Dieser Ausdruck gilt für Ableitungen aller Variablen mit Ausnahme der ersten und der letzten, die wie folgt definiert sind:

$\ frac {\ partielle f} {\ partielle x_0} = -400 x_0 (x_1 - x_0 ^ 2) - 2 (1 - x_0),$

$\ frac {\ partielle f} {\ partielle x_ {N-1}} = 200 (x_ {N-1} - x_ {N-2} ^ 2).$

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 der

Die Gradientenberechnungsfunktion wird wie unten gezeigt als Wert des jac-Parameters der Minim-Funktion angegeben.

 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)

Der Newton-Konjugatgradientenalgorithmus ist eine modifizierte Newton-Methode.
Newtons Methode basiert auf der Approximation einer Funktion in einer lokalen Region durch ein Polynom zweiten Grades:

$f \ left (\ mathbf {x} \ right) \ approx f \ left (\ mathbf {x} _ {0} \ right) + \ nabla f \ left (\ mathbf {x} _ {0} \ right) \ cdot \ left (\ mathbf {x} - \ mathbf {x} _ {0} \ right) + \ frac {1} {2} \ left (\ mathbf {x} - \ mathbf {x} _ {0} \ rechts) ^ {T} \ mathbf {H} \ links (\ mathbf {x} _ {0} \ rechts) \ links (\ mathbf {x} - \ mathbf {x} _ {0} \ rechts)$

wo $\ mathbf {H} \ left (\ mathbf {x} _ {0} \ right)$ ist eine Matrix von zweiten Ableitungen (Hessische Matrix, Hessische).
Wenn der Hessische positiv definit ist, kann das lokale Minimum dieser Funktion gefunden werden, indem der Nullgradient der quadratischen Form mit Null gleichgesetzt wird. Das Ergebnis ist ein Ausdruck:

$\ mathbf {x} _ {\ textrm {opt}} = \ mathbf {x} _ {0} - \ mathbf {H} ^ {- 1} \ nabla f$

Inverses Hessisches wird unter Verwendung der konjugierten Gradientenmethode berechnet. Ein Beispiel für die Verwendung dieser Methode zur Minimierung der Rosenbrock-Funktion ist unten angegeben. Um die Newton-CG-Methode verwenden zu können, müssen Sie eine Funktion definieren, die Hessisch auswertet.
Die hessische Funktion des Rosenbrock in analytischer Form ist gleich:

$H_ {i, j} = \ frac {\ partiell ^ 2 f} {\ partiell x_i x_j} = 200 (\ delta_ {i, j} - 2x_ {i-1} \ delta {i-1, j} - 400x_i (\ delta_ {i + 1, j} - 2x_i \ delta {i, j}) - 400 \ delta_ {i, j} (x_ {i + 1} - x_i ^ 2) + 2 \ delta_ {i, j } =$

$= (202 + 1200x_i ^ 2 - 400x_ {i + 1}) \ delta_ {i, j} - 400x_i \ delta_ {i + 1, j} - 400x_ {i-1} \ delta_ {i-1, j}$

wo $i, j \ in \ left [1, N-2 \ right]$ und $i, j \ in \ left [0, N-1 \ right]$ Bestimmen Sie die Matrix $N \ mal N$ .

Die verbleibenden Nicht-Null-Elemente der Matrix sind gleich:

$\ frac {\ partiell ^ 2 f} {\ partiell x_0 ^ 2} = 1200x_0 ^ 2 - 400x_1 +2$

$\ frac {\ partiell ^ 2 f} {\ partiell x_0 x_1} = \ frac {\ partiell ^ 2 f} {\ partiell x_1 x_0} = -400x_0$

$\ frac {\ partiell ^ 2 f} {\ partiell x_ {N-1} x_ {N-2}} = \ frac {\ partiell ^ 2 f} {\ partiell x_ {N-2} x_ {N-1 }} = -400x_ {N-2}$

$\ frac {\ partiell ^ 2 f} {\ partiell x_ {N-1} ^ 2} = 200x$

Beispielsweise hat im fünfdimensionalen Raum N = 5 die hessische Matrix für die Rosenbrock-Funktion die Bandform:

$\ tiny \ mathbf {H} = \ begin {bmatrix} 1200x_ {0} ^ {2} -400x_ {1} +2 & -400x_ {0} & 0 & 0 & 0 \\ -400x_ {0} & 202 + 1200x_ {1} ^ {2} -400x_ {2} & -400x_ {1} & 0 & 0 \\ 0 & -400x_ {1} & 202 + 1200x_ {2} ^ {2} -400x_ {3} & -400x_ {2} & 0 \\ 0 & & -400x_ {2} & 202 + 1200x_ {3} ^ {2} -400x_ {4} & -400x_ {3} \\ 0 & 0 & 0 & -400x_ { 3} & 200 \ end {bmatrix}$

Der Code, der dieses Hessische zusammen mit dem Code berechnet, um die Rosenbrock-Funktion unter Verwendung der konjugierten Gradientenmethode (Newton) zu minimieren:

 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 Funktion des Produkts des Hessischen und eines beliebigen Vektors

Bei realen Problemen kann das Berechnen und Speichern der gesamten hessischen Matrix erhebliche Zeit- und Speicherressourcen erfordern. Darüber hinaus besteht in der Tat keine Notwendigkeit, die hessische Matrix selbst anzugeben, da Das Minimierungsverfahren erfordert nur einen Vektor, der dem Produkt des Hessischen mit einem anderen beliebigen Vektor entspricht. Aus rechnerischer Sicht ist es daher sehr vorzuziehen, die Funktion, die das Ergebnis des Produkts des Hessischen mit einem beliebigen Vektor zurückgibt, sofort zu bestimmen.

Betrachten Sie die hess-Funktion, die einen Minimierungsvektor als erstes Argument und einen beliebigen Vektor als zweites Argument verwendet (zusammen mit anderen Argumenten der minimierten Funktion). In unserem Fall ist es nicht sehr schwierig, das Produkt der hessischen Rosenbrock-Funktion mit einem beliebigen Vektor zu berechnen. Wenn p ein beliebiger Vektor ist, dann das Produkt $H (x) \ cdot p$ hat die Form:

$\ mathbf {H} \ left (\ mathbf {x} \ right) \ mathbf {p} = \ begin {bmatrix} \ left (1200x_ {0} ^ {2} -400x_ {1} +2 \ right) p_ {0} -400x_ {0} p_ {1} \\ \ vdots \\ -400x_ {i-1} p_ {i-1} + \ left (202 + 1200x_ {i} ^ {2} -400x_ {i + 1} \ right) p_ {i} -400x_ {i} p_ {i + 1} \\ \ vdots \\ -400x_ {N-2} p_ {N-2} + 200p_ {N-1} \ end {bmatrix }.$

Eine Funktion, die das Produkt aus dem Hessischen 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: 66

Der Trust-Region-Algorithmus für konjugierte Gradienten (Newton)

Eine schlechte Konditionalität der hessischen Matrix und falsche Suchrichtungen können dazu führen, dass der Algorithmus für konjugierte Newton-Gradienten ineffizient sein kann. In solchen Fällen wird die Vertrauensbereichsmethode für konjugierte Newton-Gradienten bevorzugt.

Beispiel einer hessischen Matrixdefinition:

 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.]

Ein Beispiel mit der Produktfunktion des Hessischen 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)

 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 Krylovsky-Typ

Wie die Trust-NCG-Methode eignen sich Krylovsky-Methoden gut zur Lösung von Problemen im großen Maßstab, da sie nur Matrixvektorprodukte verwenden. Ihre Essenz besteht darin, das Problem im vertraulichen Bereich zu lösen, der durch den abgeschnittenen Krylov-Unterraum begrenzt ist. Für unsichere Aufgaben ist es besser, diese Methode zu verwenden, da im Vergleich zur Trust-NCG-Methode weniger nichtlineare Iterationen verwendet werden, da weniger Matrixvektorprodukte pro Unteraufgabe vorhanden sind. Darüber hinaus ist die Lösung für die quadratische Teilaufgabe genauer als die Trust-NCG-Methode.
Beispiel einer hessischen Matrixdefinition:

 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.]

Ein Beispiel mit der Produktfunktion des Hessischen und einem beliebigen Vektor:

 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.]

Vertrauensbasierter Näherungsalgorithmus

Alle Methoden (Newton-CG, trust-ncg und trust-krylov) eignen sich gut zur Lösung großer Aufgaben (mit Tausenden von Variablen). Dies liegt an der Tatsache, dass der zugrunde liegende Algorithmus für konjugierte Gradienten eine ungefähre Bestimmung der inversen hessischen Matrix impliziert. Die Lösung ist iterativ ohne explizite Zerlegung des Hessischen. Da nur die Funktion für das Produkt des Hessischen und eines beliebigen Vektors bestimmt werden muss, eignet sich dieser Algorithmus besonders für die Arbeit mit spärlichen (Banddiagonal-) Matrizen. Dies bietet niedrige Speicherkosten und erhebliche Zeiteinsparungen.

Bei mittelgroßen Problemen sind die Kosten für die Speicherung und Faktorisierung des Hessischen nicht kritisch. Dies bedeutet, dass eine Lösung in weniger Iterationen erhalten werden kann, wodurch die Unteraufgaben des Vertrauensbereichs fast genau aufgelöst werden. Dazu werden einige nichtlineare Gleichungen für jede quadratische Teilaufgabe iterativ gelöst. Eine solche Lösung erfordert normalerweise 3 oder 4 Zerlegungen der Holets-Hessischen Matrix. Infolgedessen konvergiert das Verfahren in weniger Iterationen und erfordert weniger Berechnung der Zielfunktion als andere implementierte Verfahren des Vertrauensbereichs. Dieser Algorithmus impliziert nur die Bestimmung der vollständigen hessischen Matrix und unterstützt nicht die Fähigkeit, die Produktfunktion des hessischen und eines beliebigen Vektors zu verwenden.

Ein Beispiel für die Minimierung der Rosenbrock-Funktion:

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

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 13 Function evaluations: 14 Gradient evaluations: 13 Hessian evaluations: 14 array([1., 1., 1., 1., 1.])

Dies verweilen vielleicht. 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 der Funktion einer Variablen, beliebige Minimierer und das Finden der Wurzeln des Gleichungssystems mithilfe des scipy.optimize-Pakets zu erzählen.

Quelle: https://docs.scipy.org/doc/scipy/reference/

SciPy, Optimierung