下面定义了三个Python语言编写的函数：函数表达式fun，梯度向量gfun，和海森矩阵hess。这三个表达式在后面各个算法的实现中会用到。

# 函数表达式fun
fun = lambda x:100*(x[0]**2-x[1])**2 + (x[0]-1)**2# 梯度向量 gfun
gfun = lambda x:np.array([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1), -200*(x[0]**2-x[1])])# 海森矩阵 hess
hess = lambda x:np.array([[1200*x[0]**2-400*x[1]+2, -400*x[0]],[-400*x[0],200]])

最速下降法的Python实现。

def gradient(fun,gfun,x0):#最速下降法求解无约束问题# x0是初始点，fun和gfun分别是目标函数和梯度maxk = 5000rho = 0.5sigma = 0.4k = 0epsilon = 1e-5while k<maxk:gk = gfun(x0)dk = -gkif np.linalg.norm(dk) < epsilon:breakm = 0mk = 0while m< 20:if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1x0 += rho**mk*dkk += 1return x0,fun(x0),k

4.牛顿法
牛顿法的基本思想是用迭代点xkxk处的一阶导数（梯度gkgk）和二阶倒数（海森矩阵GkGk）对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点。

牛顿法推导

函数f(x)f(x)在xkxk处的泰勒展开式的前三项得
T(f,xk,3)=fk+gTk(x−xk)+12(x−xk)TGk(x−xk)
T(f,xk,3)=fk+gkT(x−xk)+12(x−xk)TGk(x−xk)

则其稳定点
∇T=gk+Gk(x−xk)=0
∇T=gk+Gk(x−xk)=0

若 GkGk非奇异，那么解上面的线性方程（标记其解为 xk+1xk+1）即得到牛顿法的迭代公式
xk+1=xk−G−1kgk
xk+1=xk−Gk−1gk

即优化方向为 dk=−G−1kgkdk=−Gk−1gk
求稳定点是用到了以下公式：

考虑y=xTAxy=xTAx，根据矩阵求导法则，有
dydx=Ax+ATxd2ydx2=A+AT
dydx=Ax+ATxd2ydx2=A+AT
注意实际中dkdk的是通过求解线性方程组Gkd=−gkGkd=−gk获得的。

阻尼牛顿法及其Python实现
阻尼牛顿法采用了牛顿法的思想，唯一的差别是阻尼牛顿法并不直接采用xk+1=xk+dkxk+1=xk+dk，而是用线搜索技术获得合适的步长因子αkαk,用公式xk+1=xk+δmdkxk+1=xk+δmdk计算xk+1xk+1。

阻尼牛顿法的算法描述：

step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5). 初始点x0∈Rnx0∈Rn. 令k←0k←0
step 2: 计算gk=∇f(xk)gk=∇f(xk). 若||gk||≤ϵ||gk||≤ϵ，停止迭代，输出x∗≈xkx∗≈xk
step 3: 计算Gk=∇2f(xk)Gk=∇2f(xk)，并求解线性方程组得到解dkdk,
Gkd=−gk
Gkd=−gk

step 4: 记 mkmk是满足下列不等式的最小非负整数 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk

step 5: 令 αk=δmk,xk+1=xk+αkdk,k←k+1αk=δmk,xk+1=xk+αkdk,k←k+1，转 step 2
阻尼牛顿法的Python实现：

def dampnm(fun,gfun,hess,x0):# 用阻尼牛顿法求解无约束问题#x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数maxk = 500rho = 0.55sigma = 0.4k = 0epsilon = 1e-5while k < maxk:gk = gfun(x0)Gk = hess(x0)dk = -1.0*np.linalg.solve(Gk,gk)if np.linalg.norm(dk) < epsilon:breakm = 0mk = 0while m < 20:if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1x0 += rho**mk*dkk += 1return x0,fun(x0),k

性能展示

作图代码：

iters = []
for i in xrange(np.shape(data)[0]):rt = dampnm(fun,gfun,hess,data[i])if rt[2] <=200:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=50)
plt.title(u'阻尼牛顿法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})

修正牛顿法法及其Python实现
牛顿法要求目标函数的海森矩阵G(x)=∇2f(x)G(x)=∇2f(x)在每个迭代点xkxk处是正定的，否则难以保证牛顿方向dk=−G−1kgkdk=−Gk−1gk是ff在xkxk处的下降方向。为克服这一缺陷，需要对牛顿法进行修正。
下面介绍两种修正方法，分别是“牛顿-最速下降混合算法”和“修正牛顿法”。
牛顿-最速下降混合算法
该方法将牛顿法和最速下降法结合起来，基本思想是：当海森矩阵正定时，采用牛顿法确定的优化方向作为搜索方向；否则，即海森矩阵为非正定矩阵时，则采用负梯度方向作为搜索方向。

修正牛顿法
引入阻尼因子μk≥0μk≥0，即在每一迭代步适当选取参数μkμk，使得矩阵Ak=Gk+μkIAk=Gk+μkI正定，用AkAk代替GkGk来求解dkdk。

算法描述：

step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5)0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5). 初始点x0∈Rnx0∈Rn. 令k←0k←0
step 2: 计算gk=∇f(xk)gk=∇f(xk)，μk=||gk||1+τμk=||gk||1+τ. 若||gk||≤ϵ||gk||≤ϵ，停止迭代，输出x∗≈xkx∗≈xk
step 3: 计算Gk=∇2f(xk)Gk=∇2f(xk)，并求解线性方程组得到解dkdk,
(Gk+μkI)d=−gk
(Gk+μkI)d=−gk

step 4: 记 mkmk是满足下列不等式的最小非负整数 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk

step 5: 令 αk=δmk,xk+1=xk+αkdk,k←k+1αk=δmk,xk+1=xk+αkdk,k←k+1，转 step 2
修正牛顿法的Python实现代码：

def revisenm(fun,gfun,hess,x0):# 用修正牛顿法求解无约束问题#x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数maxk = 1e5n = np.shape(x0)[0]rho = 0.55sigma = 0.4tau = 0.0k = 0epsilon = 1e-5while k < maxk:gk = gfun(x0)      if  np.linalg.norm(gk) < epsilon:breakmuk = np.power(np.linalg.norm(gk),1+tau)Gk = hess(x0)Ak = Gk + muk*np.eye(n)dk = -1.0*np.linalg.solve(Ak,gk)m = 0mk = 0while m < 20:if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1x0 += rho**mk*dkk += 1return x0,fun(x0),k

5.拟牛顿法
牛顿法的优点是具有二阶收敛速度，缺点是：

但当海森矩阵G(xk)=∇2f(x)G(xk)=∇2f(x) 不正定时，不能保证所产生的方向是目标函数在xkxk处的下降方向。
特别地，当G(xk)G(xk)奇异时，算法就无法继续进行下去。尽管修正牛顿法可以克服这一缺陷，但修正参数的取值很难把握，过大或过小都会影响到收敛速度。
牛顿法的每一步迭代都需要目标函数的海森矩阵G(xk)G(xk)，对于大规模问题其计算量是惊人的。
拟牛顿法的基本思想是用海森矩阵GkGk的某个近似矩阵BkBk取代GkGk. BkBk通常具有下面三个特点：

在某种意义下有Bk≈GkBk≈Gk ,使得相应的算法产生的方向近似于牛顿方向，确保算法具有较快的收敛速度。
对所有的kk，BkBk是正定的，从而使得算法所产生的方向是函数ff在xkxk处下降方向。
矩阵BkBk更新规则比较简单
设 f:Rn→Rf:Rn→R在开集D⊂RnD⊂Rn上二次连续可微，那么ff在xk+1xk+1处的二次近似模型为：
f(x)≈f(xk+1)+gTk+1(x−xk+1)+12(x−xk+1)TGk+1(x−xk+1)
f(x)≈f(xk+1)+gk+1T(x−xk+1)+12(x−xk+1)TGk+1(x−xk+1)
对上式求导得
g(x)≈gk+1+Gk+1(x−xk+1)
g(x)≈gk+1+Gk+1(x−xk+1)

令 x=xkx=xk，位移 sk=xk+1−xksk=xk+1−xk，梯度差 yk=gk+1−gkyk=gk+1−gk，则有
Gk+1sk≈yk
Gk+1sk≈yk
.
因此，拟牛顿法中近似矩阵BkBk要满足关系式
Bk+1sk=yk
Bk+1sk=yk

令 Hk+1=B−1k+1Hk+1=Bk+1−1，得到拟牛顿法的另一形式
Hk+1yk=sk
Hk+1yk=sk

其中 Hk+1Hk+1为海森矩阵逆矩阵的近似。上面两个公式称为拟牛顿方程。
搜索方向由 dk=−Hkgkdk=−Hkgk或 Bkdk=−gkBkdk=−gk确定。根据近似矩阵的第三个特点，有
Bk+1=Bk+Ek,Hk+1=Hk+Dk
Bk+1=Bk+Ek,Hk+1=Hk+Dk

其中 Ek,DkEk,Dk为秩1或秩2矩阵。该公式称为校正规则。
通常将上面的三个式子（拟牛顿方程和校正规则）所确立的方法称为拟牛顿法。从下面的拟牛顿法的几个变种可以看出，不同的拟牛顿法的主要区别在于更新公式的不同。
DFP算法及其Python实现
DFP算法的校正公式如下：
Hk+1=Hk−HkykyTkHkyTkHkyk+sksTksTkyk
Hk+1=Hk−HkykykTHkykTHkyk+skskTskTyk
注意公式中的两个分数结构，分母yTkHkykykTHkyk和sTkykskTyk是标量，分子HkykyTkHkHkykykTHk和sksTkskskT是与HkHk同型的矩阵，而且都是对称矩阵。若HkHk正定，且sTkyk>0skTyk>0，则Hk+1Hk+1也正定。

当采用精确线搜索时，矩阵序列HkHk的正定新条件sTkyk>0skTyk>0可以被满足。但对于ArmijoArmijo搜索准则来说，不能满足这一条件，需要做如下修正：
Hk+1=⎧⎩⎨HkHk−HkykyTkHkyTkHkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬
Hk+1={HkskTyk≤0Hk−HkykykTHkykTHkyk+skskTskTykskTyk>0}
DFP算法描述：

step 1: 给定参数δ∈(0,1),σ∈(0,0.5)δ∈(0,1),σ∈(0,0.5)，初始点x0∈Rx0∈R,终止误差0≤ϵ≪10≤ϵ≪1.初始对称正定矩阵H0H0(通常取为G(x0)−1G(x0)−1或单位矩阵InIn).令k←0k←0
step 2: 计算gk=∇f(xk)gk=∇f(xk). 若||gk||≤ϵ||gk||≤ϵ，停止迭代，输出x∗≈xkx∗≈xk
step 3: 计算搜索方向
dk=−Hkgk
dk=−Hkgk

step 4: 记 mkmk是满足下列不等式的最小非负整数 mm.
f(xk+βmdk)≤f(xk)+δβmgTkdk
f(xk+βmdk)≤f(xk)+δβmgkTdk
令 αk=δmk,xk+1=xk+αkdkαk=δmk,xk+1=xk+αkdk
step 5: 由校正公式确定 Hk+1Hk+1
step 6： k←k+1k←k+1,转 step 2
DFP算法的代码实现

def dfp(fun,gfun,hess,x0):#功能：用DFP族算法求解无约束问题：min fun(x)#输入：x0是初始点，fun,gfun分别是目标函数和梯度#输出：x,val分别是近似最优点和最优解,k是迭代次数maxk = 1e5rho = 0.55sigma = 0.4epsilon = 1e-5k = 0n = np.shape(x0)[0]#海森矩阵可以初始化为单位矩阵Hk = np.linalg.inv(hess(x0)) #或者单位矩阵np.eye(n)while k < maxk :gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1.0*np.dot(Hk,gk)m=0;mk=0while m < 20: # 用Armijo搜索求步长if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1#print mk#DFP校正x = x0 + rho**mk*dksk = x - x0yk = gfun(x) - gk   if np.dot(sk,yk) > 0:Hy = np.dot(Hk,yk)print Hysy = np.dot(sk,yk) #向量的点积yHy = np.dot(np.dot(yk,Hk),yk) # yHy是标量#表达式Hy.reshape((n,1))*Hy 中Hy是向量，生成矩阵Hk = Hk - 1.0*Hy.reshape((n,1))*Hy/yHy + 1.0*sk.reshape((n,1))*sk/syk += 1x0 = xreturn x0,fun(x0),k#分别是最优点坐标，最优值，迭代次数

由以上代码可知，拟牛顿法只需要在迭代开始前计算一次海森矩阵，更一般的，根本就不计算海森矩阵，而是初始化为单位矩阵，在每次迭代过程中是不需按传统方法（二阶导数）计算海森矩阵的。

实际性能

相关代码：

n = 50
x = y = np.linspace(-10,10,n)
xx,yy = np.meshgrid(x,y)
data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]
iters = []
for i in xrange(np.shape(data)[0]):rt = dfp(fun,gfun,hess,np.array(data[i]))if rt[2] <=200: # 提出迭代次数过大的iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=50)
plt.title(u'DFP迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})

BFGS算法及其Python实现
BFGS算法与DFP步骤基本相同，区别在于更新公式的差异
Bk+1=⎧⎩⎨BkBk−BkykyTkBkyTkBkyk+sksTksTkyksTkyk≤0sTkyk>0⎫⎭⎬
Bk+1={BkskTyk≤0Bk−BkykykTBkykTBkyk+skskTskTykskTyk>0}
BFGS算法的Python实现

def bfgs(fun,gfun,hess,x0):#功能：用BFGS族算法求解无约束问题：min fun(x)#输入：x0是初始点，fun,gfun分别是目标函数和梯度#输出：x,val分别是近似最优点和最优解,k是迭代次数  maxk = 1e5rho = 0.55sigma = 0.4epsilon = 1e-5k = 0n = np.shape(x0)[0]#海森矩阵可以初始化为单位矩阵Bk = eye(n)#np.linalg.inv(hess(x0)) #或者单位矩阵np.eye(n)while k < maxk:gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1.0*np.linalg.solve(Bk,gk)m = 0mk = 0while m < 20: # 用Armijo搜索求步长if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1#BFGS校正x = x0 + rho**mk*dksk = x - x0yk = gfun(x) - gk   if np.dot(sk,yk) > 0:    Bs = np.dot(Bk,sk)ys = np.dot(yk,sk)sBs = np.dot(np.dot(sk,Bk),sk) Bk = Bk - 1.0*Bs.reshape((n,1))*Bs/sBs + 1.0*yk.reshape((n,1))*yk/ysk += 1x0 = xreturn x0,fun(x0),k#分别是最优点坐标，最优值，迭代次数

实际性能

相关代码：

iters = []
for i in xrange(np.shape(data)[0]):rt = bfgs(fun,gfun,hess,np.array(data[i]))if rt[2] <=200:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=50)
plt.title(u'BFGS迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})

Broyden族算法及其Python实现
之前的BFGS和DFP校正都是由ykyk和BkskBksk(或 sksk和HkykHkyk)组成的秩2矩阵。而Broyden族算法采用了BFGS和DFP校正公式的凸组合，如下：
Hϕk+1=ϕkHBFGSk+1+(1−ϕk)HDFPk+1=Hk−HkykyTkHkyTkHkyk+sksTksTkyk+ϕkvkvTk
Hk+1ϕ=ϕkHk+1BFGS+(1−ϕk)Hk+1DFP=Hk−HkykykTHkykTHkyk+skskTskTyk+ϕkvkvkT

其中 ϕk∈[0,1]ϕk∈[0,1]， vkvk由下式定义：
vk=yTkHkyk−−−−−−√(skyTksk−HkykyTkHkyk)
vk=ykTHkyk(skykTsk−HkykykTHkyk)
Broyden族算法Python实现

def broyden(fun,gfun,hess,x0):#功能：用Broyden族算法求解无约束问题：min fun(x)#输入：x0是初始点，fun,gfun分别是目标函数和梯度#输出：x,val分别是近似最优点和最优解,k是迭代次数x0 = np.array(x0)maxk = 1e5rho = 0.55;sigma = 0.4;epsilon = 1e-5phi = 0.5;k=0;n = np.shape(x0)[0]Hk = np.linalg.inv(hess(x0))while k<maxk :gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1*np.dot(Hk,gk)m=0;mk=0while m < 20: # 用Armijo搜索求步长if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):mk = mbreakm += 1#Broyden族校正x = x0 + rho**mk*dksk = x - x0yk = gfun(x) - gkHy = np.dot(Hk,yk)sy = np.dot(sk,yk)yHy = np.dot(np.dot(yk,Hk),yk)if(sy < 0.2 *yHy):theta = 0.8*yHy/(yHy-sy)sk = theta*sk + (1-theta)*Hysy = 0.2*yHyvk = np.sqrt(yHy)*(sk/sy-Hy/yHy)Hk = Hk - Hy.reshape((n,1))*Hy/yHy +sk.reshape((n,1))*sk/sy + phi*vk.reshape((n,1))*vkk += 1x0 = xreturn x0,fun(x0),k #分别是最优点坐标，最优值，迭代次数

实际性能

相关代码：

n = 50
x = y = np.linspace(-10,10,n)
xx,yy = np.meshgrid(x,y)
data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]iters = []
for i in xrange(np.shape(data)[0]):rt = lbfgs(fun,gfun,data[i])if rt[2] <=1000:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=100)
plt.title(u'L-BFGS法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})

L-BFGS算法及其Python实现
拟牛顿法(如BFGS算法)需要计算和存储海森矩阵，其空间复杂度是n2n2，当nn很大时，其需要的内存量是非常大的。为了解决该内存问题，有限内存BFGS（即传说中的L-BFGS算法）横空出世。

基本BFGS算法的校正公式可以改写成
Hk+1=(I−skyTksTkyk)Hk(I−yksTksTkyk)+sksTksTkyk
Hk+1=(I−skykTskTyk)Hk(I−ykskTskTyk)+skskTskTyk

记ρk=1/sTkykρk=1/skTyk，以及Vk=(I−ρkyksTk)Vk=(I−ρkykskT)，则上式可以写成
Hk+1=VTkHkVk+ρksksTk
Hk+1=VkTHkVk+ρkskskT
给定一个初始矩阵H0H0(其取值后面有提到)，则由上式的递推关系可以得到

H1=VT0H0V0+ρ0s0sT0
H1=V0TH0V0+ρ0s0s0T
H2=VT1H1V1+ρ1s1sT1=VT1(VT0H0V0+ρ0s0sT0)V1+ρ1s1sT1=VT1VT0H0V0V1+VT1ρ0s0sT0V1+ρ1s1sT1
H2=V1TH1V1+ρ1s1s1T=V1T(V0TH0V0+ρ0s0s0T)V1+ρ1s1s1T=V1TV0TH0V0V1+V1Tρ0s0s0TV1+ρ1s1s1T
⋯⋯⋯
⋯⋯⋯

更一般的，我们有
Hm+1=(VTKVTm−1⋯VT1VT0)H0(V0V1⋯Vm−1Vk)+(VTmVTm−1⋯VT2VT1)ρ0s0sT0(V1V2⋯Vm−1Vk)+(VTmVTm−1⋯VT3VT2)ρ1s1sT1(V2V3⋯Vm−1Vk)+⋯+(VTmVTm−1)ρm−2sm−2sTm−2(Vm−1Vk)VTkρm−1sm−1sTm−1Vk+ρmsmsTm
Hm+1=(VKTVm−1T⋯V1TV0T)H0(V0V1⋯Vm−1Vk)+(VmTVm−1T⋯V2TV1T)ρ0s0s0T(V1V2⋯Vm−1Vk)+(VmTVm−1T⋯V3TV2T)ρ1s1s1T(V2V3⋯Vm−1Vk)+⋯+(VmTVm−1T)ρm−2sm−2sm−2T(Vm−1Vk)VkTρm−1sm−1sm−1TVk+ρmsmsmT

在上式的右边中， H0H0是由我们设置的，其余变量可由保存的最近 mm次迭代所形成的向量序列，如位移向量序列{sk}{sk}和一阶导数差所形成的向量序列 {yk}{yk}获得。也就是说，可由最近 mm次迭代的信息计算当前的海森矩阵的近似矩阵。
补充几点：

H0=I⋅sTmymyTmymH0=I⋅smTymymTym，第一次迭代时，序列{ sksk}和{ ykyk}为空，则 H0=IH0=I
最初的若干次迭代，序列{sksk}和{ykyk}的元素个数较小，会有n<mn<m，则将Hm+1Hm+1计算公式右边的mm换成nn即可。
随着迭代次数增加，序列{sksk}和{ykyk}的元素个数增大，会有n>mn>m。由于我们只需要最近mm次迭代产生的序列元素，所以序列{sksk}和{ykyk}只需要保存最新的mm个元素即可，如果再有新的元素进入，则同时扔掉最旧的元素，以保证序列元素个数恒为mm。
这就是L-BFGS算法的思想。聪明的同学会问：你上面的公式不还是要计算海森矩阵的近似矩阵吗？那岂不是还是需要n2n2规模的内存？
其实在实际中是不计算该矩阵的, 而且不是采用上面的方法，而是直接得到HkgkHkgk，而搜索方向就是dk=−Hkgkdk=−Hkgk。下面给出了计算的函数twoloop(s,y,ρρ,gk)的伪代码，可知该函数内部的计算仅限于标量和向量，未出现矩阵。

函数参数s,ys,y即向量序列{sksk}和{ykyk}，ρρ为元素序列{ρkρk}，其元素ρk=1/sTkykρk=1/skTyk，gkgk是向量，为当前的一阶导数，输出为向量z=Hkgkz=Hkgk，即搜索方向的反方向
Functiontwoloop(s,y,ρ,gk)q=gkFori=k−1,k−2,…,k−mαi=ρisTiqq=q−αiyiHk=yTk−1sk−1/yTk−1yk−1z=HkqFori=k−m,k−m+1,…,k−1βi=ρiyTizz=z+si(αi−βi)Returnz
Functiontwoloop(s,y,ρ,gk)q=gkFori=k−1,k−2,…,k−mαi=ρisiTqq=q−αiyiHk=yk−1Tsk−1/yk−1Tyk−1z=HkqFori=k−m,k−m+1,…,k−1βi=ρiyiTzz=z+si(αi−βi)Returnz
给出L-BFGS的算法

可以看到其搜索方向dkdk是根据函数twolooptwoloop计算得到的。

L-BFGS算法Python实现

def twoloop(s, y, rho,gk):n = len(s) #向量序列的长度if np.shape(s)[0] >= 1:#h0是标量，而非矩阵h0 = 1.0*np.dot(s[-1],y[-1])/np.dot(y[-1],y[-1])else:h0 = 1a = empty((n,))q = gk.copy() for i in range(n - 1, -1, -1): a[i] = rho[i] * dot(s[i], q)q -= a[i] * y[i]z = h0*qfor i in range(n):b = rho[i] * dot(y[i], z)z += s[i] * (a[i] - b)return z   def lbfgs(fun,gfun,x0,m=5):# fun和gfun分别是目标函数及其一阶导数,x0是初值,m为储存的序列的大小maxk = 2000rou = 0.55sigma = 0.4epsilon = 1e-5k = 0n = np.shape(x0)[0] #自变量的维度s, y, rho = [], [], []while k < maxk :gk = gfun(x0)if np.linalg.norm(gk) < epsilon:breakdk = -1.0*twoloop(s, y, rho,gk)m0=0;mk=0while m0 < 20: # 用Armijo搜索求步长if fun(x0+rou**m0*dk) < fun(x0)+sigma*rou**m0*np.dot(gk,dk):mk = m0breakm0 += 1x = x0 + rou**mk*dksk = x - x0yk = gfun(x) - gk   if np.dot(sk,yk) > 0: #增加新的向量rho.append(1.0/np.dot(sk,yk))s.append(sk)y.append(yk)if np.shape(rho)[0] > m: #弃掉最旧向量rho.pop(0)s.pop(0)y.pop(0)k += 1x0 = xreturn x0,fun(x0),k#分别是最优点坐标，最优值，迭代次数

相关代码：

n = 50
x = y = np.linspace(-10,10,n)
xx,yy = np.meshgrid(x,y)
data = [[xx[i][j],yy[i][j]] for i in range(n) for j in range(n)]iters = []
for i in xrange(np.shape(data)[0]):rt = lbfgs(fun,gfun,data[i])if rt[2] <=1000:iters.append(rt[2])if i%100 == 0:print i,rt[2]plt.hist(iters,bins=100)
plt.title(u'L-BFGS法迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})