机器学习算法__1__牛顿法,拟牛顿法,DFP,BFGS,L-BFGS 原理及代码详解(4)

此上最后部分LBFGS算法的展开式有疑议

正确参考如下图所示:

补充几点：

######################## 拟牛顿_L-BFGS算法#######################

'''
拟牛顿法(如BFGS算法)需要计算和存储海森矩阵，其空间复杂度是n2，当n很大时，
其需要的内存量是非常大的。为了解决该内存问题，有限内存BFGS（即传说中的L-BFGS算法）横空出世。
H0 是由我们设置的，其余变量可由保存的最近m次迭代所形成的向量序列，
如位移向量序列{sk}和一阶导数差所形成的向量序列{yk}获得。也就是说，
可由最近m次迭代的信息计算当前的海森矩阵的近似矩阵。
'''

def twoloop(s, y, rho,gk):
    
    '''
    函数参数s,ys,y即向量序列{sk}和{yk}，
    ρ为元素序列{ρk}，其元素ρk=1/((sk)T)yk，gk是向量，
    为当前的一阶导数，输出为向量z=Hkgk，即搜索方向的反方向 

    '''

    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 = np.empty((n,))

    q = gk.copy() 
    for i in range(n - 1, -1, -1): 
        a[i] = rho[i] * np.dot(s[i], q)
        q -= a[i] * y[i]
    z = h0*q

    for i in range(n):
        b = rho[i] * np.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#非线性搜索中的B因子
    sigma = 0.4#非线性搜索中的δ因子
    epsilon = 1e-5#设定迭代终止得的阈值
    k = 0#初始化迭代次数
    n = np.shape(x0)[0] #自变量的维度

    s, y, rho = [], [], []#s是sk序列的保存列表,y是yk序列的保存列表

    while k < maxk :
        gk = gfun(x0)#计算得到gk
        if np.linalg.norm(gk) < epsilon:
            break

        dk = -1.0*twoloop(s, y, rho,gk)#由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#分别是最优点坐标，最优值，迭代次数

print(lbfgs(fun,gfun,[10,20],m=5))
#(array([1.00000016, 1.00000031]), 5.273891558003461e-14, 55)