Cet article décrit la méthode Nelder – Mead, connue sous le nom d'algorithme d'optimisation.
Le but est.
La méthode Nelder – Mead est un algorithme [^ 1] publié par John A. Nelder et Roger Mead en 1965. Recherchez la valeur minimale de la fonction pendant le déplacement. ([Wikipédia](https://ja.wikipedia.org/wiki/%E3%83%8D%E3%83%AB%E3%83%80%E3%83%BC%E2%80%93%E3% À partir de 83% 9F% E3% 83% BC% E3% 83% 89% E6% B3% 95)) Par exemple, si le déterminant est un problème dimensionnel $ 2 $, la solution optimale sera recherchée en déplaçant le triangle.
L'algorithme spécifique est le suivant. (Reportez-vous à la référence [^ 2])
Mettez à jour le point $ x_h $ qui donne la plus grande valeur de fonction parmi les sommets $ n + 1 $. À ce moment, les points candidats de mise à jour suivants sont calculés en utilisant le centre de gravité $ c $ de $ n $ sommets à l'exclusion de $ x_h $.
Si aucun de ces points candidats n'est bon, tous les points sauf $ x_ \ ell $ sont contractés plus près de $ x_ \ ell $. (RÉTRÉCIR)
https://codesachin.wordpress.com/2016/01/16/nelder-mead-optimization/ Les chiffres de ce blog sont faciles à comprendre pour Reflect, Expand, Contract et SHRINK.
En Python, spécifiez method = 'Nelder-Mead'
dans scipy.optimize.minimize Il peut être utilisé par.
Cependant, dans cet article, nous voulons utiliser tous les sommets du triangle pour créer une image GIF, nous l'avons donc implémentée comme suit.
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
Nous l'avons également comparé à l'implémentation de Scipy. ($ \ mathop {\ mathrm {minimiser}} _ {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]}')
Résultat d'exécution
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]
Une image GIF a été créée à l'aide de matplotlib.animation.FuncAnimation. Au moment de la mise en œuvre, j'ai fait référence à l'article suivant.
Tout d'abord, calculez les sommets du triangle à utiliser. Comme dans l'exemple précédent, la fonction objectif est $ 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)
Cela enregistre les sommets «maxiter + 1» dans «X».
Ensuite, créez une image GIF en utilisant FuncAnimation
.
FuncAnimation (fig, func, frames, fargs)
crée une image GIF avecfunc (frames [i], * fragments)
comme une seule image.
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')
Image GIF créée
Recommended Posts