拟牛顿法之DFP算法

DFP算法(Davidon-Fletcher-Powell algorithm)一种秩2拟牛顿法.由戴维登(Davidon, W. D.)于1959年导出，并由弗莱彻(Fletcher,R.)和鲍威尔(Powe11,M. J. D.)于1963年进行了改善.

对函数$f(x)$在$x=x_{k+1}$处进行泰勒展开到二阶：

$$\begin{align} f(x)&= f(x_{k+1})+f'(x_{k+1})(x-x_{k+1})+\frac{1}{2}f''(x_{k+1})(x-x_{k+1})^2+R_n(x)\\ & \approx f(x_{k+1})+f'(x_{k+1})(x-x_{k+1})+\frac{1}{2}f''(x_{k+1})(x-x_{k+1})^2 \end{align}$$
对上式求导并令其为0，由于 $f(x)$中的 $x$是一个向量，$f(x)$对 $x$求导意味着对$x$向量中的每个值求偏导。即， $f(x)$对 $x$的一阶导数为一个向量，对$x$的二阶导数为一个 $n*n$的矩阵
$$f'(x)=\left( \frac{\partial f(x)}{\partial x_1}， \frac{\partial f(x)}{\partial x_2},\dots, \frac{\partial f(x)}{\partial x_n}\right)\\ f''(x)=\left[\frac{\partial^2f(x)}{\partial x_i\partial x_j}\right]_{n*n}\\ $$
令 $g_{k+1}^T=\nabla f(x_{k+1})= f'(x_{k+1})$,表示的是目标函数在 $x_{k+1}$的梯度，是一个向量。 $G_{k+1}=\nabla^2f(x_{k+1})=f''(x_{k+1})(x-x_{k+1})^2$表示的是目标函数在 $x_{k+1}$处的Hesse矩阵。
求导后得:
$$\nabla f(x)=\nabla f(x_{k+1})+G_{k+1}(x-x_{k+1})$$
在基本牛顿法中，取得最值的点处的导数值为0，即上式左侧为0。则：
$$\nabla f(x_{k+1})+G_{k+1}(x-x_{k+1})=0$$
求出其中的 $x$： $$x=x_{k+1}-G_{k+1}^{-1}\nabla f(x_{k+1})$$

从上式中发现，在牛顿法中要求Hesse矩阵是可逆的。
当$x=x_k$时，上式为(拟牛顿方程)：

$$\nabla f(k)=\nabla f(x_{k+1})+G_{k+1}(x_k-x_{k+1})$$

此时，是否可以通过$x_k,x_{k+1},\nabla f(x_k),\nabla f(x_{k+1})$模拟出Hesse矩阵的构造过程？此方法便称为拟牛顿法(QuasiNewton)，上式称为拟牛顿方程。

DFP拟牛顿法

DFP拟牛顿法简介

对于拟牛顿方程：

$$\nabla f(k)=\nabla f(x_{k+1})+G_{k+1}(x_k-x_{k+1})$$
化简可得：
$$G_{k+1}^{-1}[\nabla f(x_{k+1})-\nabla f(x_k)]=x_{k+1}-x_k$$
令 $H_{k+1}\triangleq G_{k+1}^{-1}$，可以得到：
$$H_{k+1}[\nabla f(x_{k+1})-\nabla f(x_k)]=x_{k+1}-x_k$$
在DFP校正方法中，假设： $$H_{k+1}=H_k+E_k$$

DFP校正方法的推导

令：$E_k=\alpha u_ku_k^T+\beta v_kv_k^T$，其中$u_k,v_k$均为$n \times 1$的向量。其中：

$$\begin{align} y_k &=\nabla f(x_{k+1})-\nabla f(x_k),\\ s_k&=x_{k+1}-x_k \end{align}$$

则对于$H_{k+1}[\nabla f(x_{k+1})-\nabla f(x_k)]=x_{k+1}-x_k$可以简化为：

$$H_{k+1}y_k=s_k$$
将 $H_{k+1}=H_k+E_k$代入上式：
$$(H_k+E_k)y_k=s_k$$
将 $E_k=\alpha u_ku_k^T+\beta v_kv_k^T$代入上式：
$$\begin{align} (H_k+\alpha u_ku_k^T+\beta v_kv_k^T)y_k&=s_k\\ \Rightarrow \alpha (u_k^Ty_k)u_k+\beta (v_k^Ty_k)v_k&=s_k-H_ky_k \end{align}$$
已知： $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$。则：
$$E_k=\alpha r^2H_ky_ky_k^TH_k+\beta\theta^2s_ks_k^T$$
代入上式：
$$\begin{align} \Rightarrow &\alpha [(rH_ky_k)_k^Ty_k)](rH_ky_k)+\beta [(\theta s_k)_k^Ty_k)(\theta s_k)=s_k-H_ky_k\\ \Rightarrow &[\alpha r^2(y_k^TH_ky_k)+1](H_ky_k)+[\beta \theta^2( s_k^Ty_k)-1](x_k)=0 \end{align}$$
令 $\alpha r^2(y_k^TH_ky_k)+1=0,\beta\theta^2(s_k^Ty_k)-1=0$，则：
$$\begin{align} \alpha r^2 & =-\frac{1}{y_k^TH_ky_k}\\ \beta\theta^2 & =\frac{1}{s_k^Ty_k} \end{align}$$
则最终的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}$$

DFP拟牛顿法的算法流程

设$H_k$对称正定，$H_{k+1}$由上述的DFP校正公式确定，那么$H_{k+1}$对称正定的充要条件是$s_k^Ty_k\gt0$。
Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对DFP校正公式做些许改变：

$$H_{k+1}= \begin{cases} H_k, & if \quad s_k^Ty_k\le0 \\ H_k-\frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\frac{s_ks_k^T}{s_k^Ts_k} ,&if \quad s_k^Ty_k\gt0\end{cases}$$
DFP拟牛顿法的算法流程如下：
1. 给定参数 $\delta \in(0,1),\sigma \in(0,0.5)$，初始化点 $x_0 \in R^n$，终止误差 $0 \le\epsilon \ll1$,初始化对称正定阵 $H_0$，通常取为 $G(x_o)^{-1}$或单位阵 $I_n$;令 $k=0$。
2. 计算 $g_k=\nabla f(x_k)$，若 $\left \| g_k \right \| \ll \epsilon$,终止，输出 $x_k$作为近似极小点。
3. 计算搜索方向 : $d_k=-H_kg_k$.
4. 设 $m_k$是满足下列不等式的最小非负整数m:
$$f(x_k+\delta^m d_k)\le f(x_k)+\sigma \delta^m g_k^Td_k$$
令 $\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k$.
5. 由校正公式确定 $H_{k+1}$
6. 令 $k=k+1$，转向步骤“2”

代码：

dfp.py

#coding:UTF-8
'''
Created on 2017年4月25日

@author: zhangdapeng
'''
from numpy import *
from function import *

def dfp(fun, gfun, x0):
    result = []
    maxk = 500
    delta = 0.55
    sigma = 0.4
    m = shape(x0)[0]
    Hk = eye(m)
    k = 0
    epsilon=1e-10

    while (k < maxk):
        gk = mat(gfun(x0))#计算梯度
        if linalg.norm(gk,1)<epsilon:
            break
        dk = -mat(Hk)*gk
        m = 0
        mk = 0
        while (m < 20):
            newf = fun(x0 + delta ** m * dk)
            oldf = fun(x0)
            if (newf < oldf + sigma * (delta ** m) * (gk.T * dk)[0,0]):
                mk = m
                break
            m = m + 1

        #DFP校正
        x = x0 + delta ** 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

function.py

#coding:UTF-8
'''
Created on 2017年4月24日

@author: zhangdapeng
'''

from numpy import *

#fun  原始函数
def fun(x):  
    return 100 * (x[0,0] ** 2 - x[1,0]) ** 2 + (x[0,0] - 1) ** 2  

#对x1，x2求导后的函数  
def gfun(x):  
    result = zeros((2, 1))  
#     对x1求导
    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])  #对x2求导
    return result

testDFP.py

#coding:UTF-8
'''
Created on 2017年4月25日

@author: zhangdapeng
'''
from dfp import *

import matplotlib.pyplot as plt  

x0 = mat([[-1.2], [1]])
result = dfp(fun, gfun, x0)
print(result[-1])
n = len(result)
ax = plt.figure().add_subplot(111)
x = arange(0, n, 1)
y = result
ax.plot(x,y)

plt.show()

输出结果：

1.54074395551e-28

输出图片：

http://blog.csdn.net/google19890102/article/details/45848439
http://blog.csdn.net/u012102306/article/details/50534086
http://www.codelast.com/%E5%8E%9F%E5%88%9B%E6%8B%9F%E7%89%9B%E9%A1%BF%E6%B3%95quasi-newton%EF%BC%8Cdfp%E7%AE%97%E6%B3%95davidon-fletcher-powell%EF%BC%8C%E5%8F%8Abfgs%E7%AE%97%E6%B3%95broyden-fletcher-goldfarb-shanno/