拟牛顿法BFGS的一些修正公式

牛顿法 (Newton Method)
牛顿法的基本思想是在极小点附近通过对目标函数$f(x)$做二阶Taylor展开，进而找到$f(x)$的极小点的估计值[1]。一维情况下，也即令函$\phi (x)$为

$\phi (x)=f(x_k)+f^{{}^{\prime }}(x_k)(x-x_k)+\frac{1}{2}{f}^{{}^{\prime \prime }}(x_k)(x-x_k{)}^2$

则其导数${\phi }^{{}^{\prime }}(x)$满足

${\phi }^{{}^{\prime }}(x)=f^{{}^{\prime }}(x_k)+f^{{}^{\prime \prime }}(x_k)(x-x_k)=0$

因此

$x_{k+1}=x_k-\frac{f^{{}^{\prime }}(x_k)}{f^{{}^{\prime \prime }}(x_k)}$ (1)

将$x_{k+1}$作为$f(x)$极小点的一个进一步的估计值。重复上述过程，可以产生一系列的极小点估值集合 { x k } \{x_k\} {xk​}。一定条件下，这个极小点序列 { x k } \{x_k\} {xk​}收敛于$f(x)$的极值点。

将上述讨论扩展到$N$维空间，类似的，对于$N$维函数$f(\mathbf{x})$有

$f(\mathbf{x})\approx \phi (\mathbf{x})=f({\mathbf{x}}_k)+\nabla f({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^{\mathrm{T}}{\nabla }^{2}f(\mathbf{x}-{\mathbf{x}}_k)$

其中$\nabla f(\mathbf{x})$和${\nabla }^{2}f(\mathbf{x})$分别是目标函数的的一阶和二阶导数，表现为$N$维向量和$N$ $\times$ $N$矩阵，而后者又称为目标函数$f(\mathbf{x})$在${\mathbf{x}}_k$处的Hesse矩阵。 设${\nabla }^{2}f(\mathbf{x})$可逆，则可得与方程(1)类似的迭代公式：

${\mathbf{x}}_{k+1}={\mathbf{x}}_k-[{\nabla }^{2}f({\mathbf{x}}_k{]}^{-1}\nabla f({\mathbf{x}}_k)$ (2)

这就是原始牛顿法的迭代公式。

原始牛顿法虽然具有二次终止性（即用于二次凸函数时，经有限次迭代必达极小点），但是要求初始点需要尽量靠近极小点，否则有可能不收敛。因此人们又提出了阻尼牛顿法[1]。这种方法在算法形式上等同于所有流行的优化方法，即确定搜索方向，再沿此方向进行一维搜索，找出该方向上的极小点，然后在该点处重新确定搜索方向，重复上述过程，直至函数梯度小于预设判据$ϵ$。具体步骤列为算法1。

算法1：

(1) 给定初始点${\mathbf{x}}_0$，设定收敛判据$ϵ$，$k=0$.

(2) 计算$\nabla f({\mathbf{x}}_k)$和${\nabla }^{2}f({\text{x}}_k)$.

(3) 若$||\nabla f({\mathbf{x}}_k)||<ϵ$，则停止迭代，否则确定搜索方向${\mathbf{d}}_k=-[{\nabla }^{2}f({\mathbf{x}}_k){]}^{-1}\nabla f({\mathbf{x}}_k)$.

(4) 从${\mathbf{x}}_k$出发，沿${\mathbf{d}}_k$做一维搜索，

$\underset{\lambda }{\min }f({\mathbf{x}}_k+\lambda {\mathbf{d}}_k)=f({\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k)$

令${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k$.

(5) 设$k=k+1$，转步骤(2).

在一定程度上，阻尼牛顿法具有更强的稳定性

拟牛顿法

如同上一节指出，牛顿法虽然收敛速度快，但是计算过程中需要计算目标函数的二阶偏导数，难度较大。更为复杂的是目标函数的Hesse矩阵无法保持正定，从而令牛顿法失效。为了解决这两个问题，人们提出了拟牛顿法。这个方法的基本思想是不用二阶偏导数而构造出可以近似Hesse矩阵的逆的正定对称阵， 从而在"拟牛顿"的条件下优化目标函数。构造方法的不同决定了不同的拟牛顿法。

首先分析如何构造矩阵可以近似Hesse矩阵的逆：

设第k次迭代之后得到点${\mathbf{x}}_{k+1}$，将目标函数$f(\mathbf{x})$在${\mathbf{x}}_{k+1}$处展成Taylor级数，取二阶近似，得到

$f(\mathbf{x})\approx f({\mathbf{x}}_{k+1})+\nabla f({\mathbf{x}}_{k+1})(\mathbf{x}-{\mathbf{x}}_{k+1})+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_{k+1}{)}^{\mathrm{T}}{\nabla }^{2}f({\mathbf{x}}_{k+1})(\mathbf{x}-{\mathbf{x}}_{k+1})$

因此

$\nabla f(\mathbf{x})\approx \nabla f({\mathbf{x}}_{k+1})+{\nabla }^{2}f({\mathbf{x}}_{k+1})(\mathbf{x}-{\mathbf{x}}_{k+1})$

令$\mathbf{x}={\mathbf{x}}_k$，则

$\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)\approx {\nabla }^{2}f({\mathbf{x}}_{k+1})({\mathbf{x}}_k-{\mathbf{x}}_{k+1})$ (3)

记

${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k,{\mathbf{y}}_k=\nabla f({\mathbf{x}}_{k+1})-\nabla f({\mathbf{x}}_k)$

同时设Hesse矩阵${\nabla }^{2}f({\mathbf{x}}_{k+1})$可逆，则方程(3)可以表示为

${\mathbf{s}}_k\approx [{\nabla }^{2}f({\mathbf{x}}_{k+1}){]}^{-1}{\mathbf{y}}_k$ (4)

因此，只需计算目标函数的一阶导数，就可以依据方程(4)估计该处的Hesse矩阵的逆。也即，为了用不包含二阶导数的矩阵${\mathbf{H}}_{k+1}$近似牛顿法中的Hesse矩阵${\nabla }^{2}f({\mathbf{x}}_{k+1})$的逆矩 阵，${\mathbf{H}}_{k+1}$必须满足

${\mathbf{s}}_k\approx {\mathbf{H}}_{k+1}{\mathbf{y}}_k$ (5)

方程(5)也称为拟牛顿条件。不加证明的，下面给出两个最常用的${\mathbf{H}}_{k+1}$构造公式

BFGS公式

BFGS公式有时也称为DFP公式的对偶公式。这是因为其推导过程与方程(6)完全一样，只需要用矩阵${\mathbf{B}}_{k+1}$取代${\mathbf{H}}_{k+1}^{-1}$，同时将${\mathbf{s}}_k$和 ${\mathbf{y}}_k$互换，最后可以得到

${\mathbf{H}}_{k+1}={\mathbf{H}}_k+[1+\frac{{\mathbf{y}}_k^{\mathrm{T}}{\mathbf{H}}_k{\mathbf{y}}_k}{{\mathbf{s}}_k^{\mathrm{T}}{\mathbf{y}}_k}]\cdot \frac{{\mathbf{s}}_k{\mathbf{s}}_k^{\mathrm{T}}}{{\mathbf{s}}_k^{\mathrm{T}}{\mathbf{y}}_k}-\frac{{\mathbf{s}}_k{\mathbf{y}}_k^{\mathrm{T}}{\mathbf{H}}_k}{{\mathbf{s}}_k^{\mathrm{T}}{\mathbf{y}}_k}$ (7)

这个公式要优于DFP公式，因此目前得到了最为广泛的应用。

利用方程(6)或(7)的拟牛顿法计算步骤列为算法2。

算法2：

(1) 给定初始点${\mathbf{x}}_0$，设定收敛判据$ϵ$，$k=0$.

(2) 设${\mathbf{H}}_0=\mathbf{I}$，计算出目标函数$f(\mathbf{x})$在${\mathbf{x}}_k$处的梯度$g_k=\nabla f({\mathbf{x}}_k)$.

(3) 确定搜索方向${\mathbf{d}}_k$：

${\mathbf{d}}_k={\mathbf{H}}_k{\mathbf{g}}_k$.

(4) 从${\mathbf{x}}_k$出发，沿${\mathbf{d}}_k$做一维搜索，满足

$f({\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k)=\underset{\lambda \ge 0}{\min }f({\mathbf{x}}_k+\lambda {\mathbf{d}}_k)$

令${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k$.

(5) 若$||f({\mathbf{x}}_{k+1})||\le ϵ$，则停止迭代，得到最优解$\mathbf{x}={\mathbf{x}}_{k+1}$，否则进行步骤(6).

(6) 若$k=N-1$，则令${\mathbf{x}}_0={\mathbf{x}}_{k+1}$，回到步骤(2)，否则进行步骤(7).

(7) 令${\mathbf{g}}_{k+1}=f^{{}^{\prime }}({\mathbf{x}}_{k+1})$，${\mathbf{s}}_k={\mathbf{x}}_{k+1}-{\mathbf{x}}_k$，${\mathbf{y}}_k={\mathbf{g}}_{k+1}-{\mathbf{g}}_k$，利用方程(6)或(7)计算${\mathbf{H}}_{k+1}$，设$k=k+1$，回到步骤(3)。

对于二次凸函数，BFGS算法具有良好的全局收敛性。但是对于二次非凸函数，也即目标函数Hesse矩阵非正定的情况，无法保证按照BFGS算法构造的拟牛顿方向必为下降方向。为了推广BFGS公式的应用范围，很多工作提出了对BFGS公式稍作修改或变形的办法。下面举两个例子。

Li-Fukushima方法[3]
Li和Fukushima提出新的构造矩阵${\mathbf{H}}_k$的方法：

${\mathbf{H}}_{k+1}^{-1}={\mathbf{H}}_k^{-1}-\frac{{\mathbf{H}}_k^{-1}{\mathbf{s}}_k{\mathbf{s}}_k^{\mathrm{T}}{\mathbf{H}}_k^{-1}}{{\mathbf{s}}_k^{\mathrm{T}}{\mathbf{H}}_k^{-1}{\mathbf{s}}_k}+\frac{{\mathbf{y}}_k^{\ast }{\mathbf{y}}_k^{\ast \mathrm{T}}}{{\mathbf{y}}_k^{\ast \mathrm{T}}{\mathbf{s}}_k}$

${\mathbf{H}}_{k+1}=(\mathbf{I}-{\rho }_k^{\ast }{\mathbf{s}}_k{\mathbf{y}}_k^{\ast \mathrm{T}}){\mathbf{H}}_k(\mathbf{I}-{\rho }_k^{\ast }{\mathbf{y}}_k^{\ast \mathrm{T}}{\mathbf{s}}_k^{\mathrm{T}})+{\rho }_k^{\ast }{\mathbf{s}}_k{\mathbf{s}}_k^{\mathrm{T}}$ (11)

其中

${\mathbf{y}}_k^{\ast }={\mathbf{g}}_{k+1}-{\mathbf{g}}_k+t_k||{\mathbf{g}}_k||{\mathbf{s}}_k$
t k = 1 + max ⁡ { 0 , − y k T s k ∣ ∣ s k ∣ ∣ 2 } t_k=1+\max\{0, \frac{-\mathbf{y}^{\rm T}_k\mathbf{s}_k}{||\mathbf{s}_k||^2}\} tk​=1+max{0,∣∣sk​∣∣2−ykT​sk​​}

${\mathbf{y}}_k$的定义见算法2中步骤(7)，而

${\rho }_k^{\ast }=\frac{1}{{\mathbf{y}}_k^{\ast \mathrm{T}}{\mathbf{s}}_k}$

除此之外，算法2中步骤(3)的一维搜索采用如下方式：

给定两个参数$\sigma \in (0,1)$和$ϵ\in (0,1)$，找出最小的非负整数j，满足

$f({\mathbf{x}}_k+{ϵ}_j{\mathbf{d}}_k)\le f({\mathbf{x}}_k)+\sigma {ϵ}_j{\mathbf{g}}_k^{\mathrm{T}}{\mathbf{d}}_k$

取$j_k=j$，步长${\lambda }_k={ϵ}_{j_k}$。

Xiao-Wei-Wang方法[4]

Xiao、Wei和Wang提出了计入目标函数值$f(\mathbf{x})$的另一种${\mathbf{H}}_k$的构造方法：

设${\mathbf{y}}_k^{†}={\mathbf{y}}_k+{\alpha }_k{\mathbf{s}}_k$，其中

${\alpha }_k=\frac{1}{||{\mathbf{s}}_k|{|}^2}[2(f({\mathbf{x}}_k)-f({\mathbf{x}}_{k+1})+({\mathbf{g}}_{k+1}+{\mathbf{g}}_k{)}^{\mathrm{T}}{\mathbf{s}}_k]$

${\mathbf{H}}_k$的构造方法与方程(7)和(11)形式相同：

${\mathbf{H}}_{k+1}=(\mathbf{I}-{\rho }_k^{†}{\mathbf{s}}_k{\mathbf{y}}_k^{†\mathrm{T}}){\mathbf{H}}_k(\mathbf{I}-{\rho }_k^{†}{\mathbf{y}}_k^{†}){\mathbf{s}}_k^{\mathrm{T}}+{\rho }_k^{†}{\mathbf{s}}_k{\mathbf{s}}_k^{\mathrm{T}}$ (12)
相应的${\rho }_k^{†}=\frac{1}{{\mathbf{y}}_k^{†\mathrm{T}}{\mathbf{s}}_k}$
而一维搜索则采用弱Wolfe-Powell准则：

给定两个参数$\delta \in (0,1/2)$和$\sigma \in (\delta ,1)$，找出步长${\lambda }_k$，满足

$f({\mathbf{x}}_k+{\lambda }_k{\mathbf{d}}_k)\le f({\mathbf{x}}_k)+\delta {\lambda }_k{\mathbf{g}}_k^{\mathrm{T}}{\mathbf{d}}_k$ (13)

${\mathbf{g}}_{k+1}^{\mathrm{T}}{\mathbf{d}}_k\ge \sigma {\mathbf{g}}_k^{\mathrm{T}}{\mathbf{d}}_k$ (14)

如果${\lambda }_k$ = $1$满足方程(13)、(14)，则取${\lambda }_k$ = $1$。

可以看出，这两种方法只是改变了${\mathbf{y}}_k$的定义方式，其他则与标准的BFGS方法完全一样。因此将二者推广到限域形式是非常直接的，这里不再给出算法。对于二次非凸函数的拟牛顿方法还在进一步发展当中，上述的两个例子并不一定是最佳算法。