窃·格瓦拉

只有坚强起来,才能不丧失温柔 在家一个人很无聊,都没有朋友、女友玩 ,进了局子里面个个是人才,说话又好听,我超喜欢在里面。...

梯度-牛顿-拟牛顿优化算法和实现

1.要求解的问题

求解无约束非线性优化问题

minxR2f(x)=100(x21x2)2+(x11)2

该问题有精确解x=(1,1)T,f(x)=0

其梯度向量

g(x)=(400x1(x21x2)+2(x11),200(x21x2))

其海森矩阵
G(x)=[1200x21400x2+2400x1400x1200]

下面定义了三个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]])

2.线搜索技术和Armijo准则

线搜索技术是求解许多优化问题下降算法的基本组成部分,但精确线性搜索需要计算很多的函数值和梯度值,从而耗费较多资源。特别是当迭代点远离最优点时,精确线搜索技术通常不是有效和合理的。对于许多优化算法,其收敛速度并不依赖于精确搜索过程。因此既能保证目标函数具有可接受的下降量,又能使最终形成的迭代序列收敛的非精确先搜索越来越流行。

符号约定:

  • gk : f(xk),即第目标函数关于k次迭代值xk的导数。
  • Gk : G(xk)=2f(xk),即海森矩阵。
  • dk : 第k次迭代的优化方向。在最速下降算法中,有dk=gk
  • αk : 第k次迭代的步长因子,有xk+1=xk+αkdk

在精确线性搜索中,步长因子αk由下面的式子确定:

αk=argminαf(xk+αdk)

而对于非精确线性搜索,选取的αk只要使得目标函数f得到可接受的下降量,即

Δfk=f(xk)f(xk+αkdk)>0

Armijo准则用于非精确线性搜索中步长因子α的确定,内容如下:

Armijo准则:
已知当前位置xk和优化的方向dk ,参数β(0,1),δ(0,0.5).令步长因子αk=βmk,其中mk为满足下列不等式的最小非负整数m

f(xk+βmdk)f(xk)+δβmgTkdk

由此确定的下一个位置xk+1=xk+αkdk
NOTE: 上面的公式仅仅适用于梯度下降(Gradient Descent),对于梯度上升,则略有不用。首先公式变成了
f(xkβmdk)f(xk)δβmgTkdk,
然后确定下一个位置xk+1=xkαkdk

该准则在接下来的介绍的几个算法中多次使用。

3.最速下降法及其Python实现

算法描述

step 1 :选取初始点x0Rn,容许误差0ϵ1 .令 k1.
step 2 :计算gk=f(xk). 若||gk||ϵ,停止迭代,输出xk作为近似最优解。
step 3 :取方向dk=gk.
step 4 :由线搜索技术确定步长因子αk.
step 5 :令xk+1xk+αkdk,转step 2.

上述step 3中步长因子αk的确定既可以使用精确线搜索方法,也可以使用非精确线搜索方法,在理论上都能保证其全局收敛性。若采用精确线搜索方法,则αk满足

ddαf(xk+αdk)|α=αk=f(xk+αkdk)Tdk=0

从而有

f(xk+1)Tf(xk)=0

xk+1处的梯度与其前驱xk处的梯度是正交的,也就是说,迭代点序列所走的路线是锯齿形的,这造成其收敛速度很缓慢。如下图所示,其中绿色折线所显示的路线就是由最速下降法得到的,红色曲线是由共轭梯度下降法确定的,通过对比可以看出梯度下降法所走的路线是锯齿形的,经测试,其收敛速度非常慢。

来自于wiki百科

最速下降法的Python实现

def gradient(fun,gfun,x0):
    #最速下降法求解无约束问题
    # x0是初始点,fun和gfun分别是目标函数和梯度
    maxk = 5000
    rho = 0.5
    sigma = 0.4
    k = 0
    epsilon = 1e-5

    while k<maxk:
        gk = gfun(x0)
        dk = -gk
        if np.linalg.norm(dk) < epsilon:
            break
        m = 0
        mk = 0
        while m< 20:
            if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                mk = m
                break
            m += 1
        x0 += rho**mk*dk
        k += 1

    return x0,fun(x0),k

性能测试
这里写图片描述

可以看到大约有一半的样本的迭代次数要超过2000次。

4.牛顿法

牛顿法的基本思想是用迭代点xk处的一阶导数(梯度gk)和二阶倒数(海森矩阵Gk)对目标函数进行二次函数近似,然后把二次模型的极小点作为新的迭代点。

牛顿法推导

函数f(x)xk处的泰勒展开式的前三项得

T(f,xk,3)=fk+gTk(xxk)+12(xxk)TGk(xxk)

则其稳定点
T=gk+Gk(xxk)=0

Gk非奇异,那么解上面的线性方程(标记其解为xk+1)即得到牛顿法的迭代公式
xk+1=xkG1kgk

即优化方向为dk=G1kgk

求稳定点是用到了以下公式:

考虑y=xTAx,根据矩阵求导法则,有

dydx=Ax+ATxd2ydx2=A+AT

注意实际中dk的是通过求解线性方程组Gkd=gk获得的。

阻尼牛顿法及其Python实现

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

阻尼牛顿法的算法描述

step 1: 给定终止误差0ϵ1,δ(0,1),σ(0,0.5). 初始点x0Rn. 令k0
step 2: 计算gk=f(xk). 若||gk||ϵ,停止迭代,输出xxk
step 3: 计算Gk=2f(xk),并求解线性方程组得到解dk,

Gkd=gk

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

step 5: 令αk=δmk,xk+1=xk+αkdk,kk+1,转step 2

阻尼牛顿法的Python实现

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

    return 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)在每个迭代点xk处是正定的,否则难以保证牛顿方向dk=G1kgkfxk处的下降方向。为克服这一缺陷,需要对牛顿法进行修正。
下面介绍两种修正方法,分别是“牛顿-最速下降混合算法”和“修正牛顿法”。
牛顿-最速下降混合算法
该方法将牛顿法和最速下降法结合起来,基本思想是:当海森矩阵正定时,采用牛顿法确定的优化方向作为搜索方向;否则,即海森矩阵为非正定矩阵时,则采用负梯度方向作为搜索方向。

修正牛顿法
引入阻尼因子μk0,即在每一迭代步适当选取参数μk,使得矩阵Ak=Gk+μkI正定,用Ak代替Gk来求解dk

算法描述

step 1: 给定终止误差0ϵ1,δ(0,1),σ(0,0.5). 初始点x0Rn. 令k0
step 2: 计算gk=f(xk)μk=||gk||1+τ. 若||gk||ϵ,停止迭代,输出xxk
step 3: 计算Gk=2f(xk),并求解线性方程组得到解dk,

(Gk+μkI)d=gk

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

step 5: 令αk=δmk,xk+1=xk+αkdk,kk+1,转step 2

修正牛顿法的Python实现代码

def revisenm(fun,gfun,hess,x0):
    # 用修正牛顿法求解无约束问题
    #x0是初始点,fun,gfun和hess分别是目标函数值,梯度,海森矩阵的函数

    maxk = 1e5
    n = np.shape(x0)[0]
    rho = 0.55
    sigma = 0.4
    tau = 0.0
    k = 0
    epsilon = 1e-5

    while k < maxk:
        gk = gfun(x0)      
        if  np.linalg.norm(gk) < epsilon:
            break
        muk = 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 = 0
        mk = 0
        while m < 20:
            if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                mk = m
                break
            m += 1
        x0 += rho**mk*dk
        k += 1
    return x0,fun(x0),k

性能展示
这里写图片描述

5.拟牛顿法

牛顿法的优点是具有二阶收敛速度,缺点是:

  • 但当海森矩阵G(xk)=2f(x) 不正定时,不能保证所产生的方向是目标函数在xk处的下降方向。
  • 特别地,当G(xk)奇异时,算法就无法继续进行下去。尽管修正牛顿法可以克服这一缺陷,但修正参数的取值很难把握,过大或过小都会影响到收敛速度。
  • 牛顿法的每一步迭代都需要目标函数的海森矩阵G(xk),对于大规模问题其计算量是惊人的。

拟牛顿法的基本思想是用海森矩阵Gk的某个近似矩阵Bk取代Gk. Bk通常具有下面三个特点:

  • 在某种意义下有BkGk ,使得相应的算法产生的方向近似于牛顿方向,确保算法具有较快的收敛速度。
  • 对所有的kBk是正定的,从而使得算法所产生的方向是函数fxk处下降方向。
  • 矩阵Bk更新规则比较简单

f:RnR在开集DRn上二次连续可微,那么fxk+1处的二次近似模型为:

f(x)f(xk+1)+gTk+1(xxk+1)+12(xxk+1)TGk+1(xxk+1)

对上式求导得

g(x)gk+1+Gk+1(xxk+1)

x=xk,位移sk=xk+1xk,梯度差yk=gk+1gk,则有
Gk+1skyk
.

因此,拟牛顿法中近似矩阵Bk要满足关系式

Bk+1sk=yk

Hk+1=B1k+1,得到拟牛顿法的另一形式
Hk+1yk=sk

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

其中Ek,Dk为秩1或秩2矩阵。该公式称为校正规则
通常将上面的三个式子(拟牛顿方程和校正规则)所确立的方法称为拟牛顿法。从下面的拟牛顿法的几个变种可以看出,不同的拟牛顿法的主要区别在于更新公式的不同。

DFP算法及其Python实现

DFP算法的校正公式如下:

Hk+1=HkHkykyTkHkyTkHkyk+sksTksTkyk

注意公式中的两个分数结构,分母yTkHkyksTkyk是标量,分子HkykyTkHksksTk是与Hk同型的矩阵,而且都是对称矩阵。若Hk正定,且sTkyk>0,则Hk+1也正定。

当采用精确线搜索时,矩阵序列Hk的正定新条件sTkyk>0可以被满足。但对于Armijo搜索准则来说,不能满足这一条件,需要做如下修正:

Hk+1=HkHkHkykyTkHkyTkHkyk+sksTksTkyksTkyk0sTkyk>0

DFP算法描述

step 1: 给定参数δ(0,1),σ(0,0.5),初始点x0R,终止误差0ϵ1.初始对称正定矩阵H0(通常取为G(x0)1或单位矩阵In).令k0
step 2: 计算gk=f(xk). 若||gk||ϵ,停止迭代,输出xxk
step 3: 计算搜索方向

dk=Hkgk

step 4: 记mk是满足下列不等式的最小非负整数m.
f(xk+βmdk)f(xk)+δβmgTkdk
αk=δmk,xk+1=xk+αkdk
step 5: 由校正公式确定Hk+1
step 6kk+1,转 step 2

DFP算法的代码实现

def dfp(fun,gfun,hess,x0):
    #功能:用DFP族算法求解无约束问题:min fun(x)
    #输入:x0是初始点,fun,gfun分别是目标函数和梯度
    #输出:x,val分别是近似最优点和最优解,k是迭代次数
    maxk = 1e5
    rho = 0.55
    sigma = 0.4
    epsilon = 1e-5
    k = 0
    n = 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:
            break
        dk = -1.0*np.dot(Hk,gk)

        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
        sk = x - x0
        yk = gfun(x) - gk   

        if np.dot(sk,yk) > 0:
            Hy = np.dot(Hk,yk)
            print Hy
            sy = 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/sy

        k += 1
        x0 = x


    return 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=BkBkBkykyTkBkyTkBkyk+sksTksTkyksTkyk0sTkyk>0

BFGS算法的Python实现

def bfgs(fun,gfun,hess,x0):
    #功能:用BFGS族算法求解无约束问题:min fun(x)
    #输入:x0是初始点,fun,gfun分别是目标函数和梯度
    #输出:x,val分别是近似最优点和最优解,k是迭代次数  
    maxk = 1e5
    rho = 0.55
    sigma = 0.4
    epsilon = 1e-5
    k = 0
    n = 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:
            break
        dk = -1.0*np.linalg.solve(Bk,gk)
        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

        #BFGS校正
        x = x0 + rho**mk*dk
        sk = x - x0
        yk = 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/ys

        k += 1
        x0 = x

    return 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校正都是由ykBksk(或 skHkyk)组成的秩2矩阵。而Broyden族算法采用了BFGS和DFP校正公式的凸组合,如下:

Hϕk+1=ϕkHBFGSk+1+(1ϕk)HDFPk+1=HkHkykyTkHkyTkHkyk+sksTksTkyk+ϕkvkvTk

其中ϕk[0,1]vk由下式定义:
vk=yTkHkyk(skyTkskHkykyTkHkyk)

Broyden族算法Python实现

def broyden(fun,gfun,hess,x0):
    #功能:用Broyden族算法求解无约束问题:min fun(x)
    #输入:x0是初始点,fun,gfun分别是目标函数和梯度
    #输出:x,val分别是近似最优点和最优解,k是迭代次数
    x0 = np.array(x0)

    maxk = 1e5
    rho = 0.55;
    sigma = 0.4;
    epsilon = 1e-5
    phi = 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:
            break
        dk = -1*np.dot(Hk,gk)

        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
        #Broyden族校正
        x = x0 + rho**mk*dk
        sk = x - x0
        yk = gfun(x) - gk

        Hy = 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)*Hy
            sy = 0.2*yHy

        vk = 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))*vk
        k += 1
        x0 = x
    return x0,fun(x0),k #分别是最优点坐标,最优值,迭代次数

实际性能
这里写图片描述
相关代码:

iters = []
for i in xrange(np.shape(data)[0]):
    rt = broyden(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'Broyden迭代次数分布',{'fontname':'STFangsong','fontsize':18})
plt.xlabel(u'迭代次数',{'fontname':'STFangsong','fontsize':18})
plt.ylabel(u'频率分布',{'fontname':'STFangsong','fontsize':18})

L-BFGS算法及其Python实现

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

基本BFGS算法的校正公式可以改写成

Hk+1=(IskyTksTkyk)Hk(IyksTksTkyk)+sksTksTkyk

ρk=1/sTkyk,以及Vk=(IρkyksTk),则上式可以写成
Hk+1=VTkHkVk+ρksksTk

给定一个初始矩阵H0(其取值后面有提到),则由上式的递推关系可以得到

H1=VT0H0V0+ρ0s0sT0

H2=VT1H1V1+ρ1s1sT1=VT1(VT0H0V0+ρ0s0sT0)V1+ρ1s1sT1=VT1VT0H0V0V1+VT1ρ0s0sT0V1+ρ1s1sT1


更一般的,我们有

Hm+1=(VTKVTm1VT1VT0)H0(V0V1Vm1Vk)+(VTmVTm1VT2VT1)ρ0s0sT0(V1V2Vm1Vk)+(VTmVTm1VT3VT2)ρ1s1sT1(V2V3Vm1Vk)++(VTmVTm1)ρm2sm2sTm2(Vm1Vk)VTkρm1sm1sTm1Vk+ρmsmsTm

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

补充几点

  • H0=IsTmymyTmym,第一次迭代时,序列{sk}和{yk}为空,则H0=I
  • 最初的若干次迭代, 序列{sk}和{yk}的元素个数较小,会有n<m,则将Hm+1计算公式右边的m换成n即可。
  • 随着迭代次数增加, 序列{sk}和{yk}的元素个数增大,会有n>m。由于我们只需要最近m次迭代产生的序列元素,所以序列{sk}和{yk}只需要保存最新的m个元素即可,如果再有新的元素进入,则同时扔掉最旧的元素,以保证序列元素个数恒为m

这就是L-BFGS算法的思想。聪明的同学会问:你上面的公式不还是要计算海森矩阵的近似矩阵吗?那岂不是还是需要n2规模的内存?
其实在实际中是不计算该矩阵的, 而且不是采用上面的方法,而是直接得到Hkgk,而搜索方向就是dk=Hkgk。下面给出了计算的函数twoloop(s,y,ρ,gk)的伪代码,可知该函数内部的计算仅限于标量和向量,未出现矩阵。

函数参数s,y即向量序列{sk}和{yk},ρ为元素序列{ρk},其元素ρk=1/sTkykgk是向量,为当前的一阶导数,输出为向量z=Hkgk,即搜索方向的反方向

Functiontwoloop(s,y,ρ,gk)q=gkFori=k1,k2,,kmαi=ρisTiqq=qαiyiHk=yTk1sk1/yTk1yk1z=HkqFori=km,km+1,,k1βi=ρiyTizz=z+si(αiβi)Returnz

给出L-BFGS的算法

可以看到其搜索方向dk是根据函数twoloop计算得到的。
这里写图片描述

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 = 1

    a = 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*q

    for 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 = 2000
    rou = 0.55
    sigma = 0.4
    epsilon = 1e-5
    k = 0
    n = np.shape(x0)[0] #自变量的维度

    s, y, rho = [], [], []

    while k < maxk :
        gk = gfun(x0)
        if np.linalg.norm(gk) < epsilon:
            break

        dk = -1.0*twoloop(s, y, rho,gk)

        m0=0;
        mk=0
        while m0 < 20: # 用Armijo搜索求步长
            if fun(x0+rou**m0*dk) < fun(x0)+sigma*rou**m0*np.dot(gk,dk):
                mk = m0
                break
            m0 += 1


        x = x0 + rou**mk*dk
        sk = x - x0
        yk = 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 += 1
        x0 = x

    return 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})

参考文献

阅读更多
个人分类: ConvexOptimization
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭