【周末AI课堂】理解梯度下降(一)（代码篇）| 机器学习你会遇到的“坑”

唐僧不用海飞丝读芯术 2019-12-26

AI课堂开讲，就差你了！

很多人说，看了再多的文章，可是没有人手把手地教授，还是很难真正地入门AI。为了将AI知识体系以最简单的方式呈现给你，从这个星期开始，芯君邀请AI专业人士开设“周末学习课堂”——每周就AI学习中的一个重点问题进行深度分析，课程会分为理论篇和代码篇，理论与实操，一个都不能少！

来，退出让你废寝忘食的游戏页面，取消只有胡吃海塞的周末聚会吧。未来你与同龄人的差异，也许就从每周末的这堂AI课开启了！

读芯术读者交流群，请加小编微信号：zhizhizhuji。等你。后台回复“周末AI课堂”，查阅相关源代码。

全文共2326字，预计学习时长5分钟

先祝大家中秋快乐~

实际的模型在训练过程中，尤其是在深度学习中，参数会达到几千万，参数空间会变的非常庞大。为了更好帮助大家理解优化算法，我们先用一维的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 定义为一个余弦函数：

我们将参数初始值设置为-15，因为当-20的时候，刚好在极大值点，梯度为零，无法进行迭代，同时将帧数设置为10，节约我们看图的时间：

def f(x):
  return(-np.cos(np.pi*x/20))
def df(x):
  return(np.sin(np.pi*x/20)*np.pi/20)
def ddf(x):
  return((np.pi/20)**2*np.cos(np.pi*x/20))

如图，牛顿法反而执行的是上升的步骤，因为Hessian此时是负定的，对于一维情况来说，二阶导数就是负值，从梯度的更新公式来看牛顿法，在Hessian负定的时候，学习率就是负的。

但梯度下降却不会出现这个问题，因为我们在执行梯度下降时，学习率永远都是正的。

如图，梯度下降应用在余弦函数构成的Loss function上，下降的结果也能得到保证。

如何解决牛顿法不降反升的问题呢？我们可以模仿Ridge regression添加L2正则化的办法，使得样本矩阵强行满秩，同样的，我们也可以在Hessian上加上一个正则化项，使得Hessian强行正定：

代码就可以很方便的写成：

def Newton(start,alpha):

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+alpha)
  return(Newton_x,Newton_y)