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

青旅半醒 2022-05-27 10:08 308阅读 0赞

**一.牛顿法**

Taylor公式 – Maclaurin公式:

![70][]

![70 1][]

1.求根问题:

牛顿法的最初提出是用来求解方程的根的。我们假设点x∗为函数f(x)的根，那么有f(x∗)=0。现在我们把函数f(x)在点xk处一阶泰勒展开有：

![70 2][]

那么假设点xk\+1为该方程的根，则有

![70 3][]

那么就可以得到

![70 4][]

这样我们就得到了一个递归方程，我们可以通过迭代的方式不断的让x趋近于x∗从而求得方程f(x)的解。该递归式同样可以通过下方式得到：

利用三角形特征可以知道。其中，

![70 5][]

f′(xn)在三角形中表示点(xn,f(xn))切线的斜率。

2.优化问题:

对于最优化问题，其极值点处有一个特性就是在极值点处函数的一阶导数为0。因此我们可以在一阶导数处利用牛顿法通过迭代的方式来求得最优解，即相当于求一阶导数对应函数的根。   
首先，我们对函数在xkxk点处进行二阶泰勒展开

![70 6][]

因此，当x→xk时

![70 7][]

这里假设点xk\+1是一阶导数的根，那么就有

![70 8][]

依据上式可以得到

![70 9][]

这样我们就得到了一个不断更新xx迭代求得最优解的方法。这个也很好理解，假设我们上面的第一张图的曲线表示的是函数f(x)f(x)一阶导数的曲线，那么其二阶导数就是一阶导数对应函数在某点的斜率，也就是那条切线的斜率，那么该公式就和上面求根的公式本质是一样的。

我们这里讨论的都是在低维度的情形下，那么对于高维函数，其二阶导数就变为了一个海森矩阵，那么迭代公式就变为了

![70 10][]

我们可以看到，当Hk为正定(它的逆也为正定）的时候，可以保证牛顿法的搜索方向是向下搜索的。

牛顿法求最优值的步骤如下：  
1. 随机选取起始点x0

2. 计算目标函数f(x)在该点xk的一阶导数和海森矩阵

3. 依据迭代公式更新x

![70 11][]

如果 ![70 12][]则收敛返回，否则继续步骤2,3直至收敛 ,

我们可以看到，当我们的特征特别多的时候，求海森矩阵的逆的运算量是非常大且慢的，这对于在实际应用中是不可忍受的，因此我们想能否用一个矩阵来代替海森矩阵的逆呢，这就是拟牛顿法的基本思路。

E(f(xk+1)−f(xk))<ϵ

**二.拟牛顿法**

因为我们要选择一个矩阵来代替海森矩阵的逆，那么我们首先要研究一下海森矩阵需要具有什么样的特征才能保证牛顿法成功的应用。通过上面的描述我们知道

![70 13][]

上式我们称之为拟牛顿条件。

![70 14][]

下文图片转发自:https://blog.csdn.net/itplus/article/details/21896619

![70 15][]

![70 16][]

![70 17][]

![70 18][]

与之相比较的是最速下降算法,最速下降算法并没有使用到海森矩阵即f(x)的二阶导矩阵,他仅仅使用了函数f(x)的梯度向量,

首先给定需要优化的目标函数,梯度,海森矩阵在python中对应函数表达式fun，梯度向量gfun，和海森矩阵hess

![70 19][]

在此之外,我们需要给定一个步长.

线搜索技术是求解许多优化问题下降算法的基本组成部分，但精确线性搜索需要计算很多的

函数值和梯度值，从而耗费较多资源。特别是当迭代点远离最优点时，精确线搜索技术通常不是有效和合理的。对于许多优

化算法，其收敛速度并不依赖于精确搜索过程。因此既能保证目标函数具有可接受的下降量，又能使最终形成的迭代序列

收敛的非精确先搜索越来越流行。

![70 20][]

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

![70 21][]

1.最速下降算法的python实现

![70 22][]

该部分为公共部分,放于各个算法之前, 真实x\*=\[1,1\],目标函数值为0

import numpy as np
    
    
    # 函数表达式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]])

######################## 最速下降算法#######################
    
    #算法描述:
    
    #step 1 :选取初始点x0∈Rn,容许误差0≤ϵ≪1 .令 k←1. 
    #step 2 :计算gk=∇f(xk). 若||gk||≤ϵ，停止迭代，输出xk作为近似最优解。 
    #step 3 :取方向dk=−gk. 
    #step 4 :由线搜索技术确定步长因子αk. 
    #step 5 :令xk+1←xk+αkdk，转step 2.
    
    def gradient(fun,gfun,x0):
        #最速下降法求解无约束问题
        # x0是初始点，fun和gfun分别是目标函数和梯度
        maxk = 10000#最大迭代次数
        rho = 0.5#非精确线性搜索中的B因子
        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#用于存放非精确搜索后得到的参数mk
            while m< 20:#np.dot(a,b)a,b参数均为一维那么结果为内积运算,为多维则矩阵点乘运算
                if fun(x0+rho**m*dk) < fun(x0) + sigma*rho**m*np.dot(gk,dk):
                    mk = m#mk返回的是最小的m,m为非负整数
                    break
                m += 1
            x0 += rho**mk*dk#将步长因子代入,更新x的值
            k += 1#迭代次数加一,进行下一次迭代循环
    
        return x0,fun(x0),k
    
    print(gradient(fun,gfun,[20,10]))
    #输出结果如下,真实x*=[1,1],目标函数值为

#(array([1.00001095, 1.00002195]), 1.201869655811236e-10, 8872)

2. 阻尼牛顿法采用了牛顿法的思想，与牛顿法唯一的差别是阻尼牛顿法并不直接采用 xk\+1=xk\+dk

，而是用线搜索技术获得合适的步长因 子a用公式 x k \+ 1 = x k \+a\* d k计 算 x k \+ 1

阻尼牛顿法的算法描述：

![70 23][]

######################## 阻尼牛顿算法#######################
    
    #算法描述:
    
    #step 1: 给定终止误差0≤ϵ≪1,δ∈(0,1),σ∈(0,0.5) 初始点x0∈Rn. 令k←0
    #step 2: 计算gk=∇f(xk). 若||gk||≤ϵ，停止迭代，输出x∗≈xk 
    #step 3: 计算Gk=∇2f(xk)，并求解线性方程组得到解dk,
    #                        Gk d=−gk
    #step 4: 记mk是满足下列不等式的最小非负整数m.
    #f(xk+β^mdk)≤f(xk)+δ*β^m*(gk)T*dk,为方便起见*代表点乘
    
    #step 5: 令αk=δ^mk,xk+1=xk+αkdk,k←k+1，转step 2
    
    '''
    其中的np.linalg.solve()为求解函数
    如果x = np.linalg.solve(a, b)
    那么就有ax=b
    可以得到x=a^-1 b
    Solve the system of equations
    ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``:   
    >>> a = np.array([[3,1], [1,2]])
    >>> b = np.array([9,8])
    >>> x = np.linalg.solve(a, b)
    >>> x
    array([ 2.,  3.])  
    Check that the solution is correct:   
    >>> np.allclose(np.dot(a, x), b)
    True
    '''
    
    def dampnm(fun,gfun,hess,x0):
        # 用阻尼牛顿法求解无约束问题
        #x0是初始点，fun，gfun和hess分别是目标函数值，梯度，海森矩阵的函数
        maxk = 500#最大迭代次数
        rho = 0.55#非线性搜索中的B因子
        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)#相当于dk=-1.0*(Gk^-1)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#更新x的值
            k += 1#进入下一次循环
    
        return x0,fun(x0),k
    print(dampnm(fun,gfun,hess,[20,10]))
    #输出结果如下,真实结果x*=[1,1],目标函数值为0

#(array([1.00000067, 1.00000132]), 4.610695363186738e-13, 73)

[70]: /images/20220527/697bd75452ce4a04a2ba3798c8200d90.png
[70 1]: /images/20220527/83ca303d052a4dfaa8706e8a62f35407.png
[70 2]: /images/20220527/6a84b47b5fdd48a99d42d760a16251b5.png
[70 3]: /images/20220527/6d6480d8d2fb4992bee37e6c9b8883cd.png
[70 4]: /images/20220527/9995f09f4fa442429b146c9fdef4c7f6.png
[70 5]: /images/20220527/ed0a2c5dea944266a7bac6b63e132537.png
[70 6]: /images/20220527/c7d4ff4557ea436c9e6aae2beb0a3f37.png
[70 7]: /images/20220527/25c18c0409f649a1850d8a3680c310db.png
[70 8]: /images/20220527/89ca53d780cf475c9a94e606f73ea1bc.png
[70 9]: /images/20220527/9f7634f03b6f4fa8a58fd953211546db.png
[70 10]: /images/20220527/e92276d4510840c893eccd933031c3c7.png
[70 11]: /images/20220527/de3fcb9e4e7e4e21877999cfe107b4b7.png
[70 12]: /images/20220527/f4836822d6214208b73e3993ec8ddf3c.png
[70 13]: /images/20220527/cb5c4faf52b34a5fa7aa46fa85b4ace1.png
[70 14]: /images/20220527/079f5fe81d424810a5b9ced136ac0caa.png
[70 15]: /images/20220527/6d62df69444248fc9df8871e41f70460.png
[70 16]: /images/20220527/2d6590b5697b4870a3f4c32d8eb7d732.png
[70 17]: /images/20220527/cd0bafeb489b4702a84f24460f4c6797.png
[70 18]: /images/20220527/bb81e02c09d24956b09d7b06bc5f46df.png
[70 19]: /images/20220527/48c0f86ece084cdba5fccdc64b1be675.png
[70 20]: /images/20220527/707398a56c324499a490f962b1ae44ce.png
[70 21]: /images/20220527/5090f5fe7bb64b789962afa1cdc41f84.png
[70 22]: /images/20220527/b31a17f0a84340ce8533f987a36c0344.png
[70 23]: /images/20220527/2cffcf0e5e2c49b9a99795b72f5730b7.png