引言

本节将继续介绍无约束优化问题的常用求解方法，包括牛顿法、拟牛顿法。

本节是对优化算法(十八~二十)牛顿法的理论补充，其中可能出现一些定理的证明过程这里不再赘述，并在相应位置附加链接。

回顾：最速下降算法的缺陷

在上一节中介绍了最速下降法，并提到一个缺陷：梯度下降法仅使用负梯度方向作为下降方向。

• 不可否认的是：在每一次迭代过程中，负梯度方向必然是当前迭代步骤的最速下降方向，但它仅仅是局部最优解；
• 这个最速下降方向仅使用了目标函数$f(\cdot )$在当前数值解$x_k$的一阶梯度信息；而实际上二阶梯度信息同样可以参与下降方向的判别。

经典牛顿法基本介绍

相比之下，牛顿法$(\text{Newton Method})$则使用一阶与二阶梯度信息共同判别下降方向。从目标函数的角度观察，可以理解为：对$x_k$处的二次逼近函数进行最小化。

使用泰勒展开式对目标函数$f(x)$进行二阶泰勒展开：

• 对于经典牛顿法$(\text{Pure Newton Method})$,仅设置步长${\alpha }_k=1$。与最速下降法相反，在牛顿法中我们更关注迭代过程中选择的方向，而非步长。
• 其中$x-x_k$表示下降方向；
  $$\begin{array}{rl}\min f(x)& =\min f[x_k+1\cdot (x-x_k)]\\ & =\min \left\{f(x_k)+\frac{1}{1!}\cdot [\nabla f(x_k){]}^T(x-x_k)+\frac{1}{2!}(x-x_k{)}^T{\nabla }^{2}f(x_k)(x-x_k)+\mathcal{O}(\Vert x-x_k{\Vert }^2)\right\}\\ & \approx \min \left\{f(x_k)+\frac{1}{1!}\cdot [\nabla f(x_k){]}^T(x-x_k)+\frac{1}{2!}(x-x_k{)}^T{\nabla }^{2}f(x_k)(x-x_k)\right\}\end{array}$$

从而直接对上述二元函数求解最小值。首先求解梯度$\nabla f(x)$：
$$\nabla f(x)=\nabla f(x_k)+\frac{1}{2}\cdot [{\nabla }^{2}f(x_k){]}^T\cdot 2(x-x_k)$$
令$\nabla f(x)\triangleq 0$，有：
$$\{\begin{array}{l}[{\nabla }^{2}f(x_k){]}^T(x-x_k)=-\nabla f(x_k)\\ \Rightarrow x=x_k-[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)\end{array}$$
很明显，该线搜索表达方式中：${\alpha }_k=1,{\mathcal{D}}_k=-[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)$。其对应算法迭代步骤这里不再赘述，见牛顿法与正则化一节。

经典牛顿法的问题

观察上述迭代公式：$x=x_k-[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)$

• 如果该迭代公式能够正常执行，必然需要满足${\nabla }^{2}f(x_k)\succ 0$，即$x_k$位置对应的$\text{Hessian Matrix}$必然是正定矩阵。由${\nabla }^{2}f(x_k)\succ 0$可以推出：二次逼近函数$\begin{array}{r}f(x)=f(x_k)+[\nabla f(x_k){]}^T(x-x_k)+\frac{1}{2}(x-x_k{)}^T{\nabla }^{2}f(x_k)(x-x_k)\end{array}$必然是凸函数，从而${\mathcal{D}}_k=-[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)$是下降方向，并且是全局最优解。
  这里可以解决优化算法(十八)经典牛顿法中的疑惑。并不是所有情况下$f(x)$是凸函数，开口向上有最小解;只有${\nabla }^{2}f(x_k)\succ 0$时才可以。
• 相反，如果$\text{Hessian Matrix }\Rightarrow {\nabla }^{2}f(x_k)$不是正定矩阵：
  • 如果存在${\nabla }^{2}f(x_k)$特征值为$0$——$[{\nabla }^{2}f(x_k){]}^{-1}$无法求解；
  • 如果存在${\nabla }^{2}f(x_k)$特征值为负值——那么观察$[\nabla f(x_k){]}^T{\mathcal{D}}_k$：
    在上一节介绍了$[\nabla f(x_k){]}^T{\mathcal{D}}_k$的物理意义：$x_k$所在位置关于方向向量${\mathcal{D}}_k$的方向导数;当$[\nabla f(x_k){]}^T{\mathcal{D}}_k<0$时，方向向量${\mathcal{D}}_k$是下降方向。
    $$[\nabla f(x_k){]}^T{\mathcal{D}}_k=-[\nabla f(x_k){]}^T[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)$$
    很明显，该结果是三个矩阵的线性组合，如果${\nabla }^{2}f(x_k)$中存在一个负特征值，加上开始的负号，其连加项中必然存在一个结果是正值，虽然不清楚该正值的具体大小，但能够肯定的是：$[\nabla f(x_k){]}^T{\mathcal{D}}_k$不能百分之百地确认其小于$0$，从而${\mathcal{D}}_k$未必是下降方向。

经典牛顿法的优点与缺陷

• 优点：
  • 当初始点$x_0$足够接近于收敛点$x^{\ast }$时，并且${\nabla }^{2}f(x)$满足较好性质时，其收敛速度是二阶收敛。
    相关证明详见优化算法(十九)经典牛顿法的收敛性分析，这里不再赘述。
  • 该方法具有二次终止性。
    如果使用牛顿法求解凸二次函数最小化问题时，不仅存在二次终止性,甚至可以实现一步终止。因为求得的下降方向是全局最优方向。
• 缺陷：
  • 首先，在迭代过程中只有更新位置的$\text{Hessian Matrix }\Rightarrow {\nabla }^{2}f(\cdot )\succ 0$时才能使用；
  • 并且由于$\text{Hessian Matrix}$是$n$阶方阵，其计算时间复杂度是$\mathcal{O}(k^3)$，计算量大。并且适用范围较窄。

经典牛顿法示例

这里依然使用最速下降法中的示例：最小化目标函数$\begin{array}{r}\min f(x,y)=\frac{1}{2}{x}^2+2y^2\end{array}$，设置初始点$x_0=(21{)}^T$对于凸二次函数的解法可以实现一步收敛。对应代码表示如下：

import numpy as np
import math
import matplotlib.pyplot as pltdef f(x, y):return 0.5 * (x ** 2) + 2 * (y ** 2)def Derf(x,y):return np.array([x,4 * y])def DoubleDerf():return np.array([[1.0,0.0],[0.0,4.0]])def ConTourFunction(x,Contour):return math.sqrt(0.5 * (Contour - (0.5 * (x ** 2))))def NewtomMethod(epsilon=0.001):Start = np.array([2.0,1.0])LocList = list()LocList.append(Start)alpha = 0.2NextList = list()while True:D = np.linalg.inv(DoubleDerf()).dot(Derf(Start[0],Start[1]))Next = Start - alpha * DNextList.append(Next)if np.linalg.norm(Derf(Next[0],Next[1])) <= epsilon:LocList.append(Next)breakelse:Start = NextLocList.append(NextList[-1])return LocListdef DrawPicture(LocList):ContourList = [0.1,0.2,0.5,1.0]LimitParameter = 0.0001plt.figure(figsize=(10, 5))for Contour in ContourList:# 设置范围时，需要满足x的定义域描述。x = np.linspace(-1 * math.sqrt(2 * Contour) + LimitParameter, math.sqrt(2 * Contour) - LimitParameter, 200)y1 = [ConTourFunction(i, Contour) for i in x]y2 = [-1 * j for j in y1]plt.plot(x, y1, '--', c="tab:blue")plt.plot(x, y2, '--', c="tab:blue")plotList = list()for (x, y) in LocList:plotList.append((x, y))plt.scatter(x, y, s=50, facecolor="none", edgecolors="tab:red", marker='o')if len(plotList) < 2:continueelse:plt.plot([plotList[0][0], plotList[1][0]], [plotList[0][1], plotList[1][1]], c="tab:red")plotList.pop(0)plt.show()if __name__ == '__main__':LocList = NewtomMethod()DrawPicture(LocList)

当alpha=1.0时，对应图像结果表示如下：
可以看出，选择准确的方向与合适的步长,即可达到一步收敛。

当然，如果给出的步长过小，导致收敛可能无法一步到位。例如：alpha=0.2时的图像结果表示如下：
在该示例中，由于目标函数是凸二次函数，因而它的二阶逼近函数就是目标函数自身。从而每次寻找的方向都是全局最优方向。

修正牛顿法介绍

针对经典牛顿法中的缺陷：
$$\{\begin{array}{l}{\mathcal{D}}_k=-[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)\\ {\alpha }_k=1\end{array}$$
针对步长${\alpha }_k$的修正：

• 修正原因：牛顿法是对目标函数的二阶逼近函数进行优化。而不是真正对目标函数自身进行优化。如果目标函数过于复杂，其对应泰勒展开式中高阶项系数无法忽视，这导致：仅仅表示成二阶逼近函数无法对目标函数进行充分表示。

  反之，如果仅使用二阶逼近函数表示复杂函数，会导致：虽然默认的${\alpha }_k=1$能够使二阶逼近函数有效收敛，但可能无法使目标函数有效收敛。

• 修正方法：由于方向必然是下降方向，那么基于该方向使用线搜索方法(如$\text{Wolfe}$准则)重新确定${\alpha }_k$。

针对方向${\mathcal{D}}_k$的修正：

• 方法一：正则化法。如果${\nabla }^{2}f(x_k)$不是正定矩阵，使用矩阵${\mathcal{B}}_k$进行替代：${\mathcal{D}}_k=-{\mathcal{B}}_k^{-1}\nabla f(x_k)$。
  关于正则化法的详细介绍见本节内链接牛顿法与正则化。
  $${\mathcal{B}}_k={\nabla }^{2}f(x_k)+\lambda \mathcal{I}$$
  其中$\mathcal{I}$表示单位矩阵；$\lambda$为适当正数为保持${\mathcal{B}}_k$正定。
  这也称作$\text{Hessian Matrix}$主对角线扰动；
• 方法二：正则化法的优化版：
  关于方法一的缺陷：$\lambda$的取值存在约束/技巧。假设${\nabla }^{2}f(x_k)$不是正定矩阵，并且其对应的特征值为${\lambda }_i(i=1,2,\cdots ,n)$,对应$\lambda$的取值必须满足：
  λ > max ⁡ i = 1 , 2 , ⋯ , n { − λ i } \lambda > \mathop{\max}\limits_{i=1,2,\cdots,n} \{-\lambda_i\} λ>i=1,2,⋯,nmax​{−λi​}
  这意味着：满足该条件的$\lambda$值可能会很大，并且是每一个对角线元素加上$\lambda$。从而导致原始特征值被$\lambda$被分掉相应权重。
  优化版：可以将${\nabla }^{2}f(x_k)$进行特征值分解
  其中$\text{Diag}({\tau }_i)$表示由特征值${\tau }_i$构成的对角阵。
  $${\nabla }^{2}f(x_k)={\mathcal{Q}}^T\text{Diag}({\tau }_i)\mathcal{Q}$$
  对于${\tau }_i(i=1,2,\cdots ,n)$可以执行如下操作：
  其中$\delta$是一适当的正数。
  $${\tau }_i=\{\begin{array}{l}{\tau }_i\text{if }{\tau }_i\ge \delta \\ \delta \text{Otherwise}\end{array}$$
  这种操作相比于方法一的优势在于：关于$>0$的特征值，没有修改它的特征信息，从而该维度依然受到${\nabla }^{2}f(x_k)$自身特征值的主导。
• 方法三：若存在某一步骤其${\nabla }^{2}f(x_k)$不是正定矩阵，可以在该步骤中直接使用最速下降法替代牛顿法。我们不否认：最速下降法的方向可能不优秀(局部最优),但它至少必然是下降方向。

拟牛顿法

拟牛顿法$(\text{Quasi-Newton Method})$与牛顿法相似，其都是考虑：目标函数$f(\cdot )$在$x_k$位置的二阶逼近函数。记该函数为$m_k(x)$，表示如下：
$$m_k(x)=f(x_k)+[\nabla f(x_k){]}^T(x-x_k)+\frac{1}{2}(x-x_k{)}^T{\mathcal{B}}_k(x-x_k){\mathcal{B}}_k\succ 0$$
但拟牛顿法与牛顿法的核心区别在于：牛顿法使用$\text{Hessian Matrix}\Rightarrow {\nabla }^{2}f(x_k)$作为二次型位置的矩阵：

• 其中${\nabla }^{2}f(x_k)\in {\mathbb{R}}^{n\times n}$，不容易计算；
• ${\nabla }^{2}f(x_k)$性质不够稳定，它可能是一个非正定矩阵。

而拟牛顿法则直接使用正定矩阵${\mathcal{B}}_k\succ 0$替代${\nabla }^{2}f(x_k)$。当然${\mathcal{B}}_k$仅仅是一个正定矩阵虽然在计算上便捷了，但我们同样期望它能存在如下性质：

• ${\mathcal{B}}_k$能够体现$m_k(x)$的一些二阶梯度信息；
• 相比于${\nabla }^{2}f(x_k)\in {\mathbb{R}}^{n\times n}$的不易获取，期望${\mathcal{B}}_k$相比之下能够更容易获取。

后续操作与牛顿法别无二致：
$$\{\begin{array}{l}\nabla m_k(x)=\nabla f(x_k)+{\mathcal{B}}_k(x-x_k)=0\\ \\ \Rightarrow x=x_k-{\mathcal{B}}_k^{-1}\nabla f(x_k)\end{array}$$
并记${\mathcal{D}}_k=-{\mathcal{B}}_k^{-1}\nabla f(x_k)$为使$\min m_k(x)$的下降方向。
其对应的方向导数：$[\nabla f(x_k){]}^T{\mathcal{D}}_k<0$恒成立。

拟牛顿法的算法过程

关于拟牛顿法的算法过程表示如下：

• 初始化：初始点$x_0$；描述终止条件的$ϵ$；$k=0$；以及初始状态下的正定矩阵${\mathcal{B}}_0$：
  我们既希望${\mathcal{B}}_0$包含该函数在$x_0$点处的二阶梯度信息,又希望该${\mathcal{B}}_0$是稳定的正定矩阵。例如根据上述方法，以${\nabla }^{2}f(x_0)$为基础，并对该矩阵进行修正。
  $${\mathcal{B}}_0={\nabla }^{2}f(x_0)+\lambda \cdot \mathcal{I}{\mathcal{B}}_0\succ 0$$
• 判断当前点$x_k$是否为收敛点：$\text{if }\Vert \nabla f(x_k)\Vert \le ϵ$，如果是，算法终止；反之，向下执行步骤；
• 计算下降方向：${\mathcal{D}}_k=-{\mathcal{B}}_k^{-1}\nabla f(x_k)$
• 基于下降方向，通过非线性搜索方法(例如$\text{Armijo},\text{Wolfe}$)来确定步长${\alpha }_k$；
• 计算下一个更新点$x_{k+1}$的位置：$x_{k+1}=x_k+{\alpha }_k\cdot {\mathcal{D}}_k$；并确定$x_{k+1}$对应的${\mathcal{B}}_{k+1}$。返回步骤$2$，如果满足条件，算法停止；反之，继续迭代。

核心问题很明显：${\mathcal{B}}_{k+1}$怎么求$?$只要${\mathcal{B}}_{k+1}$求出来，就可以得到相应的下降方向，从而持续迭代过程。
其中${\mathcal{B}}_{k+1}$的确定方式有很多种，从而对应不同的拟牛顿法。

矩阵${\mathcal{B}}_{k+1}$的获取方法

获取矩阵${\mathcal{B}}_{k+1}$的基本要求

关于矩阵${\mathcal{B}}_{k+1}$，它的基本要求是：
该方程也被称作拟牛顿方程。
$$\nabla f(x_{k+1})-\nabla f(x_k)={\mathcal{B}}_{k+1}(x_{k+1}-x_k)$$
为什么${\mathcal{B}}_{k+1}$需要满足该条件$?$

• 根据上述算法流程，完全可以确定$x_{k+1}$的具体位置，从而可以计算出目标函数$f(\cdot )$在该位置处的梯度信息：$\nabla f(x_{k+1})$；
  如果是正常牛顿法，我们需要继续计算$\text{Hessian Matrix }\Rightarrow {\nabla }^{2}f(x_{k+1})$用于下一次迭代。

• 并且$\nabla f(x_k)$是上一次迭代位置$x_k$的梯度结果，是已知项。观察上述等式左侧，根据拉格朗日中值定理，可以表示为如下形式：
  由于没有办法确定$x_k,x_{k+1}$之间的大小关系，因而关于$\xi$的描述表示为:$\xi =\lambda \cdot x_k+(1-\lambda )\cdot x_{k+1};\lambda \in (0,1)$而不是$\xi \in (x_k,x_{k+1})$
  $$\nabla f(x_{k+1})-\nabla f(x_k)={\nabla }^{2}f(\xi )\cdot (x_{k+1}-x_k)$$

• 对比拉格朗日中值定理与拟牛顿方程，相当于使用${\mathcal{B}}_{k+1}$对${\nabla }^{2}f(\xi )$进行替换，并且拟牛顿方程内，除了${\mathcal{B}}_{k+1}$，其余项均是已知项。所以${\mathcal{B}}_{k+1}$可求。

  继续观察：关于矩阵$[{\mathcal{B}}_{k+1}{]}_{n\times n}$，首先它必然是对称矩阵，从而包含$\begin{array}{r}\frac{n(n+1)}{2}\end{array}$个变量(上三角/下三角阵元素数量)；但拟牛顿方程所描述的方程组内仅包含$n$个方程($x_k,x_{k+1}\in {\mathbb{R}}^n$)，由于$\begin{array}{r}\frac{n(n+1)}{2}\ge n;n\in {\mathbb{N}}^{+}\end{array}$恒成立，从而满足拟牛顿方程的${\mathcal{B}}_{k+1}$不唯一。

• 为表达方便，记：
  它们都是拟牛顿方程~
  $$\{\begin{array}{l}{y}_k={\mathcal{B}}_{k+1}{\mathcal{S}}_k\\ {\mathcal{S}}_k={\mathcal{H}}_{k+1}{y}_k\end{array}\Leftarrow \{\begin{array}{l}{y}_k=\nabla f(x_{k+1})-\nabla f(x_k)\\ {\mathcal{S}}_k=x_{k+1}-x_k\\ {\mathcal{H}}_k={\mathcal{B}}_k^{-1}\end{array}$$

矩阵${\mathcal{B}}_{k+1}$的选择

选择${\mathcal{B}}_{k+1}$的核心思路：通过已有信息$(y_k,{\mathcal{S}}_k,{\mathcal{B}}_k)\Rightarrow {\mathcal{B}}_{k+1}$或者已有信息$(y_k,{\mathcal{S}}_k,{\mathcal{H}}_k)\Rightarrow {\mathcal{H}}_{k+1}$。
求出那个都可以，因为最终我们需要获取下降方向：${\mathcal{D}}_{k+1}=-{\mathcal{B}}_{k+1}^{-1}\nabla f(x_{k+1})=-{\mathcal{H}}_{k+1}\nabla f(x_{k+1})$

• 方法一：找到满足拟牛顿方程并且与${\mathcal{B}}_k$相似的正定矩阵$\mathcal{B}$作为${\mathcal{B}}_{k+1}$。其数学符号表示如下：
  这里通过对矩阵差异性求解范数来表示近似关系，关于近似关系的表示不唯一。
  $${\mathcal{B}}_{k+1}\Rightarrow \mathcal{B}:\{\begin{array}{l}\min \Vert \mathcal{B}-{\mathcal{B}}_k\Vert \\ \text{s.t. }\mathcal{B}{\mathcal{S}}_k=y_k;{\mathcal{B}}^T=\mathcal{B}\end{array}$$
  同样可以通过上述思想选择矩阵$\mathcal{H}$：
  $${\mathcal{H}}_{k+1}\Rightarrow \mathcal{H}:\{\begin{array}{l}\min \Vert \mathcal{H}-{\mathcal{H}}_k\Vert \\ \text{s.t. }\mathcal{H}{y}_k={\mathcal{S}}_k;{\mathcal{H}}^T=\mathcal{H}\end{array}$$

• 方法二：其思想是：${\mathcal{B}}_{k+1}/{\mathcal{H}}_{k+1}$是基于${\mathcal{B}}_k/{\mathcal{H}}_k$的校正(优化)结果。令${\mathcal{B}}_{k+1}={\mathcal{B}}_k+\Delta \mathcal{B}$：
  无论是$\text{Rank-1}$还是$\text{Rank-2}$校正，其朴素思想是迭代过程中，避免关于${\mathcal{B}}_{k+1}$的复杂运算。

  同理，关于${\mathcal{H}}_{k+1}$也可以使用${\mathcal{H}}_{k+1}={\mathcal{H}}_k+\Delta \mathcal{H}$进行描述。

  • $\text{Rank-1}$校正：要求$\Delta \mathcal{B}$的秩为$1$。代表方法有：$\text{SR-1}$方法。
  • $\text{Rank-2}$校正： 要求$\Delta \mathcal{B}$的秩为$2$。代表方法有：$\text{DFP}$方法，$\text{BFGS}$方法。

下一节将具体介绍$\text{DFP}$以及$\text{BFGS}$方法。

$\text{Reference}$：
最优化理论与方法-第六讲-无约束优化问题（二）