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

深碍√TFBOYSˉ_ 2022-05-27 12:20 310阅读 0赞

![70][]

![70 1][]

![70 2][]

![70 3][]

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

*  但当海森矩阵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更新规则比较简单

DFP算法的校正公式如下：

![70 4][]

DFP算法描述：

![70 5][]

DFP算法的代码实现

######################## 拟牛顿_DFP算法#######################
    
    #算法描述:
    '''
    已知DFP算法的校正公式,注意公式中的两个分数结构，分母((yk)T)Hkyk和((sk)T)yk是标量，
    
    分子Hkyk((yk)T)Hk和sk(sk)T是与Hk同型的矩阵，而且都是对称矩阵。若Hk正定，且((sk)T)yk>0
    则Hk+1也正定。
    
    当采用精确线搜索时，矩阵序列Hk的正定新条件((sk)T)yk>0可以被满足。
    但对于ArmijoArmijo搜索准则来说，不能满足这一条件，需要做修正
    '''
    
    #step 1: 给定参数δ∈(0,1),σ∈(0,0.5)，初始点x0∈R,终止误差0≤ϵ≪1.
    #初始对称正定矩阵H0(通常取为G(x0)^−1或单位矩阵I).令k←0 
    #step 2: 计算gk=∇f(xk). 若||gk||≤ϵ，停止迭代，输出x∗≈xk 
    #step 3: 计算搜索方向
    #               dk=−Hkgk,其中的Hk表示海森近似矩阵的逆矩阵
    #step 4: 记mk是满足下列不等式的最小非负整数m.
    #f(xk+β^mdk)≤f(xk)+δ*β^m*(gk)T*dk,为方便起见*代表点乘
    #step 5: 由校正公式确定Hk+1 
    #step 6： k←k+1,转 step 2
    
    
    def dfp(fun,gfun,hess,x0):
        #功能：用DFP族算法求解无约束问题：min fun(x)
        #输入：x0是初始点，fun,gfun分别是目标函数和梯度
        #输出：x,val分别是近似最优点和最优解,k是迭代次数
        maxk = 1e5#最大迭代次数
        rho = 0.55#非线性搜索中的B因子
        sigma = 0.4#非线性搜索中δ的因子
        epsilon = 1e-5
        k = 0#初始化迭代次数
        n = np.shape(x0)[0]#用于初始化海森单位矩阵和校正部分
        #海森矩阵可以初始化为单位矩阵
        Hk = np.linalg.inv(hess(x0)) #输入x0的海森逆矩阵,或者单位矩阵np.eye(n)
        #Hk =np.linalg.inv(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#用于返回搜索步长最小非负整数m
            while m < 20: # 用Armijo搜索求步长
                if fun(x0+rho**m*dk) < fun(x0)+sigma*rho**m*np.dot(gk,dk):
                    mk = m
                    break
                m += 1
            #DFP校正
            x = x0 + rho**mk*dk#更新x的值
            sk = x - x0#计算更新x后的sk
            yk = gfun(x) - gk#计算更新x后的yk
    
            if np.dot(sk,yk) > 0:#均为一维则内积运算
                Hy = np.dot(Hk,yk)#其中的一组Hy[ 0.3826141  15.30415407]
                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#分别是最优点坐标，最优值，迭代次数
    print(dfp(fun,gfun,hess,[20,10]))#采用Hess(x0)的逆迭代次数更小

#(array([1., 1.]), 2.1683502509487778e-18, 104)
    #print(dfp(fun,gfun,hess,[5,6]))#采用np.eye(n)的逆迭代次数更大,不一定收敛
    #(array([1.02087426, 1.04176415]), 0.0004533843736217651, 100000)
    #采用单位矩阵来初始化海森矩阵的逆,迭代次数明显增大

[70]: /images/20220527/293386437deb475bbf26fca1daf27a0f.png
[70 1]: /images/20220527/5018a96c545f421e9040ce2382da099d.png
[70 2]: /images/20220527/bf07ea48a03a4dadae888a238248433d.png
[70 3]: /images/20220527/ddd3da2677474624b09b8b174ad99d4f.png
[70 4]: /images/20220527/8e4109581fa04fac99cf8d9b3afb11f4.png
[70 5]: /images/20220527/107fdd26ae314300933bb301ce196a7a.png