目录
- 拟牛顿法
1.1拟牛顿法的导出与优点
1.2 算法步骤与特点 - 对称秩一校正公式
- DFP算法
3.1 DFP公式推导
3.2 要求解的问题
3.3 python实现
1.拟牛顿法
1.1拟牛顿法的导出与优点
在上一文中（牛顿法公式推导与python实现），谈到说牛顿法需要计算一个Hessian矩阵的逆，才能够迭代，但在实际工程中，计算如此大型的矩阵需要很大的计算资源，因此，有人提出能否不计算Hessian矩阵，在迭代过程中，仅仅利用相邻两个迭代点以及梯度信息，产生一个对称正定矩阵，使之逐步逼近目标函数Hessian矩阵的逆阵。
其实这就是你牛顿法的基本思想，这样做，既能保存Hessian矩阵的大部分信息（曲率），也能极大的减小计算量。
考虑无约束极小化问题。假设目标函数f:Rn→Rf:Rn→R是二次连续可微的，那么∇f∇f在xk+1xk+1处的泰勒展开为：
∇f(x)=∇f(xk+1)+∇2f(xk+1)(x−xk+1)+o||x−xk+1||∇f(x)=∇f(xk+1)+∇2f(xk+1)(x−xk+1)+o||x−xk+1||
,取x:=xkx:=xk.当xk与xk+1xk与xk+1充分接近时，有：
∇2f(xk+1)(xk+1−xk)≈∇f(xk+1)−∇f(xk)∇2f(xk+1)(xk+1−xk)≈∇f(xk+1)−∇f(xk)
∇2f(xk+1)∇2f(xk+1)就是f(x)f(x)在xk+1xk+1处的Hessian矩阵，那么我们可以用它的近似矩阵Bk+1Bk+1来代替它，得到如下等式：
Bk+1(xk+1−xk)=∇f(xk+1)−∇f(xk)(1)(1)Bk+1(xk+1−xk)=∇f(xk+1)−∇f(xk)
如该矩阵存在逆矩阵有：
Hk+1(∇f(xk+1)−∇f(xk))=xk+1−xk(2)(2)Hk+1(∇f(xk+1)−∇f(xk))=xk+1−xk
以上两个方程成为拟牛顿方程（条件）。其中Hk+1=∇2f(xk+1)−1Hk+1=∇2f(xk+1)−1,为Hessian的逆阵。
1.2 算法步骤与特点
拟牛顿法的算法步骤如下：
- 给出x0∈Rn,H0∈Rnxn,0≤ϵ<1,k:=0x0∈Rn,H0∈Rnxn,0≤ϵ<1,k:=0;
- 若|∇f(xk)|≤ϵ|∇f(xk)|≤ϵ,迭代停止；否则求方向：dk=−Hk∇f(xk)dk=−Hk∇f(xk)
- 沿着方向做线性搜索αk>0αk>0,令xk+1=xk+αkdkxk+1=xk+αkdk
- 校正Hk产生Hk+1,使得牛顿条件（2）Hk产生Hk+1,使得牛顿条件（2）依然成立
- k:=k+1,转至第二步
总结一下拟牛顿法的特点：
- 这种情况下，Hk+1≈∇2f(xk+1)−1Hk+1≈∇2f(xk+1)−1,使得算法产生的方向近似于牛顿方向，从而确保算法具有比较好的收敛性。
- 对任意的k，近似矩阵Bk+1Bk+1都是正定的，使得算法选取的方向(dk=−Hk∇f(xk)dk=−Hk∇f(xk))都是下降方向。
- 仅需一阶导数，就能完整整个迭代过程
- 需要校正HkHk产生Hk+1Hk+1
2.对称秩一校正公式
前面我们说过要用Hk+1Hk+1来近似Hessian的逆阵，但不可能说一次取值，就能得到最优的Hk+1Hk+1,所以我们接下来讨论一下，如何通过迭代，不断的校正这个近似矩阵，使得：
Hk+1=Hk+Ek(3)(3)Hk+1=Hk+Ek
在秩一校正情形下，有：
Hk+1=Hk+uvT(4)(4)Hk+1=Hk+uvT
其中rank(uvT)=1(秩为1)。rank(uvT)=1(秩为1)。
它的想法是希望通过以上这个迭代公式，将u,vTu,vT换成我们可以求得的xk,∇f(x)xk,∇f(x)等，达到的迭代的效果。
令sk=xk+1−xk,yk=∇f(xk+1)−∇f(xk)sk=xk+1−xk,yk=∇f(xk+1)−∇f(xk)将Hk+1Hk+1代入（2）有：
Hk+1yk=(Hk+uvT)yk=sk(5)(5)Hk+1yk=(Hk+uvT)yk=sk
即
(vTyk)u=sk−Hkyk(6)(6)(vTyk)u=sk−Hkyk
故u必定在sk−Hkyksk−Hkyk方向上，且sk−Hkyk≠0sk−Hkyk≠0(如果等于0，则已经满足拟牛顿条件了)，则u=sk−HkykvTyku=sk−HkykvTyk，代入（4），我们有：
Hk+1=Hk+(sk−Hkyk)vTvTyk(7)(7)Hk+1=Hk+(sk−Hkyk)vTvTyk
由于要求Hess矩阵对称，故其逆也必定对称，故v=sk−Hkykv=sk−Hkyk，有
Hk+1=Hk+(sk−Hkyk)(sk−Hkyk)TvTyk(8)(8)Hk+1=Hk+(sk−Hkyk)(sk−Hkyk)TvTyk
该公式被称为对称秩一校正公式，可以用它来校正我们要校正的HkHk.
3.DFP算法
3.1 DFP公式推导
由前面的对称秩一校正公式的导出，我们发现把末尾的未知参数用已知参数代替后，就能完成校正的功能，但对称秩一校正的效果并不是太好，我们可以再加一个校正，让他们协调一下，就有了DFP算法。
DFP算法是设出一个对称秩二校正：
Hk+1=Hk+auuT+bvvT(9)(9)Hk+1=Hk+auuT+bvvT
在满足你牛顿条件的情况下，将式中所有的未知参数a,u,b,va,u,b,v都用已知条件代替，得到一个迭代公式，校正HKHK。
用同样的思想，我们有：
Hkyk+auuTyk+bvvTyk=sk(10)(10)Hkyk+auuTyk+bvvTyk=sk
这里的u,v都不是唯一确定的，但很明显，如果要让等式成立，有：
u=sk,v=Hkyk(11)(11)u=sk,v=Hkyk
与（10）联立，可得：
auTyk=1,bvTyk=−1auTyk=1,bvTyk=−1
确定出：
a=1uTyk=1sTkyk,b=−1yTkHkyka=1uTyk=1skTyk,b=−1ykTHkyk
最后得到DFP公式：
Hk+1=Hk+sksTksTkyk−HkykyTkHkyTkHkykHk+1=Hk+skskTskTyk−HkykykTHkykTHkyk
注意式中的分数结构，分子sTkyk,yTkHkykskTyk,ykTHkyk都是标量，分母sksTk,HkykyTkHkskskT,HkykykTHk则是与HkHk同型的矩阵，且都是正定矩阵。若Hk为Hk为正定矩阵，且sTkyk>0skTyk>0,则Hk+1Hk+1也正定。
当采用精确线搜索时，矩阵序列HkHk的正定新条件sTkyk>0skTyk>0可以被满足。但对于Armijo搜索准则来说，不能满足这一条件，需要做如下修正：
Hk+1=⎧⎩⎨HkHk−HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬Hk+1={HkskTyk≤0Hk−HkykykTHkykTHkyk+skskTskTykskTyk>0}
3.2 要求解的问题
求解无约束线性优化问题
minx∈R2f(x)=100(x21−x2)2+(x1−1)2minx∈R2f(x)=100(x12−x2)2+(x1−1)2
该问题有精确解x∗=(1,1)T,f(x∗)=0x∗=(1,1)T,f(x∗)=0其梯度为
g(x)=(400x1(x21−x2)+2(x1−1),−200(x21−x2))g(x)=(400x1(x12−x2)+2(x1−1),−200(x12−x2))
其Hessian矩阵为
G(x)=[1200x21−400x2+2−400x1−400x1200]G(x)=[1200x12−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),k
x0 ,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）
