[Python] Implementierung der Nelder-Mead-Methode und Speichern von GIF-Bildern durch matplotlib

In diesem Artikel wird die Nelder-Mead-Methode beschrieben, die als Optimierungsalgorithmus bezeichnet wird.

Implementieren Sie die Nelder-Mead-Methode in Python
Geben Sie das Verhalten des Algorithmus als GIF-Bild mit matplotlib aus

Der Zweck ist.

Überblick und Implementierung der Nelder-Mead-Methode

Überblick über den Algorithmus

Die Nelder-Mead-Methode ist ein Algorithmus [^ 1], der 1965 von John A. Nelder und Roger Mead veröffentlicht wurde. Dabei handelt es sich um eine $ n $ -Dimensionale Einheit (Simplex), die wie eine Amöbe aus $ n + 1 $ Eckpunkten besteht. Suchen Sie während der Bewegung nach dem Mindestwert der Funktion. ([Wikipedia](https://ja.wikipedia.org/wiki/%E3%83%8D%E3%83%AB%E3%83%80%E3%83%BC%E2%80%93%E3% Von 83% 9F% E3% 83% BC% E3% 83% 89% E6% B3% 95)) Wenn es sich bei der Determinante beispielsweise um ein $ 2 $ -Dimensionsproblem handelt, wird beim Verschieben des Dreiecks nach der optimalen Lösung gesucht.

Der spezifische Algorithmus ist wie folgt. (Siehe Referenz [^ 2])

Nelder–Mead

Aktualisieren Sie den Punkt $ x_h $, der den größten Funktionswert unter den Eckpunkten $ n + 1 $ ergibt. Zu diesem Zeitpunkt werden die folgenden Aktualisierungskandidatenpunkte unter Verwendung des Schwerpunkts $ c $ von $ n $ Eckpunkten ohne $ x_h $ berechnet.

Punkt symmetrisch um den Schwerpunkt bewegt $ c $ $ x_r $ (Reflect: $ x_r --c = - (x_h --c) $)
Punkt, der sich in symmetrischer Richtung zum Schwerpunkt ausdehnt $ c $ $ x_e $ (Erweitern: $ x_e --c = -2 (x_h --c) $)
Schrumpfpunkt in der Nähe des Schwerpunkts $ c $ $ x_c $ (Vertrag: $ x_c --c = \ pm 0,5 (x_h --c) $)

Wenn keiner dieser Kandidatenpunkte gut ist, werden alle Punkte außer $ x_ \ ell $ näher an $ x_ \ ell $ zusammengezogen. (SCHRUMPFEN)

https://codesachin.wordpress.com/2016/01/16/nelder-mead-optimization/ Die Zahlen in diesem Blog sind für Reflect, Expand, Contract und SHRINK leicht zu verstehen.

Implementierung in Python

Geben Sie in Python "method =" Nelder-Mead "in [scipy.optimize.minimize] an (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html). Es kann von verwendet werden. In diesem Artikel möchten wir jedoch alle Eckpunkte des Dreiecks verwenden, um ein GIF-Bild zu erstellen. Daher haben wir es wie folgt implementiert.

from typing import Callable, Tuple, Union

import numpy as np


def _order(x: np.ndarray, ordering: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    indices = np.argsort(ordering)
    return x[indices], ordering[indices]

def optimize(
    fun: Callable,
    x0: np.ndarray,
    maxiter: Union[int, None] = None,
    initial_simplex: Union[np.ndarray, None] = None
):
    if x0.ndim != 1:
        raise ValueError(f'Expected 1D array, got {x0.ndim}D array instead')

    # initialize simplex
    if initial_simplex is not None:
        if initial_simplex.ndim != 2:
            raise ValueError(f'Expected 2D array, got {x0.ndim}D array instead')
        x = initial_simplex.copy()
        n = x[0].size
    else:
        h = lambda x: (x[0][x[1]] != 0) * (0.05 - 0.00025) + 0.00025
        n = x0.size
        x = np.array([x0 + h([x0, i]) * e for i, e in enumerate(np.identity(n))] + [x0])

    if maxiter is None:
        maxiter = 200 * n

    # parameters
    alpha = 1.0
    gamma = 2.0
    rho = 0.5
    sigma = 0.5

    # order
    fx = np.array(list(map(fun, x)))
    x, fx = _order(x, fx)

    # centroid
    xo = np.mean(x[:-1], axis=0)
    n_inv = 1 / n

    for _ in range(maxiter):
        fx1 = fx[0]
        fxn = fx[-2]
        fxmax = fx[-1]
        xmax = x[-1]

        xr = xo + alpha * (xo - xmax)
        fxr = fun(xr)

        if fx1 <= fxr and fxr < fxn:
            # reflect
            x[-1] = xr
            fx[-1] = fun(xr)
            x, fx = _order(x, fx)
            xo = xo + n_inv * (xr - x[-1])

        elif fxr < fx1:
            xe = xo + gamma * (xo - xmax)
            fxe = fun(xe)
            if fxe < fxr:
                # expand
                x = np.append(xe.reshape(1, -1), x[:-1], axis=0)
                fx = np.append(fxe, fx[:-1])
                xo = xo + n_inv * (xe - x[-1])
            else:
                # reflect
                x = np.append(xr.reshape(1, -1), x[:-1], axis=0)
                fx = np.append(fxr, fx[:-1])
                xo = xo + n_inv * (xr - x[-1])

        else:
            if fxr > fxmax:
                xc = xo + rho * (xmax - xo)
            else: 
                xc = xo + rho * (xr - xo)
                fxmax = fxr
            if fun(xc) < fxmax:
                # contract
                x[-1] = xc
                fx[-1] = fun(xc)
                x, fx = _order(x, fx)
                xo = xo + n_inv * (xc - x[-1])
            else:
                # shrink
                x[1:] = (1 - sigma) * x[0] + sigma * x[1:]
                fx[1:] = np.array(list(map(fun, x[1:])))
                x, fx = _order(x, fx)
                xo = np.mean(x[:-1], axis=0)

    return x, fx

Wir haben es auch mit der Implementierung von Scipy verglichen. ($ \ mathop {\ mathrm {minimieren}} _ {x, y} \ quad f (x, y) = x ^ 2 + y ^ 2 $)

from scipy.optimize import minimize

maxiter = 25

fun = lambda x: x @ x
x0 = np.array([0.08, 0.08])

# scipy
%time res = minimize(fun=fun, x0=x0, options={'maxiter': maxiter}, method='Nelder-Mead')
xopt_scipy = res.x

# implemented
%time xopt, _ = optimize(fun=fun, x0=x0, maxiter=maxiter)

print('\n')
print(f'Scipy: {xopt_scipy}')
print(f'Implemented: {xopt[0]}')

Ausführungsergebnis

CPU times: user 1.49 ms, sys: 41 µs, total: 1.53 ms
Wall time: 1.54 ms
CPU times: user 1.64 ms, sys: 537 µs, total: 2.18 ms
Wall time: 1.86 ms


Scipy: [-0.00026184 -0.00030341]
Implemented: [ 2.98053651e-05 -1.26493496e-05]

Erstellen von GIF-Bildern mit matplotlib

Ein GIF-Bild wurde mit matplotlib.animation.FuncAnimation erstellt. Zum Zeitpunkt der Implementierung habe ich auf den folgenden Artikel verwiesen.

Erstellen Sie ein GIF-Bild mit matplotlib.animation (FuncAnimation)

Berechnen Sie zunächst die Eckpunkte des zu verwendenden Dreiecks. Wie im vorherigen Beispiel ist die Zielfunktion $ f (x, y) = x ^ 2 + y ^ 2 $.

maxiter = 25

fun = lambda x: x @ x
x = np.array([[0.08, 0.08], [0.13, 0.08], [0.08, 0.13]])
X = [x]
for _ in range(maxiter):
    x, fx = optimize(fun, x[0], maxiter=1, initial_simplex=x)
    X.append(x)

Dies speichert "Maxiter + 1" Eckpunkte in "X".

Erstellen Sie als Nächstes ein GIF-Bild mit "FuncAnimation".

FuncAnimation (fig, func, frame, fargs) erstellt ein GIF-Bild mitfunc (frame [i], * fragmente)als einem Frame.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as pat
import matplotlib.animation as animation

def func(x, xmin, xmax, ymin, ymax, xx, yy, vals):
    # clear the current axes
    plt.cla()
    
    # set x-axis and y-axis
    plt.xlim([xmin, xmax])
    plt.ylim([ymin, ymax])
    plt.hlines(0, xmin=xmin, xmax=xmax, colors='gray')
    plt.vlines(0, ymin=ymin, ymax=ymax, colors='gray')
    
    # set aspect
    plt.gca().set_aspect('equal', adjustable='box')
    
    # draw filled contour
    plt.contourf(xx, yy, vals, 50, cmap='Blues')
    
    # draw triangle
    plt.axes().add_patch(pat.Polygon(x, ec='k', fc='m', alpha=0.2))
    
    # draw three vertices
    plt.scatter(x[:, 0], x[:, 1], color=['r', 'g', 'b'], s=20)

n_grid=100
delta=0.005
interval=300

xmax, ymax = np.max(X, axis=(0, 1)) + delta
xmin, ymin = np.min(X, axis=(0, 1)) - delta

# function values of lattice points
xx, yy = np.meshgrid(np.linspace(xmin, xmax, n_grid), np.linspace(ymin, ymax, n_grid))
vals = np.array([fun(np.array([x, y])) for x, y in zip(xx.ravel(), yy.ravel())]).reshape(n_grid, n_grid)

fig = plt.figure(figsize=(10, 10))
ani = animation.FuncAnimation(fig=fig, func=func, frames=X, fargs=(xmin, xmax, ymin, ymax, xx, yy, vals), interval=interval)
ani.save("nelder-mead.gif", writer = 'imagemagick')

GIF-Bild erstellt