拟牛顿法公式推导以及python代码实现（一）

目录

1. 拟牛顿法
   1.1拟牛顿法的导出与优点
   1.2 算法步骤与特点
2. 对称秩一校正公式
3. DFP算法
   3.1 DFP公式推导
   3.2 要求解的问题
   3.3 python实现

1.拟牛顿法

1.1拟牛顿法的导出与优点

在上一文中（牛顿法公式推导与python实现），谈到说牛顿法需要计算一个Hessian矩阵的逆，才能够迭代，但在实际工程中，计算如此大型的矩阵需要很大的计算资源，因此，有人提出能否不计算Hessian矩阵，在迭代过程中，仅仅利用相邻两个迭代点以及梯度信息，产生一个对称正定矩阵，使之逐步逼近目标函数Hessian矩阵的逆阵。

其实这就是你牛顿法的基本思想，这样做，既能保存Hessian矩阵的大部分信息（曲率），也能极大的减小计算量。

考虑无约束极小化问题。假设目标函数f:Rn→R是二次连续可微的，那么∇f在xk+1处的泰勒展开为：

∇f(x)=∇f(xk+1)+∇2f(xk+1)(x−xk+1)+o||x−xk+1||

,取x:=xk.当xk与xk+1充分接近时，有：

∇2f(xk+1)(xk+1−xk)≈∇f(xk+1)−∇f(xk)

∇2f(xk+1)就是f(x)在xk+1处的Hessian矩阵，那么我们可以用它的近似矩阵Bk+1来代替它，得到如下等式：

Bk+1(xk+1−xk)=∇f(xk+1)−∇f(xk)(1)

如该矩阵存在逆矩阵有：

Hk+1(∇f(xk+1)−∇f(xk))=xk+1−xk(2)

以上两个方程成为拟牛顿方程（条件）。其中Hk+1=∇2f(xk+1)−1,为Hessian的逆阵。

1.2 算法步骤与特点

拟牛顿法的算法步骤如下：

1. 给出x0∈Rn,H0∈Rnxn,0≤ϵ<1,k:=0;

2. 若|∇f(xk)|≤ϵ,迭代停止；否则求方向：dk=−Hk∇f(xk)

3. 沿着方向做线性搜索αk>0,令xk+1=xk+αkdk

4. 校正Hk产生Hk+1,使得牛顿条件（2）依然成立

5. k:=k+1,转至第二步

总结一下拟牛顿法的特点：

• 这种情况下，Hk+1≈∇2f(xk+1)−1,使得算法产生的方向近似于牛顿方向，从而确保算法具有比较好的收敛性。
• 对任意的k，近似矩阵Bk+1都是正定的，使得算法选取的方向(dk=−Hk∇f(xk))都是下降方向。
• 仅需一阶导数，就能完整整个迭代过程
• 需要校正Hk产生Hk+1

2.对称秩一校正公式

前面我们说过要用Hk+1来近似Hessian的逆阵，但不可能说一次取值，就能得到最优的Hk+1,所以我们接下来讨论一下，如何通过迭代，不断的校正这个近似矩阵，使得：

Hk+1=Hk+Ek(3)

在秩一校正情形下，有：

Hk+1=Hk+uvT(4)

其中rank(uvT)=1(秩为1)。

它的想法是希望通过以上这个迭代公式，将u,vT换成我们可以求得的xk,∇f(x)等，达到的迭代的效果。

令sk=xk+1−xk,yk=∇f(xk+1)−∇f(xk)将Hk+1代入（2）有：

Hk+1yk=(Hk+uvT)yk=sk(5)

即

(vTyk)u=sk−Hkyk(6)

故u必定在sk−Hkyk方向上，且sk−Hkyk≠0(如果等于0，则已经满足拟牛顿条件了)，则u=sk−HkykvTyk，代入（4），我们有：

Hk+1=Hk+(sk−Hkyk)vTvTyk(7)

由于要求Hess矩阵对称，故其逆也必定对称，故v=sk−Hkyk，有

Hk+1=Hk+(sk−Hkyk)(sk−Hkyk)TvTyk(8)

该公式被称为对称秩一校正公式，可以用它来校正我们要校正的Hk.

3.DFP算法

3.1 DFP公式推导

由前面的对称秩一校正公式的导出，我们发现把末尾的未知参数用已知参数代替后，就能完成校正的功能，但对称秩一校正的效果并不是太好，我们可以再加一个校正，让他们协调一下，就有了DFP算法。

DFP算法是设出一个对称秩二校正：

Hk+1=Hk+auuT+bvvT(9)

在满足你牛顿条件的情况下，将式中所有的未知参数a,u,b,v都用已知条件代替，得到一个迭代公式，校正HK。

用同样的思想，我们有：

Hkyk+auuTyk+bvvTyk=sk(10)

这里的u,v都不是唯一确定的，但很明显，如果要让等式成立，有：

u=sk,v=Hkyk(11)

与（10）联立，可得：

auTyk=1,bvTyk=−1

确定出：

a=1uTyk=1sTkyk,b=−1yTkHkyk

最后得到DFP公式：

Hk+1=Hk+sksTksTkyk−HkykyTkHkyTkHkyk

注意式中的分数结构，分子sTkyk,yTkHkyk都是标量，分母sksTk,HkykyTkHk则是与Hk同型的矩阵，且都是正定矩阵。若Hk为正定矩阵，且sTkyk>0,则Hk+1也正定。

当采用精确线搜索时，矩阵序列Hk的正定新条件sTkyk>0可以被满足。但对于Armijo搜索准则来说，不能满足这一条件，需要做如下修正：

Hk+1=⎧⎩⎨⎪⎪HkHk−HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬⎪⎪

3.2 要求解的问题

求解无约束线性优化问题

minx∈ℝ2f(x)=100(x21−x2)2+(x1−1)2

该问题有精确解x∗=(1,1)T,f(x∗)=0其梯度为

g(x)=(400x1(x21−x2)+2(x1−1),−200(x21−x2))

其Hessian矩阵为

G(x)=[1200x21−400x2+2−400x1−400x1200]

3.3 python实现

由1.2的算法步骤，可得：

import numpy as np#函数表达式fun = lambda x:100*(x[0]**2 - x[1]**2)**2 +(x[0] - 1)**2#梯度向量gfun = lambda x:np.array([400*x[0]*(x[0]**2 - x[1]) + 2*(x[0] - 1),-200*(x[0]**2 - x[1])])#Hessian矩阵hess = lambda x:np.array([[1200*x[0]**2 - 400*x[1] + 2,-400*x[0]],[-400*x[0],200]])def dfp(fun,gfun,hess,x0):    #功能：用DFP算法求解无约束问题：min fun(x)    #输入：x0式初始点，fun,gfun，hess分别是目标函数和梯度,Hessian矩阵格式    #输出：x,val分别是近似最优点，最优解，k是迭代次数    maxk = 1e5    rho = 0.05    sigma = 0.4    epsilon = 1e-5 #迭代停止条件    k = 0    n = np.shape(x0)[0]    #将Hessian矩阵初始化为单位矩阵    Hk = np.linalg.inv(hess(x0))    while k < maxk:        gk = gfun(x0)        if np.linalg.norm(gk) < epsilon:            break        dk = -1.0*np.dot(Hk,gk)#         print dk        m = 0;        mk = 0        while m < 20:#用Armijo搜索步长            if fun(x0 + rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):                mk = m                break            m += 1        #print mk        #DFP校正        x = x0 + rho**mk*dk        print "第"+str(k)+"次的迭代结果为："+str(x)        sk = x - x0        yk = gfun(x) - gk        if np.dot(sk,yk) > 0:            Hy = np.dot(Hk,yk)            sy = np.dot(sk,yk) #向量的点积            yHy = np.dot(np.dot(yk,Hk),yk) #yHy是标量            Hk = Hk - 1.0*Hy.reshape((n,1))*Hy/yHy + 1.0*sk.reshape((n,1))*sk/sy        k += 1        x0 = x    return x0,fun(x0),kx0 ,fun0 ,k = dfp(fun,gfun,hess,np.array([0,0]))print x0,fun0,k

输出：

第0次的迭代结果为：[ 0.05  0.  ]第1次的迭代结果为：[ 0.08583333  0.0015    ]第2次的迭代结果为：[ 0.10536555  0.00351201]-----第53次的迭代结果为：[ 1.00007963  1.00015789]第54次的迭代结果为：[ 1.00000251  1.00000578]第55次的迭代结果为：[ 1.00000079  1.00000187]第56次的迭代结果为：[ 1.  1.][ 1.  1.] 7.69713624862e-16 57

迭代57次后得到解（1，1）

reference
梯度-牛顿-拟牛顿优化算法和实现
牛顿法与拟牛顿法学习笔记（四）BFGS 算法
最优化理论与方法-袁亚湘
我的课本-最优化选讲