36 无约束优化方法

36.1 分块松弛法

对多元目标函数\(f(\boldsymbol x), \boldsymbol x \in \mathbb R^d\), 可以反复地分别沿每个坐标轴方向进行一维搜索, 称为坐标循环下降法。 按这种方法进行引申, 可以把\(\boldsymbol x\)的分量分成若干组, 每次固定其它分量不变, 针对某一组分量进行优化, 然后轮换其它组进行优化, 这样的方法叫做分块松弛法。 这种方法常常得到比较简单易行的算法, 只不过其收敛速度不一定是最优的。

例21.2 (k均值聚类) 设在\(\mathbb R^d\)中有\(n\)个点\(\boldsymbol x^{(i)}, i=1,\dots,n\), 要把这\(n\)个点按照距离分为\(k\)个类。 设\(C = \{ C_1, \dots, C_k \}\), \(C_j\)是分到第\(j\)个类中的点的集合, \(\boldsymbol\mu = \{ \boldsymbol\mu_1, \dots, \boldsymbol\mu_k \}\)是各类的中心(该类中的点的均值), 聚类问题就是求\(C\)\(\boldsymbol\mu\)使得 \[\begin{align} f(\boldsymbol\mu, C) = \sum_{j=1}^k \sum_{\boldsymbol x^{(i)} \in C_j} \| \boldsymbol x^{(i)} - \boldsymbol\mu_j \|^2 \end{align}\] 最小。

如果已知\(\boldsymbol\mu\),则对每个点计算该点到每个\(\boldsymbol\mu_j\)的距离, 把该点归类到距离最近的类中。 如果已知\(C\),则\(\boldsymbol\mu_j\)取第\(j\)类所有点的平均值。

\(\boldsymbol\mu\)\(C\)都未知的情况下求\(f(\boldsymbol\mu, \boldsymbol C)\)的最小值点,可以使用分块松弛法。 作为初始中心,可以随机地选\(k\)个不同的\(\boldsymbol x^{(i)}\)点当作初始的\(\boldsymbol\mu_j, j=1,\dots,k\), 但是最好能让这\(k\)个初始中心相互距离远一些。 然后,轮流地,先把所有点按照到类中心的距离分配到\(k\)个类中, 再重新计算类中心\(\{ \boldsymbol\mu_j \}\), 如此重复直到结果不再变动。 这个算法简单有效。

※※※※※

36.2 最速下降法

在目标函数\(f(\boldsymbol x)\)一阶可导时, 应利用导数(梯度)的信息, 向负梯度方向搜索前进, 使得每一步的目标函数值都减小。 基本的算法为:

最速下降法
取初值\(\boldsymbol x^{(0)} \in \mathbb R^d\), 置\(t \leftarrow 0\)
until (迭代收敛) {
\(\qquad\)\(\boldsymbol g^{(t)} \leftarrow -\nabla f(\boldsymbol x^{(t)})\)
\(\qquad\) 求最优步长\(\lambda_t > 0\)使得 \(f(\boldsymbol x^{(t)} + \lambda_t \boldsymbol g^{(t)}) = \min\limits_{\lambda>0} f(\boldsymbol x^{(t)} + \lambda \boldsymbol g^{(t)})\)
\(\qquad\) \(\boldsymbol x^{(t+1)} \leftarrow \boldsymbol x^{(t)} + \lambda_t \boldsymbol g^{(t)}\)
}
输出\(\boldsymbol x^{(t+1)}\)作为极小值点

在算法中,“迭代收敛”可以用梯度向量长度小于某个预定值, 或者两次迭代间函数值变化小于某个预定值来判断。 迭代中需要用某种一维搜索方法求步长\(\lambda_t\)

在适当正则性条件下, 最速下降法可以收敛到局部极小值点并且具有线性收敛速度。 但是, 最速下降法连续两次的下降方向\(-\nabla f(\boldsymbol x^{(t)})\)\(-\nabla f(\boldsymbol x^{(t+1)})\)是正交的, 这是因为在第\(t\)步沿\(-\nabla f(\boldsymbol x^{(t)})\) 搜索时找到点\(\boldsymbol x^{(t+1)}\)时已经使得函数值在该方向上不再变化, 下一个搜索方向如果与该方向不正交就不是下降最快的方向。 这样,最速下降法收敛速度比较慢。 所以, 在每次迭代的一维搜索时, 可以不要求找到一维搜索的极小值点, 而是函数值足够降低就结束一维搜索。 见沃尔夫准则

最速下降法和牛顿法

图36.1: 最速下降法和牛顿法

例36.1 求解无约束最优化问题 \[\begin{align} \mathop{\text{argmin}}_{(x_1, x_2) \in \mathbb R^2} 4 (x_1 - 1)^2 + (x_2 - 2)^4 . \end{align}\]

显然,函数有全局最小值点\(\boldsymbol x^* = (1, 2)\)

用最速下降法求解。易见 \[\begin{align*} \nabla f(\boldsymbol x) =& \left( 8(x_1 - 1), \; 4 (x_2 - 2)^3 \right)^T, \end{align*}\] 假设从点\(\boldsymbol x^{(0)} = (0, 0)\)出发,函数的等值线图以及迭代轨迹见图36.1中的实线, 轨迹中部分点坐标如下: \[\begin{align*} & \boldsymbol x^{(1)} = (0.641, 2.565), \ \boldsymbol x^{(2)} = (1.014, 2.471), \ \boldsymbol x^{(3)} = (0.962, 2.471), \ \boldsymbol x^{(4)} = (1.002, 2.271), \\ & \boldsymbol x^{(5)} = (0.986, 2.202), \ \boldsymbol x^{(6)} = (1.001, 2.197), \ \ldots, \ \boldsymbol x^{(10)} = (1.001, 2.156), \\ & \boldsymbol x^{(20)} = (1.001, 2.127), \ \boldsymbol x^{(40)} = (1.001, 2.099), \ \boldsymbol x^{(60)} = (1.001, 2.083), \ \boldsymbol x^{(80)} = (1.001, 2.074), \\ & \boldsymbol x^{(100)} = (1.001, 2.067), \ \boldsymbol x^{(150)} = (1.001, 2.055), \ \boldsymbol x^{(200)} = (1.001, 2.048) \end{align*}\]

※※※※※

从例子可以看出,最速下降法的相邻两次搜索方向正交, 而且越到最小值点附近接近得越慢。 所以,一维搜索不必找到一维的最小值点, 而是使得函数值足够降低就可以了。

最速下降法是其它方法的一个补充, 很少单独使用。

36.3 牛顿法

牛顿法利用海色阵和梯度向量求下一步, 不需要一维搜索, 对于二次多项式函数可以一步得到最小值点。 对一个存在二阶连续偏导的\(d\)元函数\(f(\boldsymbol x)\), 有如下的泰勒近似 \[\begin{align} f(\boldsymbol x) \approx f(\boldsymbol x^{(0)}) + \nabla f(\boldsymbol x_0)^T (\boldsymbol x - \boldsymbol x^{(0)}) + \frac12 (\boldsymbol x - \boldsymbol x^{(0)})^T \nabla^2 f(\boldsymbol x^{(0)}) (\boldsymbol x - \boldsymbol x^{(0)}), \tag{36.1} \end{align}\]\(\nabla^2 f(\boldsymbol x^{(0)})\)为正定矩阵,上式右侧的二次多项式函数的最小值点在 \[\begin{align*} \boldsymbol x^* = \boldsymbol x^{(0)} - \left[\nabla^2 f(\boldsymbol x^{(0)}) \right]^{-1} \, \nabla f(\boldsymbol x^{(0)}) \end{align*}\] 处达到。所以牛顿法取初值\(\boldsymbol x^{(0)}\)后,用公式 \[\begin{align} \boldsymbol x^{(t+1)} = \boldsymbol x^{(t)} - \left[\nabla^2 f(\boldsymbol x^{(t)}) \right]^{-1} \, \nabla f(\boldsymbol x^{(t)}), \ t=0,1,2,\dots \tag{36.2} \end{align}\] 进行迭代,直至收敛。 收敛的判断准则可以取为\(\| \nabla \boldsymbol f(x) \|\)足够小, 或者取为\(| f(\boldsymbol x^{(t+1)}) - f(\boldsymbol x^{(t)}) |\)足够小。

(36.2)在算法实现时, 应将逆矩阵表示改为如下的线性方程组求解: \[\begin{align} \nabla^2 f(\boldsymbol x^{(t)}) \boldsymbol d^{(t)} = -\nabla f(\boldsymbol x^{(t)}), \quad \boldsymbol x^{(t+1)} = \boldsymbol x^{(t)} - \boldsymbol d^{(t)} \tag{36.3} \end{align}\]

如果初值接近于最小值点并且目标函数满足一定正则性条件, 牛顿法有二阶收敛速度。

牛顿法的性能对海色阵\(H(\boldsymbol x)\)十分依赖。 如果\(H(\boldsymbol x^{(t)})\)不满秩, 则(36.3)的解不唯一, 算法无法进行。 另外,从(36.1)看出, 海色阵也是函数的局域性质的重要特征。 如果\(\boldsymbol x_*\)是一个稳定点, 即\(\nabla f(\boldsymbol x_*)=0\), 海色阵与此处目标函数的曲面形状密切相关。 如果\(H(\boldsymbol x_*)\)是正定的, 则\(\boldsymbol x_*\)是局部极小值点。 如果\(H(\boldsymbol x_*)\)是负定的, 则\(\boldsymbol x_*\)是局部极大值点。 如果\(H(\boldsymbol x_*)\)非奇异但是既有正特征值又有负特征值, 则\(\boldsymbol x_*\)是目标函数的鞍点,不是局部极值点。 如果\(H(\boldsymbol x_*)\)奇异(不满秩), 则\(\boldsymbol x_*\)不是极值点也不是鞍点。

使用牛顿法求最小值时, 最好的情况是海色阵始终为正定阵, 这时函数是严格凸函数, 极小值点也是全局最小值点。 在线性最小二乘问题中, 目标函数的海色阵是\(X^T X\)矩阵。

如果目标函数非负定但是接近于不满秩, 可以考虑给\(H(\boldsymbol x^{(t)})\)加一个对角的正定阵, 例如,岭回归就是将\(X^T X\)替换成\(X^T X + \lambda I\)

例36.2 再次考虑例36.1,用牛顿法。

易见\(f(\boldsymbol x)\)的海色阵为 \[\begin{align*} \nabla^2 f(\boldsymbol x) = \left(\begin{array}{cc} 8 & 0 \\ 0 & 12 (x_2 - 2)^2 \end{array}\right), \end{align*}\] 用牛顿法,迭代公式为 \[\begin{align*} \boldsymbol x^{(t)} = \boldsymbol x^{(t-1)} - \left[\nabla^2 f(\boldsymbol x^{(t-1)}) \right]^{-1} \, \nabla f(\boldsymbol x^{(t-1)}) = \boldsymbol x^{(t-1)} - (x_1 - 1, \; \frac13 (x_2 - 2))^T. \end{align*}\] 迭代过程见图36.1中的虚线。 轨迹中部分点坐标如下: \[\begin{align*} & \boldsymbol x^{(1)} = (1, 0.667), \ \boldsymbol x^{(2)} = (1, 1.111), \ \boldsymbol x^{(3)} = (1, 1.407), \ \boldsymbol x^{(4)} = (1, 1.605), \\ & \boldsymbol x^{(5)} = (1, 1.737), \ \boldsymbol x^{(6)} = (1, 1.824), \ \ldots, \ \boldsymbol x^{(10)} = (1, 1.965), \\ & \boldsymbol x^{(15)} = (1, 1.995), \ \boldsymbol x^{(20)} = (1, 1.999) \end{align*}\] 没有发生最速下降法那样反复折向前进的问题, 而且收敛速度相对较快。

※※※※※

36.4 拟牛顿法

牛顿法在目标函数具有较好光滑性而且海色阵正定时有很好的收敛速度。 但是,牛顿法需要目标函数有二阶导数, 还需要在每步迭代都计算海色阵的逆矩阵, 在海色阵非正定时可能会失败。

在公式(36.2)中, 把其中的海色阵\(\nabla^2 f(\boldsymbol x^{(t)})\) 替换成一个保证正定的矩阵\(H_t\), 则\(-H_t^{-1} \nabla f(\boldsymbol x^{(t)})\)\(f(\boldsymbol x)\)\(\boldsymbol x^{(t)}\)处的一个下降方向, 即当\(\lambda>0\)充分小而且\(\nabla f(\boldsymbol x^{(t)}) \neq \boldsymbol 0\)时一定有 \[\begin{align*} f\left( \boldsymbol x^{(t)} - \lambda H_t^{-1} \nabla f(\boldsymbol x^{(t)}) \right) < f(\boldsymbol x^{(t)}), \end{align*}\] 于是可以用一维搜索方法求步长\(\lambda_{t+1}\)使得 \[\begin{align} \begin{cases} \lambda_{t+1} =& \mathop{\text{argmin}}_{\lambda > 0} f\left( \boldsymbol x^{(t)} - \lambda H_t^{-1} \nabla f(\boldsymbol x^{(t)}) \right), \\ \boldsymbol x^{(t+1)} =& \boldsymbol x^{(t)} - \lambda_t H_t^{-1} \nabla f(\boldsymbol x^{(t)}). \end{cases} \tag{36.4} \end{align}\] 这样得到了保证函数值每次迭代都下降的算法。 这里正定矩阵\(H_t\)的选取有各种各样的方法, 比如,取\(H_t\)为单位阵,则(36.4)变成最速下降法。

(36.4)中取\(H_t\)为海色阵\(\nabla^2 f(\boldsymbol x^{(t)})\), 算法就变成增加了一维搜索的牛顿法, 称之为阻尼牛顿法。 其中的步长\(\lambda_{t+1}\)可以用一种一维搜索方法来求, 最简单的做法是, 依次考虑\(\lambda=1, 1/2, 1/4, \dots\), 一旦\(f\left( \boldsymbol x^{(t)} - \lambda H_t^{-1} \nabla f(\boldsymbol x^{(t)}) \right) < f(\boldsymbol x^{(t)})\)就停止搜索。 阻尼牛顿法也需要海色阵正定, 所以实际使用阻尼牛顿法时, 如果某步的海色阵不可逆或者\(\left[ \nabla f(\boldsymbol x^{(t)}) \right]^T \left[\nabla^2 f(\boldsymbol x^{(t)}) \right]^{-1} \nabla f(\boldsymbol x^{(t)}) \leq 0\), 可以在此步以\(-\nabla f(\boldsymbol x^{(t)})\)为搜索方向(最速下降法), 这样修改后可以使得阻尼牛顿法每步都减小目标函数值。

最速下降法、牛顿法、阻尼牛顿法在迭代中遇到不是局部最小值点的稳定点时都不能正常运行。 这时需要一些特殊的方法求得下降方向。 参见 徐成贤, 陈志平, and 李乃成 (2002) §3.2.2。

另外一些构造\(H_t\)的方法可以不用计算二阶偏导数而仅需要计算目标函数值和一阶偏导数值, 可以迭代地构造\(H_t\), 这样的算法有很多,统称为拟牛顿法变尺度法。 这样的\(H_t\)应该满足如下条件:

  1. \(H_t\)是对称正定阵,从而\(\boldsymbol d^{(t)} \stackrel{\triangle}{=} - H_t^{-1} \nabla f(\boldsymbol x^{(t)})\)是下降方向(称\(\boldsymbol d^{(t)}\)为拟牛顿方向)。
  2. \(H_{t+1}\)\(H_t\)迭代得到: \[\begin{align*} H_{t+1} = H_t + \Delta H_t, \end{align*}\]\(\Delta H_t\)是修正矩阵。
  3. \(H_{t+1}\)必须满足如下拟牛顿条件: \[\begin{align} H_{t+1} \boldsymbol\delta^{(t)} = \boldsymbol\zeta^{(t)} . \tag{36.5} \end{align}\] 其中\(\boldsymbol\delta^{(t)} \stackrel{\triangle}{=} \boldsymbol x^{(t+1)} - \boldsymbol x^{(t)}\), \(\boldsymbol\zeta^{(t)} \stackrel{\triangle}{=} \nabla f(\boldsymbol x^{(t+1)}) - \nabla f(\boldsymbol x^{(t)})\)

前两个条件比较自然,下面解释第三个条件的理由。

设迭代已经得到了\(H_t\)\(\boldsymbol x^{(t+1)}\), 需要给出下一步的\(H_{t+1}\)。 把\(f(\boldsymbol x)\)\(\boldsymbol x^{(t+1)}\)处作二阶泰勒展开得 \[\begin{align*} f(\boldsymbol x) \approx& f(\boldsymbol x^{(t+1)}) + \nabla f(\boldsymbol x^{(t+1)})^T (\boldsymbol x - \boldsymbol x^{(t+1)}) \nonumber\\ & + \frac12 (\boldsymbol x - \boldsymbol x^{(t+1)})^T \nabla^2 f(\boldsymbol x^{(t+1)}) (\boldsymbol x - \boldsymbol x^{(t+1)}), \end{align*}\] 这意味着 \[\begin{align*} \nabla f(\boldsymbol x) \approx \nabla f(\boldsymbol x^{(t+1)}) + \nabla^2 f(\boldsymbol x^{(t+1)}) (\boldsymbol x - \boldsymbol x^{(t+1)}), \end{align*}\] 在上式中取\(\boldsymbol x = \boldsymbol x^{(t)}\)\[\begin{align*} \nabla f(\boldsymbol x^{(t)}) \approx \nabla f(\boldsymbol x^{(t+1)}) + \nabla^2 f(\boldsymbol x^{(t+1)}) (\boldsymbol x^{(t)} - \boldsymbol x^{(t+1)}) , \end{align*}\] 则有 \[\begin{align*} \nabla^2 f(\boldsymbol x^{(t+1)}) (\boldsymbol x^{(t+1)} - \boldsymbol x^{(t)}) \approx \nabla f(\boldsymbol x^{(t+1)}) - \nabla f(\boldsymbol x^{(t)}) = \boldsymbol\zeta^{(t)}, \end{align*}\] 因为\(H_{t+1}\)是用来代替海色阵的,所以根据上式需要\(H_{t+1}\)满足拟牛顿条件(36.5)

修正矩阵\(\Delta H_t\)有多种取法, 其中一种比较高效而稳定的算法称为BFGS算法 (Broyden-Fletcher-Goldfarb-Shanno), 在多元目标函数无约束优化问题中被广泛应用。 BFGS算法的修正公式为 \[\begin{align} H_{t+1} = H_t + \frac{\boldsymbol\zeta^{(t)} (\boldsymbol\zeta^{(t)})^T} {(\boldsymbol\zeta^{(t)})^T \boldsymbol\delta^{(t)}} - \frac{(H_t \boldsymbol\delta^{(t)}) \; (H_t \boldsymbol\delta^{(t)})^T}{(\boldsymbol\delta^{(t)})^T H_t \boldsymbol\delta^{(t)}}, \tag{36.6} \end{align}\] 一般取初值\(H_0\)为单位阵, 这样按(36.4)(36.6)迭代, 如果一维搜索采用精确一维搜索或沃尔夫准则, (36.6)中的分母\((\boldsymbol\zeta^{(t)})^T \boldsymbol\delta^{(t)}\)可以保证为正值, 使得\(\{ H_t \}\)总是对称正定矩阵, 搜索方向\(-H_t^{-1} \nabla f(\boldsymbol x^{(t)})\)总是下降方向。 在适当正则性条件下BFGS方法具有超线性收敛速度。 详见 徐成贤, 陈志平, and 李乃成 (2002) §3.3。

(36.4)中计算\(H_t^{-1} \nabla f(\boldsymbol x^{(t)})\)要解一个线性方程组, 需要\(O(d^3)\)阶的计算量(\(d\)\(\boldsymbol x\)的维数)。 事实上,在BFGS算法中\(H_t^{-1}\)可以递推计算而不需直接求逆或求解线性方程组。 记\(V_t = H_t^{-1}\), 这些逆矩阵有如下递推公式: \[\begin{align} V_{t+1} =& V_t + \left[ 1 + \frac{(\boldsymbol\zeta^{(t)})^T V_t \boldsymbol\zeta^{(t)}} {(\boldsymbol\zeta^{(t)})^T \boldsymbol\delta^{(t)}} \right] \frac{\boldsymbol\delta^{(t)} (\boldsymbol\delta^{(t)})^T} {(\boldsymbol\zeta^{(t)})^T \boldsymbol\delta^{(t)}} - \frac{H_t \boldsymbol\zeta^{(t)} (\boldsymbol\delta^{(t)})^T + \left[ H_t \boldsymbol\zeta^{(t)} (\boldsymbol\delta^{(t)})^T \right]^T} {(\boldsymbol\zeta^{(t)})^T \boldsymbol\delta^{(t)}} \end{align}\]

36.5 Nelder-Mead方法

在导数不存在或很难获得时, 求解多元函数无约束最小值经常使用Nelder-Mead单纯型方法 (Nelder and Mead 1965)单纯型(simplex)指在\(\mathbb R^d\)空间中的\(d+1\)个点作为顶点的图形构成的凸集, 比如,在\(\mathbb R^2\)中给定了三个顶点的一个三角形是一个单纯型, 在\(\mathbb R^3\)中给定了四个顶点的一个四面体是一个单纯型。

Nelder-Mead方法是一种直接搜索算法,不使用导数或数值导数。 算法开始时要找到\(d+1\)个点\(\boldsymbol x_0, \boldsymbol x_1, \dots, \boldsymbol x_d\)构成一个单纯型 (要求这\(d+1\)个点不能位于一个超平面内), 计算相应的目标函数值\(y_j = f(\boldsymbol x_j), j=0,1,\dots,d\)。 然后,进行反复迭代, 对当前的单纯型进行变换使得目标函数值变小。 迭代在单纯型变得足够小或者函数值的变化变得足够小时结束, 保险起见,如果迭代次数超过一个预定次数就认为算法失败而结束。 单纯型可以沿着目标函数下降的方向延伸, 在遇到一个峡谷时可以改变方向, 在最小值点附近可以收缩。

选取初值时,可以任意选取初值\(\boldsymbol x_0\), 然后其他\(d\)个点可以从\(\boldsymbol x_0\)出发分别沿\(d\)个正坐标轴方向前进一段距离得到。 初始的单纯型不要取得太小, 否则很容易停留在附近的局部极小值点。

Nelder-Mead算法反射变换图示,其它变换类似

图36.2: Nelder-Mead算法反射变换图示,其它变换类似

在算法迭代的每一步, 假设当前单纯型为\((\boldsymbol x_0, \boldsymbol x_1, \dots, \boldsymbol x_d)\), 相应的函数值为\((y_0, y_1, \dots, y_d)\), 首先从中找到最坏的、次坏的、最好的函数值\(y_h, y_s, y_l\), 如果每次得到新的单纯型都把单纯型\(d+1\)个点次序重排使其函数值由小到大排列, 则\(y_l = y_0, y_s = y_{d-1}, y_h = y_d\)。 然后,计算除最坏点之外的\(d\)个顶点的重心\(\boldsymbol c = \frac{1}{d} \sum_{j \neq h} \boldsymbol x_j\)。 然后对单纯型进行变换以减小函数值。

单纯型的变换有反射、延伸(expansion)、退缩(contraction)和缩小(shrinkage)等四种, 变化量分别由四个参数\(\alpha>0, \beta \in (0, 1), \gamma>\alpha, \delta \in (0,1)\)控制, 这四个参数一般取为\(\alpha=1, \beta=\frac12, \gamma=2, \delta=\frac12\)

首先考虑反射是否有效, 计算反射点\(\boldsymbol x_r \stackrel{\triangle}{=} \boldsymbol c + \alpha (\boldsymbol c - \boldsymbol x_h)\), 这实际是把最坏点\(\boldsymbol x_h\)与中心\(\boldsymbol c\)连线后把连线延长得到点\(\boldsymbol x_r\)(见图), 这样点\(\boldsymbol x_h, \boldsymbol c, \boldsymbol x_r\) 在一条直线上但是最坏点\(\boldsymbol x_h\)与反射点\(\boldsymbol x_r\)分别在\(\boldsymbol c\)的两侧。 这样沿最坏点的反方向搜索期望能改善目标函数值。 如果这时\(y_r \stackrel{\triangle}{=} f(\boldsymbol x_r)\)满足\(y_l \leq y_r < y_s\) (反射点比原次坏点好、但不优于原最好点), 则接受\(\boldsymbol x_r\)为本迭代步的新顶点, 替换原来的最坏点\(\boldsymbol x_h\)即可。

如果反射点\(\boldsymbol x_r\)的目标函数值比原最好点还小(\(y_r < y_l\)), 则考虑把生成\(\boldsymbol x_r\)的延长线进一步延伸, 得到延伸点\(\boldsymbol x_e \stackrel{\triangle}{=} \boldsymbol c + \gamma (\boldsymbol x_r - \boldsymbol c)\)。 令\(y_e = f(\boldsymbol x_e)\), 如果比反射点进一步改善(\(y_e < y_r\)), 则接受延伸点\(\boldsymbol x_e\)为本迭代步的新顶点, 在单纯型中替换原来的最坏点\(\boldsymbol x_h\)即可, 本步迭代结束。 如果延伸后\(\boldsymbol x_e\)的函数值不如反射点\(\boldsymbol x_r\)的函数值(\(y_e \geq y_r\)), 则接受反射点\(\boldsymbol x_r\)为本迭代步的新顶点, 在单纯型中替换原来的最坏点\(\boldsymbol x_h\)即可, 本步迭代结束。

在尝试反射\(\boldsymbol x_r\)后如果发现其不优于次坏点\(\boldsymbol x_s\)(\(y_r \geq y_s\)), 则反射点不能使用。 这时, 可以在\(\boldsymbol x_h, \boldsymbol c, \boldsymbol x_r\)构成的线段上另找一个点, 看这个点是否可以接受, 这样的变换称为退缩。 新的点可以在\(\boldsymbol c\)\(\boldsymbol x_r\)之间(单纯型外部), 也可以在\(\boldsymbol x_h\)\(\boldsymbol c\)之间(单纯型内部)。 这时, 如果反射点的函数值\(y_r\)介于原次坏点与最坏点之间(\(y_s \leq y_r \leq y_h\)), 那么单纯型外部有希望改善目标函数, 所以取\(\boldsymbol x_c = \boldsymbol c + \beta(\boldsymbol x_r - \boldsymbol c)\)(在单纯型外部); 否则,若\(y_r\)比原最坏点还差(\(y_r > y_h\)), 则只能在单纯型内部考虑, 取\(\boldsymbol x_c = \boldsymbol c - \beta(\boldsymbol c - \boldsymbol x_h)\)(在单纯型内部)。 得到\(\boldsymbol x_c\)后令\(y_c = f(\boldsymbol x_c)\), 如果比反射点有改善(\(y_c < y_r\)), 则接受退缩点\(\boldsymbol x_c\)为本迭代步的新顶点, 在单纯型中替换原来的最坏点\(\boldsymbol x_h\)即可, 本步迭代结束。 否则,就要执行第四种变换:缩小。

当反射、延伸、退缩都失败时, 执行缩小变换。 缩小变换是保持原来的最好点\(\boldsymbol x_l\)不动, 其它点都按比例向\(\boldsymbol x_l\)收缩: \(\boldsymbol x_j \leftarrow \boldsymbol x_l + \delta (\boldsymbol x_j - \boldsymbol x_l)\), \(j \neq l\)。单纯型缩小后希望其它\(d\)个点的目标函数值能有所改善。

Nelder-Mead方法的收敛性很难确定。

在R软件中,用optim函数求解多元函数无约束优化问题, 该函数提供了BFGS、Nelder-Mead、共轭梯度、模拟淬火等算法。

习题

习题1

编写牛顿法的程序, 当没有输入偏导数函数时,程序用数值微分方法近似计算偏导数。

习题2

编写BFGS方法的程序, 当没有输入偏导数函数时,程序用数值微分方法近似计算偏导数。

习题3

编写Nelder-Mead算法的程序。

References

Nelder, J., and R. Mead. 1965. “A Simplex Method for Function Minimization.” Comput. J.

徐成贤, 陈志平, and 李乃成. 2002. 近代优化方法. 科学出版社.