优化方法总结续篇：下降单纯形法（downhill simplex） 及python示例代码

优化方法总结续篇：下降单纯形法（downhill simplex） 及python示例代码

下降单纯形法(downhill simplex method)是一个广泛使用的“derivative free”的优化算法。一般来说它的效率不高，但是文献[1]提到“the downhill simplex method may frequently be the *best* method to use if the figure of merit is “get something working quickly” for a problem whose computational burden is small.”

单纯形法的基本思路是在\(N\)维空间中，构造一个非退化的初始单纯形，然后做一系列的几何操作，如反射、扩展、收缩等，逐步往极值点移动该单纯形。由于这些几何操作的目的基本上都是让单纯形往极小值移动，所以叫下降单纯形法。

假设待优化的函数为\(f(\mathbf{x})\)，\(N\)维空间里的单纯形\(Z\)的\(N+1\)个顶点按照函数值从小到大排列分别为\(\mathbf{x}_{0},\mathbf{x}_{2},\cdots,\mathbf{x}_{N}\)，定义

\(\bar{\mathbf{x}}=\sum_{i=0}^{N-1}\mathbf{x}_{i}\)

为\(Z\)中除了顶点\(\mathbf{x}_{N}\)之外其余顶点的中心点。

连接\(\bar{\mathbf{x}}\)和\(\mathbf{x}_{N}\)的直线公式可以写成：

\(\bar{\mathbf{x}}(t)=(1-t)\bar{\mathbf{x}}+t\mathbf{x}_{N}\)

下降单纯形法就是从沿着直线$\bar{\mathbf{x}}(t)$方向的几个特殊步长寻找\(\mathbf{x}_{N}\)的替代点，使该替代点处的函数值比\(\mathbf{x}_{N}\)更小，如果没有找到这种替代点，那么就把除了\(\mathbf{x}_{0}\)点之外的其余点往\(\mathbf{x}_{0}\)靠拢。

假设是要求函数的极小值，可以把对应函数值越小的点认为越好，越大的点认为越差。下降单纯形法每步迭代过程简述如下：

1. 首先计算最差点沿着直线$\bar{\mathbf{x}}(t)$关于平均点$\bar{\mathbf{x}}$的对称点
2. 如果对称点介于最好和次差点之间，那么就接受它(reflection)；
3. 如果对称点比最好点还好，那么做沿该方向更大胆的尝试，令\(t=-2\)，如果新尝试点比对称点更好则接受新尝试点(expand)，否则接受当前对称点(reflection)；
4. 如果对称点介于次差点和最差点之间，那么沿该方向做更小心的尝试，即令\(t=-0.5\)，如果新尝试点比对称点更好则接受新尝试点(outside contraction)
5. 如果对称点比最差点还差，那么沿反方向做尝试，即令\(t=0.5\)，如果新尝试点比对称点更好则接受新尝试点(inside contraction)
6. 如果4和5均失败，即对称点比次差还要差而且outside contraction与inside contraction均失败，那么把最好点之外的其他点都朝最好点收缩(shrink)

上述过程如果用区间图表示会更清晰，区间的三个分界点就是最好、 次差和最差点。对应的伪代码可以参考文献 [3] 第9.5节。

上述迭代算法需要提供一个初始单纯形，该单纯形可以参考文献[1]给出的方法得到。首先任选一点\(\mathbf{x}_{0}\)，然后利用以下公式：

\(\mathbf{x}_{i}=\mathbf{x}_{0}+\lambda\mathbf{e}_{i}\)

其中\(\mathbf{e}_{i}\)代表\(N\)维空间的单位矢量，\(\lambda\)表示步长。


另外，在实际应用中还有一个很重要的诀窍，即重启(restart)。因为单纯形在迭代更新的时候很容易就卡在某个中间位置上，这时单纯形的 *最好* 和 *最差* 点几乎相同，单纯形的体积收缩的很小，会大大减慢迭代速度。为了解决这个问题，可以合理设置初始单纯形的大小。更有效的，就是可以在单纯形卡住的时候通过重新初始化单纯形来加快收敛速度。在利用初始化公式的时候，把当前单纯形的 *最好* 点作为\(\mathbf{x}_{0}\)保留下来，这样保证重启就不会影响之前已经计算的结果。


下面给出了单纯形法求解一个简单问题的python实现，其中待优化函数为\(f(\mathbf{x})=(\mathbf{x}-\mathbf{x}_{0})(\mathbf{x}-\mathbf{x}_{1})\)

首先定义一些需要用到的函数

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns

def vertice_init(vertex_0, step_length):
    '''
    initialize vertice of the simplex
    using the following formula:
    $xi=x0+step_length*ei$
    '''

    emat = np.eye(vertex_0.size) * step_length
    vertice = [vertex_0]
    for ii in range(vertex_0.size):
        vertice.append(vertex_0 + emat[:, ii])
    return vertice


def f(v):
    '''
    Evaluation of Function $f$
    '''
    dim = v.size
    v0 = np.ones(dim) * 5
    v1 = np.ones(dim) * 3
    return 0.5 * np.dot(v - v0, v - v1)


def line(t, v1, v2):
    return (1 - t) * v1 + t * v2

接下来定义算法主函数，注意里面的restart部分：

def simplex(f, vertice, maxit=1000, step_length=100, tol=1e-3):
    vertice_max_list = []  # store the max vertex during each iteration
    vertice_min_list = []  # store the min vertex during each iteration
    for jj in range(maxit):
        y = []
        for ii in vertice:
            y.append(f(ii))
        y = np.array(y)
        #  only the highest (worst), next-highest, and lowest (best) vertice
        # are neeed
        idx = np.argsort(y)
        vertice_max_list.append(vertice[idx[-1]])
        vertice_min_list.append(vertice[idx[0]])
        
        # centroid of the best n vertice
        # NOTE: the worst vertex should be excluded, but for simplicity we don't do this
        v_mean = np.mean(vertice)

        # compute the candidate vertex and corresponding function vaule
        v_ref = line(-1, v_mean, vertice[idx[-1]])
        y_ref = f(v_ref)
        if y_ref >= y[idx[0]]