什么是拟牛顿法？

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

BFGS公式推导

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

将xx写成x=(x1,x2,…,xn)x=(x_1,x_2,\dots,x_n)。对函数f(x)f(x)在x=xk+1x=x_{k+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

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
对上式求导并令其为0，由于 f(x)f(x)中的 xx是一个向量，f(x)f(x)对 xx求导意味着对xx向量中的每个值求偏导。即， f(x)f(x)对 xx的一阶导数为一个向量，对xx的二阶导数为一个 n∗nn*n的矩阵

f′(x)=(∂f(x)∂x1，∂f(x)∂x2,…,∂f(x)∂xn)f′′(x)=[∂2f(x)∂xi∂xj]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}\\
求导后得:

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

f'(x)=f'(x_{k+1})+f''(x_{k+1})(x-x_{k+1})
即：

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

\nabla f(x_k)=\nabla f(x_{k+1})+G_{k+1}(x_k-x_{k+1})
可以化简为:

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

\nabla f(x_{k+1})-\nabla f(x_k)=G_{k+1}(x_k-x_{k+1})
令 Bk+1≜Gk+1B_{k+1}\triangleq G_{k+1},则可得: Bk+1(xk−xk+1)=∇f(xk+1)−∇f(xk)B_{k+1}(x_k-x_{k+1})=\nabla f(x_{k+1})-\nabla f(x_k)
在BFGS校正方法中，假设：

Bk+1=Bk+Ek

B_{k+1}=B_k+E_k

BFGS校正公式的推导

令Ek=αukuTk+βvkvTkE_k=\alpha u_k u_k^T+\beta v_k v_k^T，其中uk,vku_k,v_k均为n∗1n *1的向量。yk=∇f(xk+1)−∇f(xk),sk=xk+1−xky_k=\nabla f(x_{k+1})-\nabla f(x_k),s_k=x_{k+1}-x_k.

那么Bk+1(xk−xk+1)=∇f(xk+1)−∇f(xk)B_{k+1}(x_k-x_{k+1})=\nabla f(x_{k+1})-\nabla f(x_k)
可以化简为：

Bk+1sk=yk

B_{k+1}s_k=y_k
将 Bk+1=Bk+EkB_{k+1}=B_k+E_k代入上式得：

(Bk+Ek)sk=yk

(B_k+E_k)s_k=y_k
将 Ek=αukuTk+βvkvTkE_k=\alpha u_k u_k^T+\beta v_k v_k^T代入上式得：

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

(B_k+\alpha u_k u_k^T+\beta v_k v_k^T)s_k=y_k

即：

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

\alpha u_k (u_k^Ts_k)+\beta v_k (v_k^Ts_k)=y_k-B_k s_k

uTksk,vTksku_k^Ts_k,v_k^Ts_k皆为实数，yk−Bksky_k-B_k s_k为n∗1n*1的向量，上式中，参数α\alpha和β\beta解的可能性有很多，我们取特殊的情况，假设uk=rBksk,vk=θyku_k=rB_ks_k,v_k=\theta y_k。则:

Ek=αrBksTkBk+βθykyTk

E_k=\alpha r B_ks_k^TB_k+\beta\theta y_ky_k^T
代入上式：

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

\Rightarrow \alpha[(rB_ks_k)^Ts_k](rB_ks_k)+\beta[(\theta y_k)^Ts_k](\theta y_k)=y_k-B_ks_k\\\Rightarrow [\alpha r^2(s_k^TB_ks_k)+1](B_ks_k)+[\beta \theta^2( y_k^Ts_k)-1](y_k)=0
令 ⇒[αr2(sTkBksk)+1](Bksk)=0,βθ2(yTksk)−1=0\Rightarrow [\alpha r^2(s_k^TB_ks_k)+1](B_ks_k)=0,\beta \theta^2( y_k^Ts_k)-1=0,则：

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

\alpha r^2=-\frac{1}{s_k^TB_ks_k}\\ \beta\theta^2=\frac{1}{y_k^Ts_k}
最终的BFGS校正公式为：

Bk+1=Bk−BksksTkBksTkBksk+ykyTkyTksk

B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}

BFGS校正的算法流程

设BkB_k对称正定，Bk+1B_{k+1}由上述的BFGS校正公式确定，那么Bk+1B_{k+1}对称正定的充要条件是yTksk>0y_k^Ts_k\gt0。

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

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

B_{k+1}= \begin{cases} B_k, & if \quad y_k^Ts_k\le0 \\ B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k} ,&if \quad y_k^Ts_k\gt0\end{cases}

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

算法

给定参数δ∈(0,1),σ∈(0,0.5)\delta \in(0,1),\sigma \in(0,0.5)，初始化点x0∈Rnx_0 \in R^n，终止误差0≤ϵ≪10 \le\epsilon \ll1,初始化对称正定阵B0B_0，通常取为G(xo)G(x_o)或单位阵InI_n;令k=0k=0。
计算gk=∇f(xk)g_k=\nabla f(x_k)，若∥gk∥≪ϵ\left \| g_k \right \| \ll \epsilon,终止，输出xkx_k作为近似极小点。
解线性方程组得解dkd_k:Bkd=−gkB_kd=-g_k.
设mkm_k是满足下列不等式的最小非负整数m:
f(xk+δmdk)≤f(xk)+σδmgTkdk

f(x_k+\delta^m d_k)\le f(x_k)+\sigma \delta^m g_k^Td_k
令αk=δmk,xk+1=xk+αkdk\alpha_k=\delta^{m_k},x_{k+1}=x_k+\alpha_kd_k.
由BFGS校正公式确定Bk+1B_{k+1}
令k=k+1k=k+1，转向步骤“2”

求解具体优化问题

求解无约束优化问题:

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

\min f(s)=100(x_1^2-x_2)^2+(x_1-1)^2,x=(x_1,x_2)^T\in R^2

#coding:UTF-8
'''
Created on 2017年4月20日@author: zhangdapeng
'''
from numpy import *
import matplotlib.pyplot as plt
from 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 = 0epsilon=1e-10while (k < maxk):  gk = mat(gfun(x0))#计算梯度 ，mat函数将数组转化为矩阵。
#         print(gk)
#         print(linalg.norm(gk,1))#axis=0,沿着纵轴方向if linalg.norm(gk,1)<epsilon:breakdk = mat(-linalg.solve(Bk, gk))  #解矩阵方程Bk*x=gk得到xm = 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/

拟牛顿法之BFGS算法相关推荐

最优化学习笔记(十九)——拟牛顿法(5)BFGS算法
一.BFGS算法的更新公式为了推导BFGS算法,需要用到对偶或者互补的概念,前边已经讨论过hessian矩阵逆矩阵的近似矩阵需要满足以下条件: Hk+1Δg(i)=Δx(i)0≤i≤k \bolds ...
优化算法——拟牛顿法之BFGS算法
一.BFGS算法简介 BFGS算法是使用较多的一种拟牛顿方法,是由Broyden,Fletcher,Goldfarb,Shanno四个人分别提出的,故称为BFGS校正. 同DFP校正的推导公式一样,D ...
c语言dfp算法程序,拟牛顿法中的DFP算法和BFGS算法
注明:程序中调用的函数jintuifa.m golddiv.m我在之前的笔记中已贴出 DFP算法和BFGS算法不同在于H矩阵的修正公式不同 DFP算法 %拟牛顿法中DFP算法求解f = x1*x1+2 ...
L-BFGS算法/Broyden族/BFGS算法/阻尼牛顿法的Python实现代码
下面定义了三个Python语言编写的函数:函数表达式fun,梯度向量gfun,和海森矩阵hess.这三个表达式在后面各个算法的实现中会用到. # 函数表达式fun fun = lambda x:100 ...
优化算法——拟牛顿法之L-BFGS算法
一.BFGS算法在"优化算法--拟牛顿法之BFGS算法"中,我们得到了BFGS算法的校正公式: 利用Sherman-Morrison公式可对上式进行变换,得到令,则得到: 二. ...
bfgs算法c语言,机器学习算法实现解析——liblbfgs之L-BFGS算法
在博文"优化算法--拟牛顿法之L-BFGS算法"中,已经对L-BFGS的算法原理做了详细的介绍,本文主要就开源代码liblbfgs重新回顾L-BFGS的算法原理以及具体的实现过程, ...
拟牛顿法/Quasi-Newton，DFP算法/Davidon-Fletcher-Powell，及BFGS算法/Broyden-Fletcher-Goldfarb-Shanno...
拟牛顿法/Quasi-Newton,DFP算法/Davidon-Fletcher-Powell,及BFGS算法/Broyden-Fletcher-Goldfarb-Shanno 转载须注明出处:htt ...
l bfgs算法java代码_优化算法——拟牛顿法之L-BFGS算法
一.BFGS算法 BFGS算法的校正公式: 利用Sherman-Morrison公式可对上式进行变换,得到令 ,则得到: 二.BGFS算法存在的问题在BFGS算法中.每次都要存储近似Hesse矩阵 ...
dfp matlab,MATLAB拟牛顿法之DFP与BFGS算法
DFP算法原理由于博主使用WPS编辑的文本,公式无法赋值粘贴,这里以截图的方法给出了推导过程.博主会上传该DOC文档. BFGS算法原理 matlab代码(DFP) syms x1 x2 f=@(x ...

拟牛顿法之BFGS算法