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

- 日理万妓 2022-05-27 12:23 197阅读 0赞

![70][]

![70 1][]

![70 2][]

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

**正确参考如下图所示:**

![70 3][]

![70 4][]

补充几点：

![70 5][]

![70 6][]

![70 7][]

![70 8][]

######################## 拟牛顿_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)

[70]: /images/20220527/46fc70ae77904927852948f604c55101.png
[70 1]: /images/20220527/7d25685f71294540bd95fd824c53c6a1.png
[70 2]: /images/20220527/8ab99dbf8eb341c69c13597d2a1468b6.png
[70 3]: /images/20220527/358f6a3d1dbc47c1a9c711808bd9c574.png
[70 4]: /images/20220527/b16a2635b67346a8bbc2dd3069708368.png
[70 5]: /images/20220527/e3b35c1e562c44ba8b085de6b18766b3.png
[70 6]: /images/20220527/90ec2853f55f4fafbe9fe3b796feeda8.png
[70 7]: /images/20220527/870c14015a9445d793068f4b92b7343a.png
[70 8]: /images/20220527/599ddb24eff34094b419081cf74d170e.png