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

例37.1 (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 \}\), 如此重复直到结果不再变动。 这个算法简单有效。

※※※※※

在目标函数\(f(\boldsymbol x)\)一阶可导时， 应利用导数（梯度）的信息， 向负梯度方向搜索前进， 使得每一步的目标函数值都减小。 最速下降法沿负梯度方向搜索找到函数值下降最多的步长， 基本的算法为:

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

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

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

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

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

用最速下降法求解。易见 \[\begin{aligned} \nabla f(\boldsymbol x) =& \left( 8(x_1 - 1), \; 4 (x_2 - 2)^3 \right)^T, \end{aligned}\] 假设从点\(\boldsymbol x^{(0)} = (0, 0)\)出发，函数的等值线图以及迭代轨迹见图37.1中的实线, 轨迹中部分点坐标如下： \[\begin{aligned} & \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{aligned}\]

使用R函数optim求解（不是最速下降法）：

这里正定矩阵\(H_t\)的选取有各种各样的方法， 比如，取\(H_t\)为单位阵，则(37.8)变成最速下降法。

在(37.8)中取\(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。

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

第二条可以替换成\(H_{t+1}^{-1}\)可以用\(H_{t}^{-1}\)递推得到。 前两个条件比较自然，下面解释第三个条件的理由。

设迭代已经得到了\(H_t\)和\(\boldsymbol x^{(t+1)}\)， 需要给出下一步的\(H_{t+1}\)。 把\(f(\boldsymbol x)\)在\(\boldsymbol x^{(t+1)}\)处作二阶泰勒展开得 \[\begin{aligned} 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{aligned}\] 这意味着 \[\begin{aligned} \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{aligned}\] 在上式中取\(\boldsymbol x = \boldsymbol x^{(t)}\)得 \[\begin{aligned} \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{aligned}\] 则有 \[\begin{aligned} \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{aligned}\] 因为\(H_{t+1}\)是用来代替海色阵的，所以根据上式需要\(H_{t+1}\)满足拟牛顿条件(37.9)。

修正矩阵\(\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{37.10} \end{align}\] 一般取初值\(H_0\)为单位阵， 这样按(37.8)和(37.10)迭代， 如果一维搜索采用精确一维搜索或沃尔夫准则， (37.10)中的分母\((\boldsymbol\zeta^{(t)})^T \boldsymbol\delta^{(t)}\)可以保证为正值， 使得\(\{ H_t \}\)总是对称正定矩阵， 搜索方向\(-H_t^{-1} \nabla f(\boldsymbol x^{(t)})\)总是下降方向。 在适当正则性条件下BFGS方法具有超线性收敛速度。 详见 徐成贤, 陈志平, and 李乃成 (2002) §3.3。

在(37.8)中计算\(H_t^{-1} \nabla f(\boldsymbol x^{(t)})\)要解一个线性方程组， 需要\(O(d^3)\)阶的计算量(\(d\)是\(\boldsymbol x\)的维数)。 事实上，在BFGS算法中\(H_t^{-1}\)可以递推计算而不需直接求逆或求解线性方程组。 记\(V_t = H_t^{-1}\), 这些逆矩阵有如下递推公式: \[\begin{aligned} 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{aligned}\]

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

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\)个正坐标轴方向前进一段距离得到。 初始的单纯型不要取得太小， 否则很容易停留在附近的局部极小值点。

图37.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\)。 然后对单纯型进行变换以减小函数值。

单纯型的变换有反射、退缩(contraction)、延伸(expansion)和缩小(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\)(见图37.2)， 这样点\(\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\)个点的目标函数值能有所改善。

算法的停止法则可以取为最后的\(d+1\)个定点处的\(d+1\)个函数值的样本标准差小于某一预先设定的精度阈值， 或迭代次数超过某一界限后宣布算法失败。

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

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

退火是一种机械加工技术， 加热金属零件到一定温度后进行加工， 然后控制温度缓慢地下降到常温。 在高温时， 分子运动速度较快， 在缓慢降温时可以调整机械应力、减少缺陷。

模拟退火算法的名称是利用了“缓慢降温”的做法和“温度”参数， 这种算法利用随机性进行全局优化， 用类似于Metropolis-Hasting的MCMC算法那样的方法在目标函数定义域进行随机游动， 试图找到全局最小值点。 算法从当前位置\(\boldsymbol x^{(t)}\)按照某种抽样分布生成随机数\(\boldsymbol \epsilon^{(t)}\), 令\(\boldsymbol x' = \boldsymbol x^{(t)} + \boldsymbol \epsilon^{(t)}\)， 记\(\Delta y = f(\boldsymbol x') - f(\boldsymbol x^{(t)})\)， 以如下概率令\(\boldsymbol x^{(t+1)} = \boldsymbol x'\): \[ p^{(t+1)} = \begin{cases} 1, & \text{若} \Delta y \leq 0, \\ \exp(- \Delta y / T^{(t+1)}), & \text{若} \Delta y > 0 . \end{cases} \] 否则令\(\boldsymbol x^{(t+1)} = \boldsymbol x^{(t)}\)。 即随机移动一步， 如果目标函数值变小， 则确认移动； 如果目标函数值变大， 就以一定正概率移动， 否则停留不动，等待下一次随机转移。 这里\(T^{(t+1)}\)是转移概率的参数， 类似于退火时的温度， \(T^{(t+1)}\)越大， 则\(p^{(t+1)}\)越大， 逆向的随机跳转越频繁， 可以用来从局部极小值点跳开； 随着迭代的进行， \(T^{(t+1)}\)应该逐步变小， 使得逆向随机跳转变稀少， 这样迭代能最终收敛。

\(T^{(t+1)}\)应选取为趋于0的数列， 在一定条件下如果取 \[ T^{(t)} = T^{(0)} \frac{\ln 2}{\ln (t+1)}, \] 可以确保迭代收敛到全局最小值点， 但显然这个速度太慢了。 实际应用中往往用 \[ T^{(t+1)} = \gamma T^{(t)}, \] 其中\(0 < \gamma < 1\)， 或取\(T^{(t)} = T^{(0)} / t\)。

\(\boldsymbol \epsilon^{(t+1)}\)的分布可取为多元标准正态分布， 或取与各自变量的取值更匹配的一元分布合成多元分布。

参见(Kochenderfer 2019)第8章。

模拟退火算法是利用随机移动来寻找全局最优点。 遗传算法是另一类算法的典型代表， 它产生多个搜索点， 利用这些搜索点以及搜索点之间的交互去接近全局最优点。

遗传算法利用了生物学中择优配对、交叉遗传、遗传突变的思想使得整个种群（搜索点集合）不断进化（接近全局最优点）。 设要求\(\min_{\boldsymbol x \in \mathbb R^d} f(\boldsymbol x)\)， 算法的核心包括：

详见(Kochenderfer 2019)第9章。 这种利用多个搜索点的方法还有许多， 例如“粒子蜂拥”方法，记录每一个搜索点的位置、速度、历史最优点， 更新速度为倾向于自己的历史最优点和群体的历史最优点。 “萤火虫”方法， 是每个搜索点都逐次被每一个优于自己的点吸引。

这些利用群体随机变化寻找全局最优的方法有利用查找比较广泛的定义域， 但是收敛到最小值点比较慢， 所以可以配合一些局部搜索算法； 可以对每一个搜索点都作局部搜索， 并更新搜索点和目标函数值； 也可以对每一个搜索点都作局部搜索后， 仅更新其搜索得到的局部更优的目标函数值， 但并不更新搜索点本身， 这种做法更能避开局部最小点的吸引。

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

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

编写Nelder-Mead算法的程序。

设\(f(\boldsymbol x)\)偏导数连续， 证明最速下降法的两个方向\(\boldsymbol g^{(t)}\)和\(\boldsymbol g^{(t+1)}\)正交。