优化算法——拟牛顿法之DFP算法

来源：互联网发布：苏州园区网络托管编辑：程序博客网时间：2024/04/30 06:30

一、牛顿法

在博文“优化算法——牛顿法(Newton Method)”中介绍了牛顿法的思路，牛顿法具有二阶收敛性，相比较最速下降法，收敛的速度更快。在牛顿法中使用到了函数的二阶导数的信息，对于函数 $f\left ( x \right )$ ，其中 $x$ 表示向量。在牛顿法的求解过程中，首先是将函数 $f\left ( x \right )$ 在 $x=x_{k+1}$ 处展开，展开式为：

$f\left ( x \right )=f\left ( x_{k+1} \right )+g_{k+1}^T\left ( x-x_{k+1} \right )+\frac{1}{2}\left ( x-x_{k+1} \right )^TG_{k+1}\left ( x-x_{k+1} \right )+o\left ( x-x_{k+1} \right )^2$

其中， $g_{k+1}^T=\bigtriangledown f\left ( x_{k+1} \right )$ ，表示的是目标函数在 $x=x_{k+1}$ 的梯度，是一个向量。 $G_{k+1}=\bigtriangledown ^2f\left ( x_{k+1} \right )$ ，表示的是目标函数在 $x=x_{k+1}$ 处的Hesse矩阵。省略掉最后面的高阶无穷小项，即为：

$f\left ( x \right )=f\left ( x_{k+1} \right )+g_{k+1}^T\left ( x-x_{k+1} \right )+\frac{1}{2}\left ( x-x_{k+1} \right )^TG_{k+1}\left ( x-x_{k+1} \right )$

上式两边对 $x$ 求导，即为：

$\bigtriangledown f\left ( x \right )=\bigtriangledown f\left ( x_{k+1} \right )+G_{k+1}\left ( x-x_{k+1} \right )$

在基本牛顿法中，取得最值的点处的导数值为 $0$ ，即上式左侧为 $0$ 。则：

$\bigtriangledown f\left ( x_{k+1} \right )+G_{k+1}\left ( x-x_{k+1} \right )=0$

求出其中的 $x$ ：

$x=x_{k+1}-G_{k+1}^{-1}\bigtriangledown f\left ( x_{k+1} \right )$

从上式中发现，在牛顿法中要求Hesse矩阵是可逆的。

当 $x=x_k$ 时，上式为：

$\bigtriangledown f\left ( x_k \right )=\bigtriangledown f\left ( x_{k+1} \right )+G_{k+1}\left ( x_k-x_{k+1} \right )$

此时，是否可以通过 $x_{k}$ ， $x_{k+1}$ ， $\bigtriangledown f\left ( x_k \right )$ 和 $\bigtriangledown f\left ( x_{k+1} \right )$ 模拟出Hesse矩阵的构造过程？此方法便称为拟牛顿法(QuasiNewton)，上式称为拟牛顿方程。在拟牛顿法中，主要包括DFP拟牛顿法，BFGS拟牛顿法。

二、DFP拟牛顿法

1、DFP拟牛顿法简介

DFP拟牛顿法也称为DFP校正方法，DFP校正方法是第一个拟牛顿法，是有Davidon最早提出，后经Fletcher和Powell解释和改进，在命名时以三个人名字的首字母命名。

对于拟牛顿方程：

$\bigtriangledown f\left ( x_k \right )=\bigtriangledown f\left ( x_{k+1} \right )+G_{k+1}\left ( x_k-x_{k+1} \right )$

化简可得：

$G_{k+1}^{-1}\left [ \bigtriangledown f\left ( x_{k+1} \right )-\bigtriangledown f\left ( x_k \right ) \right ]=x_{k+1}-x_k$

令 $H_{k+1}\overset{\bigtriangleup }{=}G_{k+1}^{-1}$ ，可以得到：

$H_{k+1}\left [ \bigtriangledown f\left ( x_{k+1} \right )-\bigtriangledown f\left ( x_k \right ) \right ]=x_{k+1}-x_k$

在DFP校正方法中，假设：

$H_{k+1}=H_k+E_k$

2、DFP校正方法的推导

令： $E_k=\alpha u_ku_k^T+\beta v_kv_k^T$ ，其中 $u_k,v_k$ 均为 $n\times 1$ 的向量。 $y_k=\bigtriangledown f\left ( x_{k+1} \right )-\bigtriangledown f\left ( x_k \right )$ ， $s_k=x_{k+1}-x_k$ 。

则对于拟牛顿方程 $H_{k+1}\left [ \bigtriangledown f\left ( x_{k+1} \right )-\bigtriangledown f\left ( x_k \right ) \right ]=x_{k+1}-x_k$ 可以简化为：

$H_{k+1}y_k=s_k$

将 $H_{k+1}=H_k+E_k$ 代入上式：

$\left ( H_k+E_k \right )y_k=s_k$

将 $E_k=\alpha u_ku_k^T+\beta v_kv_k^T$ 代入上式：

$\left ( H_k+\alpha u_ku_k^T+\beta v_kv_k^T \right )y_k=s_k$

$\Rightarrow \alpha \left ( u_k^Ty_k \right )u_k+\beta \left ( v_k^Ty_k \right )v_k=s_k-H_ky_k$

已知： $u_k^Ty_k,v_k^Ty_k$ 为实数， $s_k-H_ky_k$ 为 $n\times 1$ 的向量。上式中，参数 $\alpha$ 和 $\beta$ 解的可能性有很多，我们取特殊的情况，假设 $u_k=rH_ky_k$ ， $v_k=\theta s_k$ 。则：

代入上式：

$\Rightarrow \alpha \left [ \left ( rH_ky_k\right )^Ty_k \right ]\left ( rH_ky_k \right )+\beta \left [ \left ( \theta s_k \right )^Ty_k \right ]\left ( \theta s_k \right )=s_k-H_ky_k$

$\Rightarrow \left [ \alpha r^2\left ( y_k^TH_ky_k\right ) +1\right ]\left ( H_ky_k \right )+\left [ \beta \theta ^2\left ( s_k^Ty_k \right ) -1\right ]\left (s_k \right )=0$

令 $\alpha r^2\left ( y_k^TH_ky_k\right ) +1=0$ ， $\beta \theta ^2\left ( s_k^Ty_k \right ) -1=0$ ，则：

$\alpha r^2=-\frac{1}{y_k^TH_ky_k}$

$\beta \theta ^2=\frac{1}{s_k^Ty_k}$

则最终的DFP校正公式为：

$H_{k+1}=H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k}$

3、DFP拟牛顿法的算法流程

设 $H_k$ 对称正定， $H_{k+1}$ 由上述的DFP校正公式确定，那么 $H_{k+1}$ 对称正定的充要条件是 $s_k^Ty_k> 0$ 。

在博文“优化算法——牛顿法(Newton Method)”中介绍了非精确的线搜索准则：Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对DFP校正公式做些许改变：

$H_{k+1}=\begin{cases} H_k & \text{ if } s_k^Ty_k\leq 0 \\ H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ty_k} & \text{ if } s_k^Ty_k>0 \end{cases}$

DFP拟牛顿法的算法流程如下：

4、求解具体的优化问题

求解无约束优化问题

$min\; f\left ( x \right )=100\left ( x_1^2-x_2 \right )^2+\left ( x_1-1 \right )^2$

其中， $x=\left ( x_1,x_2 \right )^T\in\mathbb{R}^2$ 。

python程序实现：

function.py

#coding:UTF-8'''Created on 2015年5月19日@author: zhaozhiyong'''from numpy import *#fundef fun(x):    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2#gfundef gfun(x):    result = zeros((2, 1))    result[0, 0] = 400 * x[0,0] * (x[0,0] ** 2 - x[1,0]) + 2 * (x[0,0] - 1)    result[1, 0] = -200 * (x[0,0] ** 2 - x[1,0])    return result

dfp.py

#coding:UTF-8'''Created on 2015年5月19日@author: zhaozhiyong'''from numpy import *from function import *def dfp(fun, gfun, x0):    result = []    maxk = 500    rho = 0.55    sigma = 0.4    m = shape(x0)[0]    Hk = eye(m)    k = 0    while (k < maxk):        gk = mat(gfun(x0))#计算梯度        dk = -mat(Hk)*gk        m = 0        mk = 0        while (m < 20):            newf = fun(x0 + rho ** m * dk)            oldf = fun(x0)            if (newf < oldf + sigma * (rho ** m) * (gk.T * dk)[0,0]):                mk = m                break            m = m + 1                #DFP校正        x = x0 + rho ** mk * dk        sk = x - x0        yk = gfun(x) - gk        if (sk.T * yk > 0):            Hk = Hk - (Hk * yk * yk.T * Hk) / (yk.T * Hk * yk) + (sk * sk.T) / (sk.T * yk)                k = k + 1        x0 = x        result.append(fun(x0))        return result

testDFP.py

#coding:UTF-8'''Created on 2015年5月19日@author: zhaozhiyong'''from bfgs import *from dfp import dfpimport matplotlib.pyplot as plt  x0 = mat([[-1.2], [1]])result = dfp(fun, gfun, x0)n = len(result)ax = plt.figure().add_subplot(111)x = arange(0, n, 1)y = resultax.plot(x,y)plt.show()

5、实验结果

0 0