引言

上一节介绍了牛顿法、拟牛顿法。本节将继续以拟牛顿法为基础，介绍$\text{DFP},\text{BFGS}$方法。

回顾：经典牛顿法的缺陷与拟牛顿法思想

经典牛顿法缺陷与修正牛顿法

关于经典牛顿法中关于下降方向${\mathcal{D}}_k(k=1,2,\cdots ,\infty )$的数学符号表示如下：
$${\mathcal{D}}_k=-[{\nabla }^{2}f(x_k){]}^{-1}\nabla f(x_k)$$
其中$\nabla f(x_k)$表示目标函数$f(\cdot )$在$x_k$位置的梯度向量结果；${\nabla }^{2}f(x_k)$表示目标函数在$x_k$位置的$\text{Hessian Matrix}$。问题在于：${\nabla }^{2}f(x_k)$可能不是正定矩阵，从而无法求解$[{\nabla }^{2}f(x_k){]}^{-1}$，最终无法执行迭代过程。

关于这类问题，可以使用正则化法对${\nabla }^{2}f(x_k)$进行修正：
$${\nabla }^{2}f(x_k):={\nabla }^{2}f(x_k)+\lambda \mathcal{I}$$
其中$\mathcal{I}$表示单位矩阵。执行该操作的目的是：保持${\nabla }^{2}f(x_k)$是正定矩阵状态。但这种方法同样存在弊端：
λ > max ⁡ i = 1 , 2 , ⋯ , n { − λ i } \lambda > \mathop{\max}\limits_{i=1,2,\cdots,n} \{- \lambda_i\} λ>i=1,2,⋯,nmax​{−λi​}
如果$\lambda$数值过大，可能会发生原始${\nabla }^{2}f(x_k)$中各特征值被$\lambda$分掉相应权重，从而导致修正后的${\nabla }^{2}f(x_k)$中关于$x_k$的二阶梯度信息减少，甚至无效。当然，也可以基于正则化法的思想，对${\nabla }^{2}f(x_k)$进行优化：
实际上，正则化法中$\lambda$过大最终影响当前迭代步骤的下降方向,并使其收敛到$\begin{array}{r}\frac{\nabla f(x_k)}{\lambda }\end{array}$。
$$\begin{array}{rl}& {\nabla }^{2}f(x_k)={\mathcal{Q}}^T\text{Diag}({\tau }_i)\mathcal{Q}\\ & {\tau }_i=\{\begin{array}{l}{\tau }_i\text{if }{\tau }_i\ge \delta \\ \delta \text{Otherwise}\end{array}\end{array}$$
其中$\delta$是一个适当正数；虽然该方式相比正则化法要缓和不少——仅调整非正特征值的结果，其余正特征值保持不变。但该方法依然存在逻辑上的缺失：通过强行修改二阶梯度信息的方式使其收敛。

拟牛顿法与矩阵${\mathcal{B}}_{k+1}$的选择

而拟牛顿法的思想是：选择一个既包含$x_{k+1}$处的二阶梯度信息，并且容易获取的正定矩阵${\mathcal{B}}_{k+1}$来替代${\nabla }^{2}f(x_{k+1})$。
由于$[{\nabla }^{2}f(x_{k+1}){]}_{n\times n}$自身计算量较大:$\mathcal{O}(n^3)$,从而不容易获取。

关于矩阵${\mathcal{B}}_{k+1}$满足的基本要求表示如下：
$$\nabla f(x_{k+1})-\nabla f(x_k)={\mathcal{B}}_{k+1}(x_{k+1}-x_k)$$
可以发现：该式子是关于$n$个方程构成的方程组；而未知量包含$\begin{array}{r}\frac{n(n+1)}{2}\end{array}$个(${\mathcal{B}}_{k+1}$上/下三角阵元素数量)，并且：$\begin{array}{r}\frac{n(n+1)}{2}\ge n;n\in {\mathbb{N}}^{+}\end{array}$。这意味着拟牛顿方程的解${\mathcal{B}}_{k+1}$不唯一。

既然满足基本要求的解不唯一，可以尝试从这些解中选择与${\mathcal{B}}_k/{\mathcal{H}}_k$相似的矩阵作为${\mathcal{B}}_{k+1}/{\mathcal{H}}_{k+1}$：

• 其中:$\{\begin{array}{l}{\mathcal{S}}_k=x_{k+1}-x_k\\ y_k=\nabla f(x_{k+1})-\nabla f(x_k)\\ {\mathcal{H}}_k={\mathcal{B}}_k^{-1}\end{array}$
• 通过这种相似性来保证二阶梯度信息的有效性。
• 无论是${\mathcal{B}}_{k+1}$还是${\mathcal{H}}_{k+1}$都可以作为我们的求解目标。因为最终都可以对下降方向${\mathcal{D}}_{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})$。
  $$\{\begin{array}{l}{\mathcal{B}}_{k+1}\Rightarrow \mathcal{B}:\{\begin{array}{l}\min \Vert \mathcal{B}-{\mathcal{B}}_k\Vert \\ \text{s.t. }\mathcal{B}\cdot {\mathcal{S}}_k=y_k;{\mathcal{B}}^T=\mathcal{B}\end{array}\\ \\ {\mathcal{H}}_{k+1}\Rightarrow \mathcal{H}:\{\begin{array}{l}\min \Vert \mathcal{H}-{\mathcal{H}}_k\Vert \\ \text{s.t. }\mathcal{H}\cdot y_k={\mathcal{S}}_k;{\mathcal{H}}^T=\mathcal{H}\end{array}\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}$或者${\mathcal{H}}_{k+1}={\mathcal{H}}_k+\Delta \mathcal{H}$，其中：

• $\text{SR-1}$方法是$\text{Rank-1}$校正：$\Delta \mathcal{B}、\Delta \mathcal{H}$的秩为$1$的代表方法；
• $\text{DFP,BFGS}$方法是$\text{Rank-2}$校正：$\Delta \mathcal{B}、\Delta \mathcal{H}$的秩为$2$的代表方法。

拟牛顿法之$\text{DFP}$方法

关于$\text{DFP(Davidon-Fletcher-Power)}$方法可看做是对${\mathcal{H}}_k$进行$\text{Rank-2}$校正。对应迭代公式表示如下：
$${\mathcal{H}}_{k+1}={\mathcal{H}}_k-\frac{{\mathcal{H}}_k{y}_k{y}_k^T{\mathcal{H}}_k}{y_k^T{\mathcal{H}}_k{y}_k}+\frac{{\mathcal{S}}_k{\mathcal{S}}_k^T}{y_k^T{\mathcal{S}}_k}$$

$\text{DFP}$迭代公式的推导过程

$\text{DFP}$是一个$\text{Rank-2}$校正方法，那么如何表示一个秩为$2$的矩阵$?$首先，先观察秩为$1$的矩阵如何表示：某矩阵${\mathcal{A}}_{n\times n}$可表示为如下形式：
该矩阵的所有行均相同。
$$\mathcal{A}=\mathcal{U}{\mathcal{V}}^T\mathcal{U},\mathcal{V}\in {\mathbb{R}}^n;\mathcal{U},\mathcal{V}\ne 0$$
此时$\mathcal{A}$就是一个秩为$1$的矩阵。但由于${\mathcal{H}}_k$必然是一个对称矩阵，相比于上式，$\Delta \mathcal{H}$想满足是秩为$1$仅需要满足：
$$\Delta \mathcal{H}=\mathcal{U}{\mathcal{U}}^T\mathcal{U}\in {\mathbb{R}}^n;\mathcal{U}\ne 0$$
这是秩为$1$的情况。那秩为$2$呢$?$只需要满足：
$$\Delta \mathcal{H}=\mathcal{U}{\mathcal{U}}^T+\mathcal{V}{\mathcal{V}}^T\{\begin{array}{l}\mathcal{U},\mathcal{V}\in {\mathbb{R}}^n\\ \mathcal{U},\mathcal{V}\ne 0\\ \mathcal{U}\ne \mathcal{V}\end{array}$$
综上，将迭代关系：${\mathcal{H}}_{k+1}={\mathcal{H}}_k+\Delta \mathcal{H}$表示为如下形式：
其中$a,b$是系数，均是标量~
$${\mathcal{H}}_{k+1}={\mathcal{H}}_k+a\cdot \mathcal{U}{\mathcal{U}}^T+b\cdot \mathcal{V}{\mathcal{V}}^T$$
由于${\mathcal{H}}_{k+1}$需要满足基本要求：${\mathcal{H}}_{k+1}\cdot y_k={\mathcal{S}}_k$，因而将上式带入。有：
$${\mathcal{H}}_k{y}_k+a\cdot \mathcal{U}{\mathcal{U}}^T{y}_k+b\cdot \mathcal{V}{\mathcal{V}}^T{y}_k-{\mathcal{S}}_k=0$$
其中：

• 由于${\mathcal{H}}_k\in {\mathbb{R}}^{n\times n},y_k\in {\mathbb{R}}^n$，因而${\mathcal{H}}_k{y}_k\in {\mathbb{R}}^n$，是一个$n$维向量；
• 由于${\mathcal{U}}^T{y}_k\in \mathbb{R}$，因而$a\cdot \mathcal{U}({\mathcal{U}}^T{y}_k)\in {\mathbb{R}}^n$，可看做向量$\mathcal{U}$的$a\cdot ({\mathcal{U}}^T{y}_k)$倍；
• 同理，$b\cdot \mathcal{V}({\mathcal{V}}^T{y}_k)\in {\mathbb{R}}^n$，可看做向量$\mathcal{V}$的$b\cdot ({\mathcal{V}}^T{y}_k)$倍。
• ${\mathcal{S}}_k=x_{k+1}-x_k\in {\mathbb{R}}^n$

对$\mathcal{U},\mathcal{V}$进行取值。将项${\mathcal{H}}_k,a\cdot \mathcal{U}({\mathcal{U}}^T{y}_k)$关联在一起；项$b\cdot \mathcal{V}({\mathcal{V}}^T{y}_k),{\mathcal{S}}_k$关联在一起：
$$\underset{=0}{\underbrace{\left[{\mathcal{H}}_k{y}_k+a\cdot \mathcal{U}({\mathcal{U}}^T{y}_k)\right]}}+\underset{=0}{\underbrace{\left[b\cdot \mathcal{V}({\mathcal{V}}^T{y}_k)-{\mathcal{S}}_k\right]}}=0$$
观察第一项：令$\mathcal{U}={\mathcal{H}}_k{y}_k$，带入有：
$$\begin{array}{rl}{\mathcal{H}}_k{y}_k+a\cdot \mathcal{U}({\mathcal{U}}^T{y}_k)& ={\mathcal{H}}_k{y}_k+a\cdot {\mathcal{H}}_k{y}_k[({\mathcal{H}}_k{y}_k{)}^T{y}_k]\\ & =({\mathcal{H}}_k{y}_k)[1+a\cdot ({\mathcal{H}}_k{y}_k{)}^T{y}_k]\\ & =0\\ & \Rightarrow 1+a\cdot ({\mathcal{H}}_k{y}_k{)}^T{y}_k=0\end{array}$$
整理得：$\begin{array}{r}a=-\frac{1}{y_k^T{\mathcal{H}}_k^T{y}_k}\end{array}$。
同理，观察第二项：令$\mathcal{V}={\mathcal{S}}_k$，带入有：
$$b\cdot {\mathcal{S}}_k^T{y}_k-1=0\Rightarrow b=\frac{1}{{\mathcal{S}}_k^T{y}_k}$$
至此，关于向量$\mathcal{U},\mathcal{V}$，系数$a,b$均已取值完毕，将该结果带入${\mathcal{H}}_{k+1}={\mathcal{H}}_k+a\cdot \mathcal{U}{\mathcal{U}}^T+b\cdot \mathcal{V}{\mathcal{V}}^T$，即可得到$\text{DFP}$公式中${\mathcal{H}}_{k+1}$与${\mathcal{H}}_k$之间的迭代关系。

小插曲：$\text{DFP}$方法与最小范数方法

关于最小范数方法：${\mathcal{B}}_{k+1}\Rightarrow \mathcal{B}:\{\begin{array}{l}\min \Vert \mathcal{B}-{\mathcal{B}}_k\Vert \\ \text{s.t. }\mathcal{B}\cdot {\mathcal{S}}_k=y_k;{\mathcal{B}}^T=\mathcal{B}\end{array}$，如果使用$\text{Frobenius}$范数对$\Vert \mathcal{B}-{\mathcal{B}}_k\Vert$进行表示：
可以看成是关于矩阵的$L_2$范数。
$$\Vert \mathcal{B}-{\mathcal{B}}_k{\Vert }_F=\sqrt{\sum _{i=1}^n\sum _{j=1}^n{\left[b_{ij}-b_{ij}^{(k)}\right]}^2}$$
通过该范数求解出的${\mathcal{B}}_{k+1}$，它的逆：${\mathcal{B}}_{k+1}^{-1}$就是$\text{DFP}$方法求解出的${\mathcal{H}}_{k+1}$。
世界真奇妙~

拟牛顿法之$\text{BFGS}$方法

关于$\text{BFGS(Broyden-Fletch-Goldfarb-Shannon)}$方法可看做是对${\mathcal{B}}_k$进行$\text{Rank-2}$校正。对应迭代公式表示如下：
$${\mathcal{B}}_{k+1}={\mathcal{B}}_k-\frac{{\mathcal{B}}_k{\mathcal{S}}_k{\mathcal{S}}_k^T{\mathcal{B}}_k}{{\mathcal{S}}_k^T{\mathcal{B}}_k{\mathcal{S}}_k}+\frac{y_k{y}_k^T}{y_k^T{\mathcal{S}}_k}$$

关于$\text{BFGS}$公式的推导，它与$\text{DFP}$公式的推导完全对称。只不过它使用的基本要求是：${\mathcal{B}}_{k+1}\cdot {\mathcal{S}}_k=y_k$。

• 对比$\text{DFP}$公式：仅需要将第一项中的$y_k$改成${\mathcal{S}}_k$，${\mathcal{H}}_k$改成${\mathcal{B}}_k$；第二项将分子中的${\mathcal{S}}_k$改成$y_k$即可。
• 关于$\text{BFGS}$公式的推导不再赘述。

新的疑问：在使用$\text{BFGS}$求解出${\mathcal{B}}_{k+1}$后，在后续求解下降方向${\mathcal{D}}_k=-{\mathcal{B}}_{k+1}^{-1}\nabla f(x_{k+1})$中，依然不可避免地需要求解逆：${\mathcal{B}}_{k+1}^{-1}$。而求逆同样是一个非常麻烦的操作，为什么还会使用$\text{BFGS}$方法$?$主要有两点原因：

• 具备${\mathcal{B}}_{k+1}={\mathcal{B}}_k+a\cdot \mathcal{U}{\mathcal{U}}^T+b\cdot \mathcal{V}{\mathcal{V}}^T$格式的逆可以使用$\text{Sherman-Morrison}$公式直接求解：
  可以看出，求逆操作自身并不麻烦。
  $$(\mathcal{A}+\mathcal{U}{\mathcal{V}}^T{)}^{-1}={\mathcal{A}}^{-1}-\frac{{\mathcal{A}}^{-1}\mathcal{U}{\mathcal{V}}^T{\mathcal{A}}^{-1}}{1+{\mathcal{V}}^T{\mathcal{A}}^{-1}\mathcal{U}}$$
• $\text{DFP}$方法求解，其结果稳定性较差。在迭代过程中可能出现${\mathcal{H}}_{k+1}$变成奇异矩阵。相反，$\text{BFGS}$迭代过程中的数值稳定性更强。并且$\text{BFGS}$被认为是最有效的拟牛顿法，它的收敛速度可达到超线性收敛。

相比于牛顿法中直接求解$\text{Hessian Matrix}\Rightarrow {\nabla }^{2}f(x_k)$，$\text{DFP,BFGS}$方法需要求解梯度$\nabla f(x_k),\nabla f(x_{k+1})$，以及套用求逆公式。其计算量远小于求解$\text{Hessian Matrix}$。

$\text{Broyden}$族

假设使用$\text{DFP}$方法求解出${\mathcal{H}}_{k+1}$，将该结果求逆，将其还原：
$${\mathcal{B}}_{\text{DFP};k+1}={\mathcal{H}}_{k+1}^{-1}$$
然后通过$\text{BFGS}$方法直接求解出${\mathcal{B}}_{k+1}$。对这两个矩阵进行线性组合：
{ λ ⋅ B DFP ; k + 1 + ( 1 − λ ) ⋅ B k + 1 } λ ∈ [ 0 , 1 ] \{\lambda \cdot \mathcal B_{\text{DFP};k+1} + (1 - \lambda) \cdot \mathcal B_{k+1}\} \quad \lambda \in [0,1] {λ⋅BDFP;k+1​+(1−λ)⋅Bk+1​}λ∈[0,1]
这明显是一个集合。如果迭代过程中，矩阵${\mathcal{B}}_{k+1}$落在集合内，对应的方法被称作$\text{Broyden}$族。

拟牛顿法之$\text{SR-1}$方法

关于$\text{SR-1}$方法可看做是对${\mathcal{B}}_k$进行$\text{Rank-1}$校正。对应迭代公式表示如下：
$${\mathcal{B}}_{k+1}={\mathcal{B}}_k+\frac{(y_k-{\mathcal{B}}_k{\mathcal{S}}_k)(y_k-{\mathcal{B}}_k{\mathcal{S}}_k{)}^T}{(y_k-{\mathcal{B}}_k{\mathcal{S}}_k{)}^T{\mathcal{S}}_k}$$

$\text{SR-1}$迭代公式的推导过程

与$\text{DFP}$方法的推导过程类似。将迭代关系：${\mathcal{B}}_{k+1}={\mathcal{B}}_k+\Delta \mathcal{B}$表示为如下形式：
$${\mathcal{B}}_{k+1}={\mathcal{B}}_k+a\cdot \mathcal{U}{\mathcal{U}}^T$$
由于${\mathcal{B}}_{k+1}$需要满足基本要求：${\mathcal{B}}_{k+1}\cdot {\mathcal{S}}_k=y_k$。因而将上式带入，有：
$${\mathcal{B}}_k{\mathcal{S}}_k+a\cdot \mathcal{U}({\mathcal{U}}^T{\mathcal{S}}_k)=y_k\Rightarrow a\cdot \mathcal{U}({\mathcal{U}}^T{\mathcal{S}}_k)=y_k-{\mathcal{B}}_k{\mathcal{S}}_k$$
令$\mathcal{U}=y_k-{\mathcal{B}}_k{\mathcal{S}}_k$，有：系数$a\cdot ({\mathcal{U}}^T{\mathcal{S}}_k)=1$，最终可求出$a$：
$$a=\frac{1}{{\mathcal{U}}^T{\mathcal{S}}_k}=\frac{1}{(y_k-{\mathcal{B}}_k{\mathcal{S}}_k{)}^T{\mathcal{S}}_k}$$
将$a,\mathcal{U}$带回${\mathcal{B}}_{k+1}={\mathcal{B}}_k+a\cdot \mathcal{U}{\mathcal{U}}^T$，就有$\text{SR-1}$迭代公式。

不可否认：$\text{SR-1}$方法的迭代公式更加简便，但它不能保证迭代过程中${\mathcal{B}}_{k+1}$的正定性。在适当条件下，$\text{SR-1}$算法可达到$n$步超线性收敛。
这里的$n$步超线性收敛是指：当前步骤与执行$n$步之后的结果呈超线性收敛趋势。对比超线性收敛，其数学符号表示如下：
$$\{\begin{array}{l}\begin{array}{r}\underset{k\to \infty }{\lim }\frac{\Vert x_{k+1}-x^{\ast }\Vert }{\Vert x_k-x^{\ast }\Vert }=0\\ \underset{k\Rightarrow \infty }{\lim }\frac{\Vert x_{k+n}-x^{\ast }\Vert }{\Vert x_k-x^{\ast }\Vert }=0\end{array}\end{array}$$

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