理解梯度下降（一）（代码篇）｜机器学习你会遇到的“坑”

先祝大家中秋快乐~

实际的模型在训练过程中，尤其是在深度学习中，参数会达到几千万，参数空间会变的非常庞大。为了更好帮助大家理解优化算法，我们先用一维的Loss function，在代码的演示过程中，我将不会展示坐标下降，不仅因为一维的坐标下降和基于梯度的更新都是在同一个方向上进行搜索，更重要的是，此代码篇希望大家通过这一节可以掌握GIF图的绘制和可视化优化算法，它会在我们的后续的课程中被广泛使用。

我们假设Loss function对于参数来说是一个非常简单的二次函数：

那么它的梯度就是：

def f(x):

return(x**2)

def df(x):

return(2*x)

同时根据上一节理论篇所说的梯度下降的公式：

我们可以定义一个用于梯度下降的函数，它接受参数的初始值和学习率，返回每次更新以后的Loss值和参数值，注意到我们此时只迭代了100次：

def GD(lr,start):

x = start

GD_x, GD_y = [], []

for it in range(100):

GD_x.append(x)

GD_y.append(f(x))

dx = df(x)

x = x - lr * dx

return(GD_x,GD_y)

在实际使用过程中，我们不需要保存每一迭代的结果，当参数非常多时，对计算机的内存消耗也会非常大，但是我们需要利用每一步的结果直观地给大家展示下降的过程。关于代码很多人都能快速的理解甚至应用，但是展示下降的过程，需要我们画出动态图，我们需要用到matplotlib中的animation，此外如果我们需要保存GIF，我们还需要安装ffmpeg。

在实现GIF的代码实操中，我们需要明确两步：

不会被更新的图形部分不断更新的图形部分和图形的初始化

对于第一步，我们可以很方便的写出：

import matplotlib.pyplot as plt

import seaborn as sns

sns.set(style='darkgrid')

FFwriter =animation.FFMpegWriter()

fig, ax = plt.subplots()

fig.set_tight_layout(True)

points_x = np.linspace(-20, 20, 1000)

points_y = f(points_x)

ax.plot(points_x,points_y, c="k", alpha=0.9, linestyle="-")

对于第二步，我们首先要获得梯度下降的结果：

GD_x,GD_y=GD(lr=pow(2,-10)*16,start=-20)

然后初始化，定义一个更新公式，将我们上述获得结果依次更新：

point_line,=ax.plot(GD_x[0],GD_y[0],'or')

def update(i):

label = 'timestep {0}'.format(i)

print(label)

point_line.set_xdata(GD_x[i])

point_line.set_ydata(GD_y[i])

ax.set_xlabel(label)

return point_line, ax

将fig，update作为参数传入到我们的animation的FuncAnimation类中，并设置帧数为60帧，每帧间隔200ms：

anim = FuncAnimation(fig, update, frames=np.arange(0, 60), interval=200)

我们的代码总结如下：

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import animation

from matplotlib.animation import FuncAnimation

import seaborn as sns

sns.set(style='darkgrid')

FFwriter =animation.FFMpegWriter()

fig, ax = plt.subplots()

fig.set_tight_layout(True)

def f(x):

return(x**2)

def df(x):

return(2*x)

points_x = np.linspace(-20, 20, 1000)

points_y = f(points_x)

ax.plot(points_x,points_y, c="k", alpha=0.9, linestyle="-")

def GD(lr,start):

x = start

GD_x, GD_y = [], []

for it in range(100):

GD_x.append(x)

GD_y.append(f(x))

dx = df(x)

x = x - lr * dx

return(GD_x,GD_y)

GD_x,GD_y=GD(lr=pow(2,-10)*16,start=-20)

print('fig size: {0} DPI, size in inches {1}'.format(

fig.get_dpi(), fig.get_size_inches()))

point_line,=ax.plot(GD_x[0],GD_y[0],'or')

def update(i):

label ='timestep {0}'.format(i)

print(label)

point_line.set_xdata(GD_x[i])

point_line.set_ydata(GD_y[i])

ax.set_xlabel(label)

return point_line, ax

if__name__=='__main__':

anim = FuncAnimation(fig, update, frames=np.arange(0, 100), interval=200)

anim.save('GD.gif', writer=FFwriter)

如图，我们实现了Loss function上实现了梯度下降法，同时，可以看出，梯度在越变越小，在固定好学习率的情形下，下降的也越来越慢，迭代60次之后仍然没有到达最低点。

我们可以尝试将学习率设置的大一点，看看是否会更快的迭代到最小值：

......

GD_x,GD_y=GD(lr=pow(2,-7)*16,start=-20)

......

如图，我们更快的到达了极小值，迭代十次就达到了，后面的迭代是不必要的。

如果继续增大学习率会发生什么情况呢？我们继续增加学习率:

......

GD_x,GD_y=GD(lr=pow(2,-4)*16,start=-20)

......

如图，当学习率继续增大的时候,参数点会在极小值附近来回震荡，可想而知，如果我们继续增大学习率，更新点会超出图形的范围。

接下来尝试牛顿法的时候，需要注意到此时的Loss function的二阶导数是一个常数，如果我们直接利用牛顿法的更新公式：

就会发现此时的牛顿法不过是另一个学习率下的梯度下降，学习率固定为二阶导数的倒数，我们要想在牛顿法和梯度下降法中看到差别，就需要更换Loss function，为达到此目的，我们原来Loss function中添加一个小的三角函数项：

三角函数好处在于它的无穷阶不仅可导，而且不会成为常数，同时也不会对原来函数的形状有太大的影响，我们可以很方便的写出：

def f(x):

return(-np.cos(np.pi*x/20)+x**2)

def df(x):

return(np.sin(np.pi*x/20)*np.pi/20+2*x)

def ddf(x):

return((np.pi/20)**2*np.cos(np.pi*x/20)+2)

这样就可以保证二阶导数不是一个常数，接下来我们构建牛顿法的更新函数：

def Newton(start):

x = start

Newton_x, Newton_y = [], []

for it in range(100):

Newton_x.append(x), Newton_y.append(f(x))

g = df(x)

h = ddf(x)

x = x -g/h

return(Newton_x,Newton_y)

然后我们利用该函数获得更新点，然后更新到GIF上：

......

newton_x,newton_y=Newton(start=-20)

point_line,=ax.plot(newton_x[0],newton_y[0],'or')

def update(i):

label ='timestep {0}'.format(i)

print(label)

point_line.set_xdata(newton_x[i])

point_line.set_ydata(newton_y[i])

ax.set_xlabel(label)

return point_line, ax

......

如图，参数只用了一次迭代就到了最低点，后面的迭代就不起作用，因为梯度为零。

牛顿法看起来很快，但我们在不清楚loss function性质的情况下却很少使用它，尤其是在深度学习中，这不仅是因为需要每一步需要计算Hessian，还因为我们如果真的希望牛顿法执行下降步骤，必须保证Hessian的正定，这在某些情况下是不成立的，比如，我们直接将Loss function 定义为一个余弦函数：