学生作业 制作一个网站,中国建设造价工程协会网站,两栏式设计网站,淄博网站建设排行榜下面定义了三个Python语言编写的函数#xff1a;函数表达式fun#xff0c;梯度向量gfun#xff0c;和海森矩阵hess。这三个表达式在后面各个算法的实现中会用到。
# 函数表达式fun
fun lambda x:100*(x[0]**2-x[1])**2 (x[0]-1)**2# 梯度向量 gfun
gfun lambda x:np.arr…下面定义了三个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 kmaxk:gk gfun(x0)dk -gkif np.linalg.norm(dk) epsilon:breakm 0mk 0while m 20:if fun(x0rho**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)fkgTk(x−xk)12(x−xk)TGk(x−xk) T(f,xk,3)fkgkT(x−xk)12(x−xk)TGk(x−xk)
则其稳定点 ∇TgkGk(x−xk)0 ∇TgkGk(x−xk)0
若 GkGk非奇异那么解上面的线性方程标记其解为 xk1xk1即得到牛顿法的迭代公式 xk1xk−G−1kgk xk1xk−Gk−1gk
即优化方向为 dk−G−1kgkdk−Gk−1gk 求稳定点是用到了以下公式
考虑yxTAxyxTAx根据矩阵求导法则有 dydxAxATxd2ydx2AAT dydxAxATxd2ydx2AAT 注意实际中dkdk的是通过求解线性方程组Gkd−gkGkd−gk获得的。
阻尼牛顿法及其Python实现 阻尼牛顿法采用了牛顿法的思想唯一的差别是阻尼牛顿法并不直接采用xk1xkdkxk1xkdk而是用线搜索技术获得合适的步长因子αkαk,用公式xk1xkδmdkxk1xkδmdk计算xk1xk1。
阻尼牛顿法的算法描述
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,xk1xkαkdk,k←k1αkδmk,xk1xkαkdk,k←k1转 step 2 阻尼牛顿法的Python实现
def dampnm(fun,gfun,hess,x0):# 用阻尼牛顿法求解无约束问题#x0是初始点fungfun和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(x0rho**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,bins50)
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使得矩阵AkGkμkIAkGkμ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,xk1xkαkdk,k←k1αkδmk,xk1xkαkdk,k←k1转 step 2 修正牛顿法的Python实现代码
def revisenm(fun,gfun,hess,x0):# 用修正牛顿法求解无约束问题#x0是初始点fungfun和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),1tau)Gk hess(x0)Ak Gk muk*np.eye(n)dk -1.0*np.linalg.solve(Ak,gk)m 0mk 0while m 20:if fun(x0rho**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 ,使得相应的算法产生的方向近似于牛顿方向确保算法具有较快的收敛速度。 对所有的kkBkBk是正定的从而使得算法所产生的方向是函数ff在xkxk处下降方向。 矩阵BkBk更新规则比较简单 设 f:Rn→Rf:Rn→R在开集D⊂RnD⊂Rn上二次连续可微那么ff在xk1xk1处的二次近似模型为 f(x)≈f(xk1)gTk1(x−xk1)12(x−xk1)TGk1(x−xk1) f(x)≈f(xk1)gk1T(x−xk1)12(x−xk1)TGk1(x−xk1) 对上式求导得 g(x)≈gk1Gk1(x−xk1) g(x)≈gk1Gk1(x−xk1)
令 xxkxxk位移 skxk1−xkskxk1−xk梯度差 ykgk1−gkykgk1−gk则有 Gk1sk≈yk Gk1sk≈yk . 因此拟牛顿法中近似矩阵BkBk要满足关系式 Bk1skyk Bk1skyk
令 Hk1B−1k1Hk1Bk1−1得到拟牛顿法的另一形式 Hk1yksk Hk1yksk
其中 Hk1Hk1为海森矩阵逆矩阵的近似。上面两个公式称为 拟牛顿方程。 搜索方向由 dk−Hkgkdk−Hkgk或 Bkdk−gkBkdk−gk确定。根据近似矩阵的第三个特点有 Bk1BkEk,Hk1HkDk Bk1BkEk,Hk1HkDk
其中 Ek,DkEk,Dk为秩1或秩2矩阵。该公式称为 校正规则。 通常将上面的三个式子拟牛顿方程和校正规则所确立的方法称为拟牛顿法。从下面的拟牛顿法的几个变种可以看出不同的拟牛顿法的主要区别在于更新公式的不同。 DFP算法及其Python实现 DFP算法的校正公式如下 Hk1Hk−HkykyTkHkyTkHkyksksTksTkyk Hk1Hk−HkykykTHkykTHkykskskTskTyk 注意公式中的两个分数结构分母yTkHkykykTHkyk和sTkykskTyk是标量分子HkykyTkHkHkykykTHk和sksTkskskT是与HkHk同型的矩阵而且都是对称矩阵。若HkHk正定且sTkyk0skTyk0则Hk1Hk1也正定。
当采用精确线搜索时矩阵序列HkHk的正定新条件sTkyk0skTyk0可以被满足。但对于ArmijoArmijo搜索准则来说不能满足这一条件需要做如下修正 Hk1⎧⎩⎨HkHk−HkykyTkHkyTkHkyksksTksTkyksTkyk≤0sTkyk0⎫⎭⎬ Hk1{HkskTyk≤0Hk−HkykykTHkykTHkykskskTskTykskTyk0} 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,xk1xkαkdkαkδmk,xk1xkαkdk step 5: 由校正公式确定 Hk1Hk1 step 6 k←k1k←k1,转 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)m0;mk0while m 20: # 用Armijo搜索求步长if fun(x0rho**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,bins50)
plt.title(uDFP迭代次数分布,{fontname:STFangsong,fontsize:18})
plt.xlabel(u迭代次数,{fontname:STFangsong,fontsize:18})
plt.ylabel(u频率分布,{fontname:STFangsong,fontsize:18}) BFGS算法及其Python实现 BFGS算法与DFP步骤基本相同区别在于更新公式的差异 Bk1⎧⎩⎨BkBk−BkykyTkBkyTkBkyksksTksTkyksTkyk≤0sTkyk0⎫⎭⎬ Bk1{BkskTyk≤0Bk−BkykykTBkykTBkykskskTskTykskTyk0} 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(x0rho**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,bins50)
plt.title(uBFGS迭代次数分布,{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ϕk1ϕkHBFGSk1(1−ϕk)HDFPk1Hk−HkykyTkHkyTkHkyksksTksTkykϕkvkvTk Hk1ϕϕkHk1BFGS(1−ϕk)Hk1DFPHk−HkykykTHkykTHkykskskTskTykϕkvkvkT
其中 ϕk∈[0,1]ϕk∈[0,1] vkvk由下式定义 vkyTkHkyk−−−−−−√(skyTksk−HkykyTkHkyk) vkykTHkyk(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;k0;n np.shape(x0)[0]Hk np.linalg.inv(hess(x0))while kmaxk :gk gfun(x0)if np.linalg.norm(gk) epsilon:breakdk -1*np.dot(Hk,gk)m0;mk0while m 20: # 用Armijo搜索求步长if fun(x0rho**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,bins100)
plt.title(uL-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算法的校正公式可以改写成 Hk1(I−skyTksTkyk)Hk(I−yksTksTkyk)sksTksTkyk Hk1(I−skykTskTyk)Hk(I−ykskTskTyk)skskTskTyk
记ρk1/sTkykρk1/skTyk以及Vk(I−ρkyksTk)Vk(I−ρkykskT)则上式可以写成 Hk1VTkHkVkρksksTk Hk1VkTHkVkρkskskT 给定一个初始矩阵H0H0(其取值后面有提到)则由上式的递推关系可以得到
H1VT0H0V0ρ0s0sT0 H1V0TH0V0ρ0s0s0T H2VT1H1V1ρ1s1sT1VT1(VT0H0V0ρ0s0sT0)V1ρ1s1sT1VT1VT0H0V0V1VT1ρ0s0sT0V1ρ1s1sT1 H2V1TH1V1ρ1s1s1TV1T(V0TH0V0ρ0s0s0T)V1ρ1s1s1TV1TV0TH0V0V1V1Tρ0s0s0TV1ρ1s1s1T ⋯⋯⋯ ⋯⋯⋯
更一般的我们有 Hm1(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 Hm1(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次迭代的信息计算当前的海森矩阵的近似矩阵。 补充几点
H0I⋅sTmymyTmymH0I⋅smTymymTym第一次迭代时序列{ sksk}和{ ykyk}为空则 H0IH0I 最初的若干次迭代 序列{sksk}和{ykyk}的元素个数较小会有nmnm则将Hm1Hm1计算公式右边的mm换成nn即可。 随着迭代次数增加 序列{sksk}和{ykyk}的元素个数增大会有nmnm。由于我们只需要最近mm次迭代产生的序列元素所以序列{sksk}和{ykyk}只需要保存最新的mm个元素即可如果再有新的元素进入则同时扔掉最旧的元素以保证序列元素个数恒为mm。 这就是L-BFGS算法的思想。聪明的同学会问你上面的公式不还是要计算海森矩阵的近似矩阵吗那岂不是还是需要n2n2规模的内存 其实在实际中是不计算该矩阵的, 而且不是采用上面的方法而是直接得到HkgkHkgk而搜索方向就是dk−Hkgkdk−Hkgk。下面给出了计算的函数twoloop(s,y,ρρ,gk)的伪代码可知该函数内部的计算仅限于标量和向量未出现矩阵。
函数参数s,ys,y即向量序列{sksk}和{ykyk}ρρ为元素序列{ρkρk}其元素ρk1/sTkykρk1/skTykgkgk是向量为当前的一阶导数输出为向量zHkgkzHkgk即搜索方向的反方向 Functiontwoloop(s,y,ρ,gk)qgkForik−1,k−2,…,k−mαiρisTiqqq−αiyiHkyTk−1sk−1/yTk−1yk−1zHkqForik−m,k−m1,…,k−1βiρiyTizzzsi(αi−βi)Returnz Functiontwoloop(s,y,ρ,gk)qgkForik−1,k−2,…,k−mαiρisiTqqq−αiyiHkyk−1Tsk−1/yk−1Tyk−1zHkqForik−m,k−m1,…,k−1βiρiyiTzzzsi(α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,m5):# 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)m00;mk0while m0 20: # 用Armijo搜索求步长if fun(x0rou**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,bins100)
plt.title(uL-BFGS法迭代次数分布,{fontname:STFangsong,fontsize:18})
plt.xlabel(u迭代次数,{fontname:STFangsong,fontsize:18})
plt.ylabel(u频率分布,{fontname:STFangsong,fontsize:18})
如果对代码有疑问欢迎添加微信 18565453898
