拟牛顿法之BFGS算法

什么是拟牛顿法？

拟牛顿法是在牛顿法的基础上引入了Hessian矩阵的近似矩阵，避免每次迭代都计算Hessian矩阵的逆，它的收敛速度介于梯度下降法和牛顿法之间。拟牛顿法跟牛顿法一样，也是不能处理太大规模的数据，因为计算量和存储空间会开销很多。
拟牛顿法虽然每次迭代不像牛顿法那样保证是最优化的方向，但是近似矩阵始终是正定的，因此算法始终是朝着最优化的方向在搜索。具有全局收敛性和超线性收敛速度

BFGS公式推导

BFGS(Broyden，Fletcher，Goldfarb，Shanno四个人)算法是使用较多的一种拟牛顿方法，故称为BFGS校正。

将x写成x=(x1,x2,…,xn)。对函数f(x)在x=xk+1处进行泰勒展开到二阶：

f(x)=f(xk+1)+f′(xk+1)(x−xk+1)+12f′′(xk+1)(x−xk+1)2+Rn(x)≈f(xk+1)+f′(xk+1)(x−xk+1)+12f′′(xk+1)(x−xk+1)2

对上式求导并令其为0，由于f(x)中的x是一个向量，f(x)对x求导意味着对x向量中的每个值求偏导。即，f(x)对x的一阶导数为一个向量，对x的二阶导数为一个n∗n的矩阵

f′(x)=(∂f(x)∂x1，∂f(x)∂x2,…,∂f(x)∂xn)f′′(x)=[∂2f(x)∂xi∂xj]n∗n

求导后得:

f′(x)=f′(xk+1)+f′′(xk+1)(x−xk+1)

即：

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

可以化简为:

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

令Bk+1≜Gk+1,则可得:Bk+1(xk−xk+1)=∇f(xk+1)−∇f(xk)
在BFGS校正方法中，假设：

Bk+1=Bk+Ek

BFGS校正公式的推导

令Ek=αukuTk+βvkvTk，其中uk,vk均为n∗1的向量。yk=∇f(xk+1)−∇f(xk),sk=xk+1−xk.

那么Bk+1(xk−xk+1)=∇f(xk+1)−∇f(xk)
可以化简为：

Bk+1sk=yk

将Bk+1=Bk+Ek代入上式得：

(Bk+Ek)sk=yk

将Ek=αukuTk+βvkvTk代入上式得：

(Bk+αukuTk+βvkvTk)sk=yk

即：

αuk(uTksk)+βvk(vTksk)=yk−Bksk

uTksk,vTksk皆为实数，yk−Bksk为n∗1的向量，上式中，参数α和β解的可能性有很多，我们取特殊的情况，假设uk=rBksk,vk=θyk。则:

Ek=αrBksTkBk+βθykyTk

代入上式：

⇒α[(rBksk)Tsk](rBksk)+β[(θyk)Tsk](θyk)=yk−Bksk⇒[αr2(sTkBksk)+1](Bksk)+[βθ2(yTksk)−1](yk)=0

令⇒[αr2(sTkBksk)+1](Bksk)=0,βθ2(yTksk)−1=0,则：

αr2=−1sTkBkskβθ2=1yTksk

最终的BFGS校正公式为：

Bk+1=Bk−BksksTkBksTkBksk+ykyTkyTksk

BFGS校正的算法流程

设Bk对称正定，Bk+1由上述的BFGS校正公式确定，那么Bk+1对称正定的充要条件是yTksk>0。

非精确的一维搜索（线搜索）准则：Armijo搜索准则，搜索准则的目的是为了帮助我们确定学习率，还有其他的一些准则，如Wolfe准则以及精确线搜索等。在利用Armijo搜索准则时并不是都满足上述的充要条件，此时可以对BFGS校正公式做些许改变：

Bk+1=⎧⎩⎨Bk,Bk−BksksTkBksTkBksk+ykyTkyTksk,ifyTksk≤0ifyTksk>0

注：在李航写的那本《统计学习方法》中说是正定的，但是并没有说上述情况下会怎么样

算法

1. 给定参数δ∈(0,1),σ∈(0,0.5)，初始化点x0∈Rn，终止误差0≤ϵ≪1,初始化对称正定阵B0，通常取为G(xo)或单位阵In;令k=0。
2. 计算gk=∇f(xk)，若∥gk∥≪ϵ,终止，输出xk作为近似极小点。
3. 解线性方程组得解dk:Bkd=−gk.
4. 设mk是满足下列不等式的最小非负整数m:
   
   f(xk+δmdk)≤f(xk)+σδmgTkdk
   
   令αk=δmk,xk+1=xk+αkdk.
5. 由BFGS校正公式确定Bk+1
6. 令k=k+1，转向步骤“2”

求解具体优化问题

求解无约束优化问题:

minf(s)=100(x21−x2)2+(x1−1)2,x=(x1,x2)T∈R2

#coding:UTF-8  '''Created on 2017年4月20日@author: zhangdapeng'''from numpy import *    import matplotlib.pyplot as pltfrom numpy.matrixlib.defmatrix import mat#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  def bfgs(fun, gfun, x0):      result = []      maxk = 500      delta = 0.55      sigma = 0.4      m = shape(x0)[0]      Bk = eye(m)      k = 0    epsilon=1e-10    while (k < maxk):          gk = mat(gfun(x0))#计算梯度 ，mat函数将数组转化为矩阵。#         print(gk) #         print(linalg.norm(gk,1))        #axis=0,沿着纵轴方向        if linalg.norm(gk,1)<epsilon:            break        dk = mat(-linalg.solve(Bk, gk))  #解矩阵方程Bk*x=gk得到x        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          #BFGS校正          x = x0 + delta ** mk * dk          sk = x - x0          yk = gfun(x) - gk  #         print(math.isnan(yk.T * sk))        if (yk.T * sk > 0):             Bk = Bk - (Bk * sk * sk.T * Bk) / (sk.T * Bk * sk) + (yk * yk.T) / (yk.T * sk)          k = k + 1          x0 = x          result.append(fun(x0))      return result  #初始化x0  x0 = mat([[-1.2], [1]])  result = bfgs(fun, gfun, x0)  print("result:",result[-1])n = len(result)  ax = plt.figure().add_subplot(111)  x = arange(0, n, 1)  y = result  ax.plot(x,y)  plt.show()  

输出：

result: 2.68262011582e-28

输出图片：

http://blog.csdn.net/google19890102/article/details/45867789

http://blog.csdn.net/acdreamers/article/details/44664941
http://www.codelast.com/%E5%8E%9F%E5%88%9B%E7%94%A8%E4%BA%BA%E8%AF%9D%E8%A7%A3%E9%87%8A%E4%B8%8D%E7%B2%BE%E7%A1%AE%E7%BA%BF%E6%90%9C%E7%B4%A2%E4%B8%AD%E7%9A%84armijo-goldstein%E5%87%86%E5%88%99%E5%8F%8Awo/