Nelder-Mead 算法

　　 Nelder-Mead 算法是一种求多元函数局部最小值的算法，其优点是不需要函数可导并能较快收敛到局部最小值．Matlab 自带的 fminsearch 函数就是使用该算法．对 $N$ 元函数 $f( \boldsymbol{\mathbf{x}} )$（这里把函数自变量用 $N$ 维矢量来表示），该算法需要提供函数自变量空间中的一个初始点，$ \boldsymbol{\mathbf{x}} _1$，算法从该点出发寻找局部最小值，以下是具体步骤．

　　 我们先根据初始点另外生成 $N$ 个初始点 $ \boldsymbol{\mathbf{x}} _2\dots \boldsymbol{\mathbf{x}} _{N + 1}$，使 $ \boldsymbol{\mathbf{x}} _{1 + i}$ 在第 $i$ 个分量比 $ \boldsymbol{\mathbf{x}} _1$ 的大 5%，其他分量保持相同．如果 $ \boldsymbol{\mathbf{x}} _1$ 的第 $i$ 个分量为零，那么 $ \boldsymbol{\mathbf{x}} _{1 + i}$ 的第 $i$ 个分量设为 $0.00025$．得到 $N+1$ 个初始点后，开始按照以下步骤进行循环，直到满足特定的精度条件时退出循环．

1. 先给这些点按照 $f( \boldsymbol{\mathbf{x}} _i)$ 从小到大的顺序重新排序并重命名，使 $i$ 越大 $f( \boldsymbol{\mathbf{x}} _i)$ 越大．
2. 计算前 $N$ 个点的平均位置为
   \begin{equation} \boldsymbol{\mathbf{m}} = \frac1N \sum_{i=1}^N \boldsymbol{\mathbf{x}} _i \end{equation}
3. 计算 $ \boldsymbol{\mathbf{x}} _{N + 1}$ 关于点 $ \boldsymbol{\mathbf{m}} $ 的反射点为
   \begin{equation} \boldsymbol{\mathbf{r}} = 2 \boldsymbol{\mathbf{m}} - \boldsymbol{\mathbf{x}} _{N + 1} \end{equation}
4. 如果 $f( \boldsymbol{\mathbf{x}} _1) \leqslant f( \boldsymbol{\mathbf{r}} ) < f( \boldsymbol{\mathbf{x}} _N)$，令 $ \boldsymbol{\mathbf{x}} _{N+1} = \boldsymbol{\mathbf{r}} $，并进入下一个循环．
5. 如果 $f( \boldsymbol{\mathbf{r}} ) < f( \boldsymbol{\mathbf{x}} _1)$，计算拓展点为
   \begin{equation} \boldsymbol{\mathbf{s}} = \boldsymbol{\mathbf{m}} + 2( \boldsymbol{\mathbf{m}} - \boldsymbol{\mathbf{x}} _{N+1}) \end{equation}
   如果 $f( \boldsymbol{\mathbf{s}} ) < f( \boldsymbol{\mathbf{r}} )$，令 $ \boldsymbol{\mathbf{x}} _{N+1} = \boldsymbol{\mathbf{s}} $ 并进入下一个循环，否则令 $ \boldsymbol{\mathbf{x}} _{N+1} = \boldsymbol{\mathbf{r}} $．并进入下一循环．
6. 如果 $f( \boldsymbol{\mathbf{x}} _N) \leqslant f( \boldsymbol{\mathbf{r}} ) < f( \boldsymbol{\mathbf{x}} _{N+1})$，令
   \begin{equation} \boldsymbol{\mathbf{c}} _1 = \boldsymbol{\mathbf{m}} + ( \boldsymbol{\mathbf{r}} - \boldsymbol{\mathbf{m}} )/2 \end{equation}
   如果 $f( \boldsymbol{\mathbf{c}} _1) < f( \boldsymbol{\mathbf{r}} )$，令 $ \boldsymbol{\mathbf{x}} _{N + 1} = \boldsymbol{\mathbf{c}} _1$ 并进入下一循环，否则执行最后一步．
7. 如果 $f( \boldsymbol{\mathbf{x}} _{N+1}) \leqslant f( \boldsymbol{\mathbf{r}} )$ 令
   \begin{equation} \boldsymbol{\mathbf{c}} _2 = \boldsymbol{\mathbf{m}} + ( \boldsymbol{\mathbf{x}} _{N+1} - \boldsymbol{\mathbf{m}} )/2 \end{equation}
   如果 $f( \boldsymbol{\mathbf{c}} _2) < f( \boldsymbol{\mathbf{x}} _{N+1})$，令 $ \boldsymbol{\mathbf{x}} _{N+1} = \boldsymbol{\mathbf{c}} _2$ 并进入下一循环，否则执行最后一步．
8. 令
   \begin{equation} \boldsymbol{\mathbf{v}} _i = \boldsymbol{\mathbf{x}} _1 + ( \boldsymbol{\mathbf{x}} _i - \boldsymbol{\mathbf{x}} _1)/2 \qquad (i = 2\dots N+1) \end{equation}
   并用 $ \boldsymbol{\mathbf{v}} _i$ 赋值给 $ \boldsymbol{\mathbf{x}} _i$，进入下一循环．

　　 观察以上步骤可知，当局部最小值的位置在 $N+1$ 个围成的图形以外时，图形倾向于变大且加速向最小值移动．当最小值的位置在图形内部时，图形倾向于缩小．随着循环次数增加，这 $N+1$ 个点最终将向局部最小值聚拢．

　　 我们可以在每个循环的第一步之后计算 $ \boldsymbol{\mathbf{x}} _1$ 和 $ \boldsymbol{\mathbf{x}} _{N+1}$ 的距离来估算自变量的误差，如果该误差小于某个值，即可结束循环并使用 $ \boldsymbol{\mathbf{x}} _1$ 作为最终结果．作为另一种方法，我们也可以在每个循环的第一步之后计算 $f( \boldsymbol{\mathbf{x}} _{N+1}) - f( \boldsymbol{\mathbf{x}} _1)$ 来估算最小值的误差．

　　 以下是该算法的 Matlab 代码．

% f 是函数句柄，只接受一个 N 维行矢量作为输入变量， 并返回一个函数值
% x0 是 N 维行矢量， xerr 是一个标量
function [xmin, fmin] = NelderMead(f, x0, xerr)
N = numel(x0); % f 是 N 元函数
x = zeros(N+1,N); % 预赋值
y = zeros(1,N+1);
% 计算 N+1 个初始点
x(1,:) = x0;
for ii = 1:N
    x(ii+1,:) = x(1,:);
    if x(1,ii) == 0
        x(ii+1,ii) = 0.00025;
    else
        x(ii+1,ii) = 1.05 * x(1,ii);
    end
end
% 主循环
for kk = 1:10000
    % 求值并排序
    for ii = 1 : N + 1
        y(ii) = f(x(ii,:));
    end
    [y, order] = sort(y);
    x = x(order,:);
    if norm(x(N+1,:) - x(1,:)) < xerr % 判断误差
        break;
    end
    m = mean(x(1:N,:)); % 平均位置
    r = 2*m - x(N+1,:); % 反射点
    f_r = f(r);
    if y(1) <= f_r && f_r < y(N) % 第 4 步
        x(N+1,:) = r; continue;
    elseif f_r < y(1) % 第 5 步
        s = m + 2*(m - x(N+1,:));
        if f(s) < f_r
            x(N+1,:) = s;
        else
            x(N+1,:) = r;
        end
        continue;
    elseif f_r < y(N+1) % 第 6 步
        c1 = m + (r - m)*0.5;
        if f(c1) < f_r
            x(N+1,:) = c1; continue;
        end
    else % 第 7 步
        c2 = m + (x(N+1,:) - m)*0.5;
        if f(c2) < y(N+1)
            x(N+1,:) = c2; continue;
        end
    end
    for jj = 2:N+1 % 第 8 步
        x(jj,:) = x(1,:) + (x(jj,:) - x(1,:))*0.5;
    end
end
% 输出变量
xmin = x(1,:);
fmin = f(xmin);
end

　　 该程序中有几个需要注意的地方．首先，主循环并没有用 while 1，而是用 for 循环指定一个最大循环次数．这是为了避免少数情况下可能发生的死循环（例如 $f( \boldsymbol{\mathbf{x}} )$ 在某个区域中的值处处相等时）．第二，4-7 步中对 $f( \boldsymbol{\mathbf{r}} )$ 的判断有且仅有一个成立，所以我们可以用 if...elseif...else 结构来选择．最后，4-5 步的情况下程序必定会执行 continue 语句而跳过第 8 步，只有 6-7 步中的 if 判断失败程序才会执行第 8 步．