拟牛顿法及其相关解法 -- 缪缪人生-- 编程爱好者博客

最近在看条件随机场中的优化算法。其中就设计到了无约束化的最优化方法，也就是牛顿法。在CRF（conditional random field）中，使用的是L-BFGS法。费了好大的劲把算法的原理及推导算是看明白了，可是到了具体实现上，又碰到问题了，比如在求搜索方向的时候，使用

$\mathbf{H}_{k+1}=(\mathbf{I}-\rho_k\mathbf{s}_k\mathbf{y}^{\rm T}_k)\mathbf{H}_k(\mathbf{I}-\rho_k\mathbf{y}_k\mathbf{s}^{\rm T}_k)+\rho_k\mathbf{s}_k\mathbf{s}^{\rm T}_k=\mathbf{v}^{\rm T}_k\mathbf{H}_k\mathbf{v}_k+\rho_k\mathbf{s}_k\mathbf{s}^{\rm T}_k$

但是程序中如何实现呢？

现在转载一篇文章，看过之后，会非常受益。

使用导数的最优化算法中，拟牛顿法是目前为止最为行之有效的一种算法，具有收敛速度快、算法稳定性强、编写程序容易等优点。在现今的大型计算程序中有着广泛的应用。本文试图介绍拟牛顿法的基础理论和若干进展。

牛顿法 (Newton Method)

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

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

则其导数 $\varphi^{'}(x)$ 满足

$\varphi^{'}(x) =f^{'}(x_k)+f^{''}(x_k)(x-x_k)=0$

因此

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

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

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

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

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

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

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

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

算法1：

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

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

(3) 若 $||\nabla f(\mathbf{x}_k)||$ < $\epsilon$ ，则停止迭代，否则确定搜索方向 $\mathbf{d}_k=-[\nabla^2f(\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).

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

拟牛顿法 (Quasi-Newton Method)

如同上一节指出，牛顿法虽然收敛速度快，但是计算过程中需要计算目标函数的二阶偏导数，难度较大。更为复杂的是目标函数的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})^{\rm T}\nabla^2f(\mathbf{x}_{k+1})(\mathbf{x}-\mathbf{x}_{k+1})$

因此

$\nabla f(\mathbf{x}) \approx \nabla f(\mathbf{x}_{k+1})+\nabla^2f(\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^2f(\mathbf{x}_{k+1})(\mathbf{x}_k-\mathbf{x}_{k+1})$ (3)

记

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

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

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

因此，只需计算目标函数的一阶导数，就可以依据方程(4)估计该处的Hesse矩阵的逆。也即，为了用不包含二阶导数的矩阵 $\mathbf{H}_{k+1}$ 近似牛顿法中的Hesse矩阵 $\nabla^2f(\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}$ 构造公式

DFP公式

设初始的矩阵 $\mathbf{H}_0$ 为单位矩阵 $\textbf{I}$ ，然后通过修正 $\mathbf{H}_k$ 给出 $\mathbf{H}_{k+1}$ ，即

$\mathbf{H}_{k+1}=\mathbf{H}_k+\Delta \mathbf{H}_k$

DFP算法中定义校正矩阵为

$\Delta \mathbf{H}_k=\frac{\mathbf{s}_k\mathbf{s}_k^{\rm T}}{\mathbf{s}_k^{\rm T}\mathbf{y}_k}-\frac{\mathbf{H}_k\mathbf{y}_k\mathbf{y}_k^{\rm T}\mathbf{H}_k}{\mathbf{y}_k^{\rm T}\mathbf{H}_k\mathbf{y}_k}$

因此

$\mathbf{H}_{k+1}=\mathbf{H}_k+\frac{\mathbf{s}_k\mathbf{s}_k^{\rm T}}{\mathbf{s}_k^{\rm T}\mathbf{y}_k}-\frac{\mathbf{H}_k\mathbf{y}_k\mathbf{y}_k^{\rm T}\mathbf{H}_k}{\mathbf{y}_k^{\rm T}\mathbf{H}_k\mathbf{y}_k}$ (6)

可以验证，这样产生的 $\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^{\rm T}\mathbf{H}_k\mathbf{y}_k}{\mathbf{s}_k^{\rm T}\mathbf{y}_k}]\cdot\frac{\mathbf{s}_k\mathbf{s}_k^{\rm T}}{\mathbf{s}_k^{\rm T}\mathbf{y}_k}-\frac{\mathbf{s}_k\mathbf{y}_k^{\rm T}\mathbf{H}_k}{\mathbf{s}_k^{\rm T}\mathbf{y}_k}$ (7)

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

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

算法2：

(1) 给定初始点 $\mathbf{x}_0$ ，设定收敛判据 $\epsilon$ ， $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$ ：

$\quad \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\geq 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})|| \leq \epsilon$ ，则停止迭代，得到最优解 $\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^{'}(\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)。

限域拟牛顿法(Limited Storage Quasi-Newton Method)

算法2的步骤(3)中，为了确定第 $k$ 次搜索方向，需要知道对称正定矩阵 $\mathbf{H}_k$ ，因此对于 $N$ 维的问题，存储空间至少是 $N(N+1)/2$ ，对于大型计算而言，这显然是一个极大的缺点。作为比较，共轭梯度法只需要存储3个 $N$ 维向量。为了解决这个问题，Nocedal首次提出了基于BFGS公式的限域拟牛顿法(L-BFGS)[2]。

L-BFGS的基本思想是存储有限次数（如 $m$ 次）的更新矩阵 $\Delta\mathbf{H}_k$ ，如果 $k$ > $m$ 的话就舍弃 $m$ 次以前的 $\Delta\mathbf{H}_{k-m+1}$ ，也即L-BFGS的记忆只有 $m$ 次。如果 $m=N$ ，则L- BFGS等价于标准的BFGS公式。

首先将方程(7)写为乘法形式：

其中 $\rho_k=\frac{1}{\mathbf{y}^{\rm T}_k\mathbf{s}_k}$ ， $\mathbf{v}_k$ 是 $N$ $\times$ $N$ 的矩阵。乘法形式下"舍弃"等价于置 $\mathbf{v}_k=\mathbf{I}$ ， $\rho_k=0$ 。容易得出，给定 $m$ 后，BFGS的矩阵更新可以写为

若 $k+1\leq m$ ，则

$\mathbf{H}_{k+1}=\mathbf{v}^{\rm T}_k\mathbf{v}^{\rm T}_{k-1}\cdots \mathbf{H}_0\mathbf{v}_0\cdots \mathbf{v}_{k-1}\mathbf{v}_k+\mathbf{v}^{\rm T}_k\cdots\mathbf{v}^{\rm T}_1\rho_0\mathbf{s}_0\mathbf{s}^{\rm T}_0\mathbf{v}_1\cdots\mathbf{v}_k$
       $+\mathbf{v}^{\rm T}_k\mathbf{v}^{\rm T}_{k-1}\cdots\mathbf{v}^{\rm T}_2\rho_1\mathbf{s}_1\mathbf{s}^{\rm T}_1\mathbf{v}_2\cdots\mathbf{v}_k$
         $\cdots$
       $+\mathbf{v}^{\rm T}_k\mathbf{v}^{\rm T}_{k-1}\rho_{k-2}\mathbf{s}_{k-2}\mathbf{s}^{\rm T}_{k-2}\mathbf{v}_{k-1}\mathbf{v}_k$
   $+\mathbf{v}^{\rm T}_k\rho_{k-1}\mathbf{s}_{k-1}\mathbf{s}^{\rm T}_{k-1}\mathbf{v}_k$
       $+\rho_k\mathbf{s}_k\mathbf{s}^{\rm T}_k$           (9)

若 $k+1$ > $m$ ，则

$\mathbf{H}_{k+1}=\mathbf{v}^{\rm T}_k\mathbf{v}^{\rm T}_{k-1}\cdots\mathbf{v}^{\rm T}_{k-m+1}\mathbf{H}_0\mathbf{v}_{k-m+1}\cdots\mathbf{v}_{k-1}\mathbf{v}_k$
         $+\mathbf{v}^{\rm T}_k\cdots\mathbf{v}^{\rm T}_{k-m+2}\rho_{k-m+1}\mathbf{s}_{k-m+1}\mathbf{s}^{\rm T}_{k-m+1}\mathbf{v}_{k-m+2}\cdots\mathbf{v}_k$
       $\cdots$
       $+\mathbf{v}^{\rm T}_k\rho_{k-1}\mathbf{s}_{k-1}\mathbf{s}^{\rm T}_{k-1}\mathbf{v}_k$
   $+\rho_k\mathbf{s}_k\mathbf{s}^{\rm T}_k$          (10)

方程(9)和(10)称为狭义BFGS矩阵(special BFGS matrices)。仔细分析这两个方程以及 $\rho_k$ 和 $\mathbf{v}_k$ 的定义，可以发现L-BFGS方法中构造 $\mathbf{H}_{k+1}$ 只需要保留 $2m+1$ 个 $N$ 维向量： $m$ 个 $\mathbf{s}_k$ 、 $m$ 个 $\mathbf{y}_k$ 以及 $\mathbf{H}_0$ （对角阵）。

快速计算 $\mathbf{H}\cdot\mathbf{g}$

L-BFGS算法中确定搜索方向 $\mathbf{d}_k$ 需要计算 $\mathbf{H}\cdot\mathbf{g}$ ，下列算法可以高效地完成计算任务：

算法3:

IF $k \leq m$ Then

${\rm incr}$ = 0; ${\rm BOUND}$ = $k$

ELSE

${\rm incr}$ = $k-m$ ; ${\rm BOUND}$ = $m$

ENDIF

$\mathbf{q}_{\rm BOUND} = \mathbf{g}_k$ ;

DO $i$ = ( ${\rm BOUND}$ -1),0,-1

$\quad j = i+{\rm incr}$ ;

$\quad \alpha_i = \rho_j\mathbf{s}^{\rm T}_j\mathbf{q}_{i+1}$ ;

$\quad$ 储存 $\alpha_i$ ;

$\mathbf{q}_i = \mathbf{q}_{i+1}-\alpha_i\mathbf{y}_j$ ;

ENDDO

$\mathbf{r}_0 = \mathbf{H}_0\cdot\mathbf{g}_0$ ;

DO $i$ = 0, ( ${\rm BOUND}$ -1)

$\quad j = i+{\rm incr}$ ;

$\quad \beta_j = \rho_j\mathbf{y}^{\rm T}_j\mathbf{r}_i$ ;

$\quad \mathbf{r}_{i+1} = \mathbf{r}_i+\mathbf{s}_j(\alpha_i-\beta_i)$ ;

ENDDO

完整的程序包可从下列网址下载：

http://www.ece.northwestern.edu/~nocedal/software.html

针对二次非凸函数的若干变形

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

Li-Fukushima方法[3]

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

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

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

其中

$\mathbf{y}^{\ast}_k=\mathbf{g}_{k+1}-\mathbf{g}_k+t_k||\mathbf{g}_k||\mathbf{s}_k$
$t_k=1+\max\{0, \frac{-\mathbf{y}^{\rm T}_k\mathbf{s}_k}{||\mathbf{s}_k||^2}\}$

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

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

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

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

$f(\mathbf{x}_k+\epsilon_j\mathbf{d}_k)\leq f(\mathbf{x}_k)+\sigma\epsilon_j\mathbf{g}^{\rm T}_k\mathbf{d}_k$

取 $j_k=j$ ，步长 $\lambda_k=\epsilon_{j_k}$ 。

Xiao-Wei-Wang方法[4]

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

设 $\mathbf{y}^{\dagger}_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)^{\rm T}\mathbf{s}_k]$

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

$\mathbf{H}_{k+1}=(\mathbf{I}-\rho^{\dagger}_k\mathbf{s}_k\mathbf{y}^{\dagger\rm T}_k)\mathbf{H}_k(\mathbf{I}-\rho^{\dagger}_k\mathbf{y}^{\dagger}_k)\mathbf{s}^{\rm T}_k+\rho^{\dagger}_k\mathbf{s}_k\mathbf{s}^{\rm T}_k$ (12)
相应的 $\rho^{\dagger}_k=\frac{1}{\mathbf{y}^{\dagger\rm T}_k\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) \leq f(\mathbf{x}_k)+\delta\lambda_k\mathbf{g}^{\rm T}_k\mathbf{d}_k$ (13)

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

如果 $\lambda_k$ = $1$ 满足方程(13)、(14)，则取 $\lambda_k$ = $1$ 。

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

一维搜索

使用导数的优化算法都涉及到沿优化方向 $\mathbf{d}_k$ 的一维搜索。事实上一维搜索算法本身就一个非常重要的课题，分为精确一维搜索以及非精确一维搜索。标准的拟牛顿法或L-BFGS均采用精确一维搜索。与前者相比，非精确一维搜索虽然牺牲了部分精度，但是效率更高，调用函数的次数更少。因此 Li-Fukushima方法和Xiao-Wei-Wang方法中均采用了这类算法。不加证明的，本节分别给出两类范畴中各自的一个应用最为广泛的例子，分别是二点三次插值方法和Wolfe-Powell准则。

二点三次插值方法

在精确一维搜索各种算法中，这种方法得到的评价最高。其基本思想是选取两个初始点 $x_1$ 和 $x_2$ ，满足 $x_1$ < $x_2$ ; $f^{'}(x_1)$ < $0$ ; $f^{'}(x_2)$ > $0$ 。这样的初始条件保证了在区间 $(x_1, x_2)$ 中存在极小点。利用这两点处的函数值 $f(x_1)$ 、 $f(x_2)$ （记为 $f_1$ 、 $f_2$ ）和导数值 $f^{'} (x_1)$ 、 $f^{'}(x_2)$ （记为 $f^{'}_1$ 、 $f^{'}_2$ ）构造一个三次多项式 $\varphi(x)$ ，使得 $\varphi(x)$ 在 $x_1$ 和 $x_2$ 处与目标函数 $f(x)$ 有相同的函数值和导数值，则设 $\varphi(x)=a(x- x_1)^3+b(x-x_1)^2+c(x-x_1)+d$ ，那么通过4个边界条件可以完全确定 $a$ 、 $b$ 、 $c$ 、 $d$ 4个参数。之后找出 $\varphi^{'}(x)$ 的零点 $x^{'}$ ，作为极小点的一个进一步的估计。可以证明，由 $x_1$ 出发，最佳估计值的计算公式为

$x^{'}=x_1+\frac{-c}{b+\sqrt{b^2-3ac}}$ (15)

为了避免每次都要求解4维线性方程组的麻烦，整个搜索过程可以采用算法4：

算法4：

(1) 给定初始点 $x_1$ 和 $x_2$ ，满足 $x_1$ < $x_2$ ，计算函数值 $f_1$ 、 $f_2$ 和导数值 $f^{'}_1$ 、 $f^{'}_2$ ，并且 $f^{'}_1$ < $0$ ， $f^{'}_2$ > $0$ ，给定允许误差 $\delta$ 。

(2) 计算如下方程，得到最佳估计值 $x^{'}$ ：

$s=\frac{3(f_2-f_1)}{x_2-x_1},\quad z = s-f^{'}_1-f^{'}_2\quad w^2 = z^2-f^{'}_1f^{'}_2$ (16)

$x^{'} = x_1+(x_2-x_1)\frac{1-f^{'}_2+w+z}{f^{'}_2-f^{'}_1+2w}$ (17)

(3) 若 $|x_2-x_1|\leq \delta$ ，则停止计算，得到点 $x^{'}$ ；否则进行步骤(4)。

(4) 计算 $f(x^{'})$ 和 $f^{'}(x^{'})$ 。若 $f^{'}(x^{'}) = 0$ ，则停止计算，得到点 $x^{'}$ ,

若 $f^{'}(x^{'})$ < $0$ ，则令 $x_1 = x^{'}; f_1 = f(x^{'}); f^{'}_1 = f^{'}(x^{'})$ ，转步骤(2)；

若 $f^{'}(x^{'})$ > $0$ ，则令 $x_2 = x^{'}; f_2 = f(x^{'}); f^{'}_2 = f^{'}(x^{'})$ ，转步骤(2)。

利用三次函数插值，方程(16)、(17)并不是唯一的方法，也可以利用下式计算 $a$ 、 $b$ 、 $c$ 三个参数：

$a = \frac{f^{'}_1+f^{'}_2}{(x_2-x_1)^2}-\frac{2(f_2-f_1)}{(x_2-x_1)^3};$
$b = \frac{2f^{'}_1-f^{'}_2}{x_2-x_1}+\frac{3}{2}\frac{f_2-f_1}{(x_2-x_1)^2};$ (18)
$c = f^{'}_1.$

然后利用(15)式寻找最佳点 $x^{'}$ [5]。此外，即使 $f^{'}(x_2)$ < $0$ ，一般而言也可以用(15)式外推寻找 $x^{'}$ （参见[5]）。

Wolfe-Powell准则

不等式(13)、(14)给出了这种非精确一维搜索算法。如果我们将不等式(14)用下式替换：

$|\mathbf{g}^{\rm T}_{k+1}\mathbf{d}_k|\leq -\sigma\mathbf{g}^{\rm T}_k\mathbf{d}_k$ (19)

也即

$\sigma\mathbf{g}^{\rm T}_k\mathbf{d}_k \leq \mathbf{g}^{\rm T}_{k+1}\mathbf{d}_k \leq -\sigma\mathbf{g}^{\rm T}_k\mathbf{d}_k$

则不等式(13)、(19)称为强Wolfe-Powell准则。其重要性在于当 $\sigma\rightarrow 0$ 时，该方法过渡为精确一维搜索算法[6]。算法如下[5]

算法5：

(1) 给定两个参数 $\delta\in (0,1/2)$ 和 $\sigma\in (\delta,1)$ ， $x_1$ 为初始点（相应于 $\lambda_k = 0$ ）， $x_2$ 为猜想点（可设为1）。计算两点处的函数值 $f_1$ 、 $f_2$ 和导数值 $f^{'}_1$ 、 $f^{'}_2$ 。给定最大循环次数 $N_{\rm max}$ ，设 $k = 0; i = 0$ ；

(2) 若 $f_2$ 和 $f^{'}_2$ 违反不等式(13)或者不等式(19)的右半段，则缩小搜索范围的上限 $x_{\rm upper} = x_2$ ，否则转步骤(5)；

(3) 若 $f_2$ > $f_1$ ，利用二次插值方法寻找最佳点 $x_{\rm min}$ ：

$x_{\rm min} = x_1-\frac{1}{2}\frac{f^{'}_1(x_2-x_1)^2}{f_2-f_1-f^{'}_1(x_2-x_1)}$

设 $x_2 = x_{\rm min}$ ，计算 $f_2$ 和 $f^{'}_2$ ；设 $k = k+1$ ，若 $k \leq N_{\rm max}$ 转步骤(2)，否则转步骤(5);

若 $f_2\leq f_1$ ，转步骤(4)；

(4) 利用式(16)、(17)（或者式(15)、(18)）寻找最佳点 $x_{\rm min}$ 。令 $x_2 = x_{\rm min}$ ，计算 $f_2$ 和 $f^{'}_2$ ；设 $k$ = $k+1$ ，若 $k\leq N_{\rm max}$ ，转步骤(2)，否则转步骤(5)；

(5) 若 $f^{'}_2$ 满足不等式(19)的左半段，则停止计算，得到最佳点 $x_2$ ；否则转步骤(6)；

(6) 利用式(16)、(17)（或者式(15)、(18)）寻找最佳点 $x_{\rm min}$ ，并计算 $f_2$ 以及 $f^{'}_2$ ；设 $i = i+1$ ，若 $i \leq N_{\rm max}$ 转步骤(2)，否则转步骤(7)；

(7) 停止计算得到目前最佳估计值 $x_2$ 。

需要补充说明的是步骤(4)可以有不同的估算方法，例如

$x_{\rm min} = x_2-\frac{f^{'}_2(x_2-x_1)}{f^{'}_2-f^{'}_1}$ (20)

此外，点 $x$ 处的导数值 $f^{'}(x) = \mathbf{g}^{\rm T}\mathbf{d}$ ，因为在一维搜索中， $x$ 相当于待求步长 $\lambda$ 。在大多数情况下， $\sigma$ = $0.4$ 以及 $\rho$ = $0.1$ 可以取得很好的效果。Wolfe-Powell准则的几何意义可以参考文献[5][6]。

参考文献

【1】陈宝林《最优化理论与算法（第二版）》（清华大学出版社 2005）第9-10章.
【2】 J. Nocedal, Math. Comput. 35, 773 (1980).
【3】 D.H. Li and M. Fukushima, J. Comput. Appl. Math. 129, 15 (2001).
【4】 Y. Xiao, Z. Wei and Z. Wang, Comput. Math. Appl. 56, 1001 (2008).
【5】 www.cs.toronto.edu/~delve/methods/mlp-ese-1/minimize.ps.gz
【6】 W. Sun and Y.X. Yuan, Optimization theory and methods (Springer 2006), p. 104.

博客介绍

正文

拟牛顿法及其相关解法2011-11-15 11:41:00