牛顿法详解（Python实现）

除了前面说的梯度下降法，牛顿法也是机器学习中用的比较多的一种优化算法。

牛顿法的基本思想是利用迭代点 x_k 处的一阶导数（梯度）和二阶导数（Hessen矩阵）对目标函数进行二次函数近似，然后把二次模型的极小点作为新的迭代点，并不断重复这一过程，直至求得满足精度的近似极小值。

牛顿法下降的速度比梯度下降的快，而且能高度逼近最优值。牛顿法分为基本牛顿法和全局牛顿法。

基本牛顿法的原理

基本牛顿法是一种基于导数的算法，它每一步的迭代方向都是沿着当前点函数值下降的方向。对于一维的情形，对于一个需要求解的优化函数 f(x)，求函数的极值的问题可以转化为求导函数 f′(x)=0。对函数 f(x) 进行泰勒展开到二阶，得到:
对上式求导并令其为 0，则为：
即得到： 这就是牛顿法的更新公式。

基本牛顿法的流程

1. 给定终止误差值 0≤ε<<1，初始点 x₀∈Rⁿ ，令 k=0；
2. 计算 g_k=▽f（x_k），若 ||g_k||≤ε，则停止，输出 x≈x_k；
3. 计算 G_k =▽²f(x_k)，并求解线性方程组 G_kd=-g_k 得解 d_k；
4. 令 x_k+1=x_k+d_k ，k=k+1，并转 2；

全局牛顿法

牛顿法最突出的优点是收敛速度快，具有局部二阶收敛性，但是，基本牛顿法初始点需要足够“靠近”极小点，否则，有可能导致算法不收敛，此时就引入了全局牛顿法。

全局牛顿法的流程为：

• 给定终止误差值 0≤ε<<1，δ∈（0，1），σ∈（0，0.5），初始点 x₀∈Rⁿ ，令 k=0；
• 计算 g_k=▽f(x_k)，若 ||g_k||≤ε，则停止，输出 x*≈x_k；
• 计算 G_k=▽²f(x_k)，并求解线性方程组 G_kd=-g_k 得解 d_k；
• 设 m_k 是不满足下列不等式的最小非负整数 m：

• 令 α_k=δm_k ，x_k+1=x_k+α_kd_k，k=k+1，并转 2；

全局牛顿法的具体实现如下程序所示：

def newton(feature, label, iterMax, sigma, delta):
    '''牛顿法
    input:  feature(mat):特征
            label(mat):标签
            iterMax(int):最大迭代次数
            sigma(float), delta(float):牛顿法中的参数
    output: w(mat):回归系数
    '''
    n = np.shape(feature)[1]
    w = np.mat(np.zeros((n, 1)))
    it = 0
    while it <= iterMax:
        # print it
        g = first_derivativ(feature, label, w)  # 一阶导数
        G = second_derivative(feature)  # 二阶导数
        d = -G.I * g
        m = get_min_m(feature, label, sigma, delta, d, w, g)  # 得到最小的m
        w = w + pow(sigma, m) * d
        if it % 10 == 0:
            print "	---- itration: ", it, " , error: ", get_error(feature, label , w)[0, 0]
        it += 1      
    return w

程序中，函数 newton 利用全局牛顿法对线性回归模型中的参数进行学习，函数 newton 的输入为训练特征 feature、训练的目标值 label、全局牛顿法的最大迭代次数 iterMax 以及全局牛顿法的两个参数 sigma 和 delta。输出是线性回归模型的参数 w。

在函数 newton 中需要计算损失函数的一阶导数，二阶导数，同时需要计算最小的 m 值，最终根据上述的值更新权重。

求最小 m 值的函数 get_min_m 的具体实现如下所示：

def get_min_m(feature, label, sigma, delta, d, w, g):
    '''计算步长中最小的值m
    input:  feature(mat):特征
            label(mat):标签
            sigma(float),delta(float):全局牛顿法的参数
            d(mat):负的一阶导数除以二阶导数值
            g(mat):一阶导数值
    output: m(int):最小m值
    '''
    m = 0
    while True:
        w_new = w + pow(sigma, m) * d
        left = get_error(feature, label , w_new)
        right = get_error(feature, label , w) + delta * pow(sigma, m) * g.T * d
        if left <= right:
            break
        else:
            m += 1
    return m

在计算的过程中，计算损失函数值时使用到了 get_error 函数，其具体实现如下所示：

def get_error(feature, label, w):
    '''计算误差
    input:  feature(mat):特征
            label(mat):标签
            w(mat):线性回归模型的参数
    output: 损失函数值
    '''
    return (label - feature * w).T * (label - feature * w) / 2

求一阶导数的函数 first_derivativ 的具体实现如下所示：

def first_derivativ(feature, label, w):
    '''计算一阶导函数的值
    input:  feature(mat):特征
            label(mat):标签
    output: g(mat):一阶导数值
    '''
    m, n = np.shape(feature)
    g = np.mat(np.zeros((n, 1)))
    for i in xrange(m):
        err = label[i, 0] - feature[i, ] * w
        for j in xrange(n):
            g[j, ] -= err * feature[i, j]
    return g    

求二阶导数的函数 second_derivative 的具体实现如下所示：

def second_derivative(feature):
    '''计算二阶导函数的值
    input:  feature(mat):特征
    output: G(mat):二阶导数值
    '''
    m, n = np.shape(feature)
    G = np.mat(np.zeros((n, n)))
    for i in xrange(m):
        x_left = feature[i, ].T
        x_right = feature[i, ]
        G += x_left * x_right
    return G

Armijo搜索

全局牛顿法是基于 Armijo 的搜索，满足Armijo准则：给定 β∈(0，1)，σ∈(0，0.5)，令步长因子 α_k =βm_k ，其中 m_k 是满足下列不等式的最小非负整数：

利用全局牛顿法求解线性回归模型

假设有 m 个训练样本，其中，每个样本有 n-1 个特征，则线性回归模型的损失函数为：

若是利用全局牛顿法求解线性回归模型，需要计算线性回归模型损失函数的一阶导数和二阶导数，其一阶导数为：

损失函数的二阶导数为：
具体实现分别如 first_derivativ 和 second_derivative 所示。

本节可运行代码下载地址：链接: https://pan.baidu.com/s/1psgU0qBQG9xToralDwaTdQ 提取码: bycg