梯度下降与牛顿法（Python）

最近两天一直在复习李航老师的《统计学习方法》这本书上面的知识，看到了优化算法。推导了梯度下降与牛顿法的计算公式，并最终实现了相应python代码。想着记录下来，所以我就用这篇文章来记录一些要点。

梯度下降

关于梯度下降没什么好说的了，就主要是利用函数的一阶导。代码如下：

def gradient_descent_ld(grad, cur_x=0.1, learning_rate=0.01, precision=0.0001, max_iters=1000):
    """
    一维问题的梯度下降法
    :param grad: 目标函数的梯度
    :param cur_x: 当前x值，通过参数可以提供初始值
    :param learning_rate: 学习率，也相当于设置的步长
    :param precision: 设置收敛精度
    :param max_iters: 最大迭代次数
    :return: 局部最小值x
    """
    for i in range(max_iters):
        grad_cur = grad(cur_x)
        # 当梯度趋近于0时，视为收敛
        if abs(grad_cur) < precision:
            break
        cur_x = cur_x - grad_cur * learning_rate
        print("第%d次迭代的x值为：%f" % (i, cur_x))
    print("局部最小值 x=", cur_x)
    return cur_x

上面这段代码里的learning_rate步长是固定的，在[0,1]直接。但是根据定义来说，这应该是随着迭代的次数增加而变小的。这个应该不难理解，所以我在下面又实现了变步长的方法。

def gradient_descent_ld_decay(grad, cur_x=0.1, learning_rate=0.01, precision=0.0001, max_iters=1000, decay=0.5):
    """
    一维问题的梯度下降法，变步长
    :param grad: 目标函数的梯度
    :param cur_x: 当前x值，通过参数可以提供初始值
    :param learning_rate: 学习率，也相当于设置的步长
    :param precision: 设置收敛精度
    :param max_iters: 最大迭代次数
    :param decay: 学习率衰减因子
    :return:
    """
    for i in range(max_iters):
        # 新的步长
        learning_rate = learning_rate * 1.0 / (1.0 + decay * i)
        grad_cur = grad(cur_x)
        # 当梯度趋近于0时，视为收敛
        if abs(grad_cur) < precision:
            break
        cur_x = cur_x - grad_cur * learning_rate
        print("第%d次迭代的x值为：%f" % (i, cur_x))
    print("局部最小值 x=", cur_x)
    return cur_x

其实步长变化公式就是这一步：

牛顿法

牛顿法其实就是求函数二阶导，通过二阶导的正负性来进一步判断极值的存在。

 

所以，这时候要判断函数的二阶导数。如果函数是多元函数的话，这时候就要引入一个新的概念了，叫海森矩阵（Hessian）.具体怎么求，我这儿就不多说了，我们只需要记住一点。如果海森矩阵的行列式的值大于0，说明该函数有极值。我这儿要解释的一点就是关于向量对向量的求导了。具体是来自于这两步公式：

等式两边对x向量求导，可得：

这个当时有点让我迷糊了，其实可以转化为以下不等式：

其中，假设A为二阶矩阵，x为二维列向量。

其实，暴力解是可以得到答案的，比如下图所示：

但总觉得这样有失美观，所以找了一位数学系的同学，帮我重新整理了思路。她首先是把泰勒展开那两步公式给换了，如下图所示：

这才恍然大悟，代码如下：

def newton(f, x, iters):
    """
    实现牛顿法
    :param f: 原函数
    :param x: 初始值
    :param iters: 遍历的最大epoch
    :return:
    """
    Hessian_T = np.linalg.inv(hessian(f, x))
    H_G = np.matmul(Hessian_T, jacobian(f, x))
    x_new = x - H_G
    print("第1次迭代后的结果为:", x_new)
    for i in range(1, iters):
        Hessian_T = np.linalg.inv(hessian(f, x_new))
        H_G = np.matmul(Hessian_T, jacobian(f, x_new))
        x_new = x_new - H_G
        print("第"+str(i+1)+"次迭代后的结果为:", x_new)
    return x_new

 完整版代码此处可以下载。