python牛顿法与拟牛顿法_梯度下降与牛顿法（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

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

牛顿法

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

比如：

,该函数的一阶导为：

.在x=0处是可以得到

.但是很明显x=0并不是该函数的极小值。如下图所示：

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

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

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

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

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

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

这才恍然大悟，这其中有一点我之前是忘了，就是以下公式：

，其中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

完整版代码

2019.11.20 修改

这两天在看拟牛顿法相关内容的时候，突然发现了牛顿法和梯度下降之间的相关联的地方。当时是越看越奇怪，后来仔细比对了公式。才发现两者之间最大的区别。不多说，先看我的手写稿子。

其实一阶导的牛顿法就是相当于直线方程求零点。梯度下降则是x值逐步往函数最小值方向移动，所以大多数情况下是牛顿法迈的步子大些，也就意味着迭代的次数所需更少。梯度下降公式不是由原有的泰勒展开公式严格推导出来的，牛顿法则是可严格推导出来。我当时就是在这两步犯了迷糊。

这张图则是很形象的描述了它们的迭代轨迹。