Nelder Mead SIMPLEX Algorithm

Nelder-Mead (Downhill Simplex Method) 算法最早由 Jone Nelder 和 Roger Mead 于 1965 年提出，是一种基于启发式规则的优化算法，类似常见的遗传算法（Generic Algorithm，GA）和粒子群算法（Particle Swarm Optimization，PSO），通过人为设计的一系列规则，从初始值出发，迭代寻找最优解。像众多启发式算法一样，Nelder-Mead 不需要了解函数的具体形式，不利用梯度信息，同样也无法保证结果的最优性，只能找到“较”好的可行解.

单纯形

Nelder-Mead 算法基于单纯性的概念构造解的迭代优化策略。对于 N 维的优化问题，初始构造一个 N+1 维的单纯形，计算单纯形顶点的函数值，然后分析比较顶点函数值，构造新的顶点和单纯形，直至达到收敛条件。

单纯形（SIMPLEX）是几何学中的概念，是2维中的三角形、3维中的四面体向任意维度的扩展。之所以叫 SIMPLEX，是因为它定义了一定维度空间中最简单的（Simplest）多面体形式。

0-SIMPLEX 是一个点；
1-SIMPLEX 是一根直线；
2-SIMPLEX 是一个三角形；
3-SIMPLEX 是一个四面体；
更严格地说，k-SIMPLEX 是一个 k 维多面体，是组成它的 k+1 个顶点构成的凸包（Convex hull）。

算法

对于$n$维最小化问题:

$$\min f(x)$$
利用 Nelder-Mead 算法搜索最小值的方法如下：

STEP-1 初始化：初始化 $n+1$个点 $x_1,...,x_{n+1}$ ，作为n-SIMPLEX的顶点；

STEP-2 排序（Order）：根据$f(x)$ 值对顶点进行重排序， $f(x_1)\le f(x_2)\le ...\le f(x_{n+1})$；检查是否截止；

截止条件可以选择点的方差作为标准，选择一定的 tolerance 作为界限；

STEP-3 重心（Centroid）：抛弃最差点 $x_{n+1}$ ，计算前 $n$ 个点的重心 $x_o=\frac{1}{n}\sum _{i=1}^n{x}_i$；

STEP-4 反射（Reflection）：计算反射点$x_r=x_o+\rho (x_o-x_{n+1})$ ；如果$f(x_r)$优于 $f(x_n)$但是差于 $f(x_1)$ ，i.e. $f(x_1)\le f(x_r)\le f(x_n)$ ，则用 $x_r$ 替换$x_{n+1}$ 构建新的 n-SIMPLEX ，继续 STEP-2；
因为 $x_{n+1}$ 是当前最差的点，因此我们有可能在 $x_{n+1}$和 $x_o$ 的反向延长线上找到一个优于 $x_{n+1}$的点；
STEP-5 扩展（Expansion）：如果反射点是最优点，i.e. $f(x_r)<f(x_1)$ ，则计算扩展点 $x_e=x_o+\gamma (x_r-x_o)$ ；如果扩展点优于反射点，i.e.$f(x_e)<f(x_r)$，将 $x_{n+1}$ 替换为 $x_e$ ，然后继续 STEP-2；否则，将 $x_{n+1}$替换为 $x_r$ ，然后继续 STEP-2；

因为$x_r$ 是当前最优点，因此我们有可能在$x_r$和$x_o$之间找到一个更好的点；

STEP-6 收缩（Contraction）：如果​$f(x_n)<f(x_r)<f(x_{n+1})$，计算收缩点 $x_c=x_o+\alpha (x_r-x_o)$​。如果​ $f(x_c)\le f(x_{n+1})$ ，则将​ $x_{n+1}$ 替换为 $x_c$ ​，然后继续 STEP-3；否则，进入 STEP-7。如果​ $f(x_c)\ge f(x_{n+1})$ ，计算内收缩点​ $x_{cc}=x_o+\alpha (x_{n+1}-x_o)$ 。如果内收缩点优于最差点，则用内收缩点 $x_{cc}$​替代最差点；否则，进入STEP-7；

如果反射点差于第二差点但是优于最差点，我们希望能够在单纯形中找到一个更好的点；如果反射点差于最差点，我们希望能够在最差点和重心之间找到一个更好的点；
STEP-7 回退（Shrink）：将除了当前最优点以外的点全部用 $x_i:=x_1+\sigma (x_i-x_1)$ 替换掉，然后继续STEP-2；

如果从重心向最差的点收缩后，点依然变差，保留最优点作为收缩中心，将所有的点向最优点收缩；
以上算法中， $\rho ,\gamma ,\alpha ,\sigma$分别为反射、扩展、收缩、回退系数，取值一般为 $\rho =1,\gamma =2,\alpha =1/2,\sigma =1/2$ 。

Nelder-Mead 算法的简单实现

基于 Scipy 中 Nelder-Mead 的算法实现，进行了简化，仅保留主算法部分作为算法展示

# !/usr/lib/env python3
import numpy as np
from numpy import asfarray


def minimize_neldermead(func, x0, maxiter=100, initial_simplex=None,
                        xatol=1e-4, fatol=1e-4, adaptive=False, bounds=None):
    """
    Args:
        bounds: Sequence of ``(min, max)`` pairs for each element in `x`. None
                is used to specify no bound.
    """

    x0 = asfarray(x0).flatten()

    # 设定反射 rho、扩展 gamma、收缩 psi、回退 sigma 系数
    # 系数可以为固定 / 自适应
    if adaptive:
        dim = float(len(x0))
        rho = 1
        gamma = 1 + 2/dim
        alpha = 0.75 - 1/(2*dim)
        sigma = 1 - 1/dim
    else:
        rho = 1
        gamma = 2
        alpha = 0.5
        sigma = 0.5

    # 设定解的上下界
    if bounds is not None:
        dim = len(x0)
        lower_bound = [bounds[i][0] for i in range(dim)]
        upper_bound = [bounds[i][1] for i in range(dim)]
    
    # 将初始解根据上下界进行截断
    if bounds is not None:
        x0 = np.clip(x0, lower_bound, upper_bound)

    # 设置初始单纯形
    # 如果未给定初始单纯形：在初始解x0的基础上在各个维度伸缩 nonzdelt = 0.05
    nonzdelt = 0.05
    zdelt = 0.00025
    if initial_simplex is None:
        N = len(x0)
        sim = np.empty((N + 1, N), dtype=x0.dtype)
        sim[0] = x0
        for k in range(N):
            y = np.array(x0, copy=True)
            if y[k] != 0:
                y[k] = (1 + nonzdelt) * y[k]
            else:
                y[k] = zdelt
            sim[k + 1] = y
    else:
        sim = np.asfarray(initial_simplex).copy()
        N = sim.shape[1]

    # 将初始单纯形根据上下界截断
    if bounds is not None:
        sim = np.clip(sim, lower_bound, upper_bound)

    one2np1 = list(range(1, N + 1))
    fsim = np.empty((N + 1,), float)

    for k in range(N + 1):
        fsim[k] = func(sim[k])

    # 根据func(x)值排序，sim[0,:]为最优值
    ind = np.argsort(fsim)
    fsim = np.take(fsim, ind, 0)
    sim = np.take(sim, ind, 0)

    iterations = 1

    while iterations < maxiter:

        # 当顶点距离最大值或者目标值最大值小于一定值时，结束迭代
        if (np.max(np.ravel(np.abs(sim[1:] - sim[0]))) <= xatol and
                np.max(np.abs(fsim[0] - fsim[1:])) <= fatol):
            break
        
        # 重心 := 平均值
        xbar = np.add.reduce(sim[:-1], 0) / N

        # 反射值
        xr = (1 + rho) * xbar - rho * sim[-1]
        if bounds is not None:
            xr = np.clip(xr, lower_bound, upper_bound)
        fxr = func(xr)

        # 控制是否回退
        doshrink = 0

        if fxr < fsim[0]:  # 如果反射点为最优点，计算扩展点
            # 扩展点
            xe = (1 + rho * gamma) * xbar - rho * gamma * sim[-1]
            # 截断
            if bounds is not None:
                xe = np.clip(xe, lower_bound, upper_bound)
            fxe = func(xe)

            if fxe < fxr:  # 如果扩展点优于反射点，将最差点替换为扩展点
                sim[-1] = xe
                fsim[-1] = fxe
            else:  # 如果扩展点不是最优点，将最差点替换为反射点
                sim[-1] = xr
                fsim[-1] = fxr
        else:  # fsim[0] <= fxr, 反射点不是最优点
            if fxr < fsim[-2]:  # 如果反射点不是最优点，但是优于第二差点，替换最差点为反射点
                sim[-1] = xr 
                fsim[-1] = fxr
            else:  # fxr >= fsim[-2]，反射点不是最优点，但是优于第二差点
                if fxr < fsim[-1]:  # 反射点优于最差点
                    # 收缩点
                    xc = (1 + alpha * rho) * xbar - alpha * rho * sim[-1]
                    if bounds is not None:
                        xc = np.clip(xc, lower_bound, upper_bound)
                    fxc = func(xc)

                    if fxc <= fxr: # 如果收缩点优于反射点
                        sim[-1] = xc 
                        fsim[-1] = fxc
                    else:
                        doshrink = 1
                else:  # 反射点差于最差点
                    # 内收缩点
                    xcc = (1 - alpha) * xbar + alpha * sim[-1]
                    if bounds is not None:
                        xcc = np.clip(xcc, lower_bound, upper_bound)
                    fxcc = func(xcc)

                    if fxcc < fsim[-1]:  # 如果内收缩点优于最差点，替代最差点
                        sim[-1] = xcc
                        fsim[-1] = fxcc
                    else:
                        doshrink = 1

                if doshrink:  # 回退
                    for j in one2np1:
                        sim[j] = sim[0] + sigma * (sim[j] - sim[0])
                        if bounds is not None:
                            sim[j] = np.clip(sim[j], lower_bound, upper_bound)
                        fsim[j] = func(sim[j])

        ind = np.argsort(fsim)
        sim = np.take(sim, ind, 0)
        fsim = np.take(fsim, ind, 0)
        iterations += 1

    x = sim[0]
    fval = np.min(fsim)

    return x, fval

示例

Rosenbrock函数

$$f(\mathbf{x})=\sum _{i=1}^{N-1}100(x_{i+1}-x_i^2{)}^2+(1-x_i{)}^2$$
$x_i=1$时，该函数取得最小值0。

用Nelder-Mead Simplex算法优化

import numpy as np
from scipy.optimize import minimize
def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2]) #初始值
res = minimize(rosen, x0, method='nelder-mead',
               options={'xatol': 1e-8, 'disp': True})

最后

Nelder Mead 是一个非常简单有效的启发式算法，一般而言效率低于基于一阶梯度或者二阶梯度的算法，对于低维度优化问题是一个不错的初始尝试，对于高维度的问题，效果可能不是很好。

参考 & 扩展阅读

Gao, F. and Han, L. Implementing the Nelder-Mead simplex algorithm with adaptive parameters. 2012. Computational Optimization and Applications. 51:1, pp. 259-277

Nelder, J.A. and Mead, R. (1965), “A simplex method for function minimization”, The Computer Journal, 7, pp. 308-313

Wright, M.H. (1996), “Direct Search Methods: Once Scorned, Now Respectable”, in Numerical Analysis 1995, Proceedings of the 1995 Dundee Biennial Conference in Numerical Analysis, D.F. Griffiths and G.A. Watson (Eds.), Addison Wesley Longman, Harlow, UK, pp. 191-208.