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

一、牛顿法在博文“优化算法——牛顿法(NewtonMethod)”中介绍了牛顿法的思路，牛顿法具有二阶收敛性，相比较最速下降法，收敛的速度更快。在牛顿法中使用到了函数的二阶导数的信息，对于函数，其中表示向量。在牛顿法的求解过程中，首先是将函数在处展开，展开式为：其中，，表示的是目标函数在的梯度，是一个向量。，表示的是目标函数在处的Hesse矩阵。省略掉

zhiyong_will

36205人浏览 · 2015-05-19 22:11:04

zhiyong_will · 2015-05-19 22:11:04 发布

一、牛顿法

在博文“ 优化算法——牛顿法(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$ 。

在博文“ 优化算法——牛顿法(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$

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 *

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

#gfun
def 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 dfp

import 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 = result
ax.plot(x,y)

plt.show()