Nelder-Mead 算法

贡献者： addis

预备知识　Matlab 的函数

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

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

先给 $x_{i}$ 点按照 $f (x_{i})$ 从小到大的顺序重新排序，使 $i$ 越大 $f (x_{i})$ 越大。
计算前 $N$ 个点的平均位置为
$\begin{matrix} (1) & m = \frac{1}{N} \sum_{i = 1}^{N} x_{i} . \end{matrix}$
计算 $x_{N + 1}$ 关于点 $m$ 的对称点为
$\begin{matrix} (2) & r = 2 m - x_{N + 1} . \end{matrix}$
如果 $f (x_{1}) ⩽ f (r) < f (x_{N})$ ，令 $x_{N + 1} = r$ ，并进入下一个循环。
如果 $f (r) < f (x_{1})$ ，计算拓展点为
$\begin{matrix} (3) & s = m + 2 (m - x_{N + 1}) . \end{matrix}$
如果 $f (s) < f (r)$ ，令 $x_{N + 1} = s$ 并进入下一个循环，否则令 $x_{N + 1} = r$ 。并进入下一循环。
如果 $f (x_{N}) ⩽ f (r) < f (x_{N + 1})$ ，令
$\begin{matrix} (4) & c_{1} = m + (r - m) / 2 . \end{matrix}$
如果 $f (c_{1}) < f (r)$ ，令 $x_{N + 1} = c_{1}$ 并进入下一循环，否则执行最后一步。
如果 $f (x_{N + 1}) ⩽ f (r)$ 令
$\begin{matrix} (5) & c_{2} = m + (x_{N + 1} - m) / 2 . \end{matrix}$
如果 $f (c_{2}) < f (x_{N + 1})$ ，令 $x_{N + 1} = c_{2}$ 并进入下一循环，否则执行最后一步。
令
$\begin{matrix} (6) & v_{i} = x_{1} + (x_{i} - x_{1}) / 2 (i = 2 \dots N + 1) . \end{matrix}$
并用 $v_{i}$ 赋值给 $x_{i}$ ，进入下一循环。

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

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

　　以下是该算法的 Matlab 代码。

代码 1：NelderMead.m

% f 是函数句柄，只接受一个 N 维行矢量作为输入变量， 并返回一个函数值
% x0 是 N 维行矢量， xerr 是 xmin 各个元素的绝对误差
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
% 主循环
x_last = x*0;
mask = true(1, N+1); % 改变的顶点
while true
    if max(max(abs(x(2:end,:) - x(1,:)))) < xerr % 判断误差
        break;
    elseif all(x(:) == x_last(:))
        warning('NelderMead: abs err too small, machine precision reached');
        break;
    else
        x_last = x;
    end
    % 求值并排序
    for ii = find(mask)
        y(ii) = f(x(ii,:));
    end
    [y, order] = sort(y);
    x = x(order,:);
    m = mean(x(1:N,:)); % 平均位置
    r = 2*m - x(N+1,:); % 反射点
    f_r = f(r);
    mask(:) = false;
    mask(end) = true;
    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;
        mask(jj) = true;
    end
end
% 输出变量
xmin = x(1,:);
fmin = f(xmin);
end

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

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。