「Machine Learning」梯度下降

一、什么是梯度下降？

1.1 口语化描述

一个风和日丽的周末，你成功登顶了泰山之巅，然而此时的喜悦还未尽兴。你却突然感觉肚子一阵隐痛，大事不妙💩。然后，坏消息是最近的厕所也在山下。

这个时候你会绞尽脑汁🤔去想，我如何才能找到一条最快的下山路径呢？

梯度下降——就是你此时此刻的救星。当你想要尽快下山，在不考虑你人体功能限制的情况下，肯定是沿着最陡峭的路下山最快。然而事实情况是，山上怪石嶙峋，并不是沿着一条直线就可以直接下山。你需要每走一步就判断下，以当前位置来看，哪条路最陡峭（下山最快）。

于是乎，就这样你决定走一步算一步，也就是每走到一个位置时，就站在这里稍微停顿一会儿，看接下来最陡峭的一条路该怎么走，……

就这样，你以最快的速度下了山，解决了燃眉之急😊。

这个每走一步就不断寻找最陡峭路径下山的过程，就是梯度下降。

下面来看看梯度下降更加正规的表达吧~

1.2 介绍

梯度下降法（Gradient Descent）并不是一个机器学习的算法，而是一种基于搜索的最优化方法。其作用就是可以用来最小化一个损失函数。与之相对应的还有一个叫做梯度上升法，其作用是用来最大化一个效用函数。

梯度下降其基本思想在于不断地逼近最优点，每一步的优化方向就是梯度的负方向。相反，梯度上升法中，进行优化的方向应该为梯度的方向。

梯度的本意是一个向量（矢量），表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。

在单变量的实值函数中，梯度可简单理解为只是导数，或者说对于一个线性函数而言，梯度就是曲线在某点的斜率。对于函数的某个特定点，它的梯度就表示从该点出发，函数值增长最为迅猛的方向（direction of greatest increase of a function）。

在一个开口向下的一元二次函数图像中，最低点左侧导数小于0，“下山”方向（$\theta - \eta \frac{d J}{d \theta}$）应该是沿$\theta$轴向右➡️（即$\theta$增大方向）；而在最低点右侧导数大于0，“下山”方向应该是沿$\theta$轴向左⬅️（即$\theta$减小方向）。

[photos]

[/photos]

求解当前位置的梯度，沿着梯度的负方向，也就是当前最陡峭的位置向下走，这样一直走下去。如果你每走一步，就计算一下当前位置的梯度（即当前这个位置最陡峭的方向），那么你所走过的路径将是下山最快的一条。

1.3 学习率

$-\eta \frac{d J}{d \theta}$中关于参数$\eta$的一些概念：

$\eta$称为学习率（Learning Rate）；
$\eta$的取值影响获得最优解的速度；
$\eta$取值不合适，甚至得不到最优解；
$\eta$是梯度下降法的一个超参数。

显而易见，参数$\eta$就是用来调节梯度（导数的），作为梯度下降的重要超参数，其有着重要的作用。可能遇到$\eta$太小或者太大的问题。

上面的四个图像，我们分别设置学习率为eta_list = [0.01, 0.1, 0.8, 1.1]，可以发现，当eta=0.01时，有点过小，导致“下山”过程中，“步长”太短，需要走很多步才可以下山；当eta=0.1时，“步长”非常合适；当eta=0.8时，“步长”太大，导致一步直接又“上山”了，不过幸好还是又下来了；当eta=1.8时，“步长”巨大，导致“下山”过程不收敛，反而越下越高了。

1.4 局部最优解与全局最优解

并不是所有函数都有唯一的极值点：

解决方案：

多次运行，随机化初始点；
梯度下降法的初始点也是一个超参数。

$$ \sum_{i=1}^{m}(y^{(i)} - \hat y^{(i)})^2 $$

线性回归法的损失函数是一个开口向下的二次函数，具有唯一的最优解。

1.5 梯度下降的数学概念

对于单变量连续函数$y=f(x)$而言，其导数定义为：

$$ f^{\prime}\left( x \right) =\lim_{\varDelta x\rightarrow 0} \frac{f\left( x+\varDelta x \right) -f\left( x \right)}{\varDelta x} $$

这个导数的定义也有另一层含义，如果$x$作微小的$\varDelta x$变化，那么$y$会如何变化呢？

$$ y=f\left( x+\varDelta x \right) \fallingdotseq f\left( x \right) +\varDelta x\cdot f^{\prime}\left( x \right) $$

将上式推广到三位空间就是：

$$ z\left( x, y \right) =f\left( x+\varDelta x, y+\varDelta y \right) \fallingdotseq f\left( x, y \right) +\varDelta x\cdot \frac{\partial f\left( x, y \right)}{\partial x}+\varDelta y\cdot \frac{\partial f\left( x, y \right)}{\partial y} $$

所以说，此时在三维空间中的微小变化$\varDelta z$就可以表示为：

$$ \begin{equation*} \begin{split} \varDelta z&=f\left( x+\varDelta x, y+\varDelta y \right) -f\left( x, y \right) \\ &\fallingdotseq \varDelta x\cdot \frac{\partial f\left( x, y \right)}{\partial x}+\varDelta y\cdot \frac{\partial f\left( x, y \right)}{\partial y} \\ &=\left[ \varDelta x, \varDelta y \right] \cdot \left[ \frac{\partial f\left( x, y \right)}{\partial x}, \frac{\partial f\left( x, y \right)}{\partial y} \right] \end{split} \end{equation*} $$

经过这样的等式变化，可以发现$\varDelta z$其实就是$\left[ \varDelta x, \varDelta y \right]$与$\left[ \frac{\partial f\left( x, y \right)}{\partial x}, \frac{\partial f\left( x, y \right)}{\partial y} \right]$两个向量的内积。

而$\left[ \frac{\partial f\left( x, y \right)}{\partial x}, \frac{\partial f\left( x, y \right)}{\partial y} \right]$被称为函数$f(x, y)$在点$(x, y)$处的梯度。

由向量内积的定义可以知道，对于两个向量$a$、$b$：

$$ a\cdot b=\left| a \right|\left| b \right|\cos \theta $$

当两者的夹角$\theta=0\degree$时，内积取得最大值（即$a$、$b$同向）；两者的夹角$\theta=180\degree$时，内积取得最小值（即$a$、$b$反向）。

所以说，如果我们的上山路径$\left[ \varDelta x, \varDelta y \right]$如果与梯度$\left[ \frac{\partial f\left( x, y \right)}{\partial x}, \frac{\partial f\left( x, y \right)}{\partial y} \right]$的方向一致，我们走的距离$\varDelta z$才会最大；如果我们下山，需要保证$\left[ \varDelta x, \varDelta y \right]$与梯度的负方向$-\left[ \frac{\partial f\left( x, y \right)}{\partial x}, \frac{\partial f\left( x, y \right)}{\partial y} \right]$一致，才可以最短距离下山。

1.6 实例：使用Excel计算梯度下降

以$z=x^2+y^2$为例，该函数是一个开口向上的抛物线曲面，最低点显然在$(x,y)=(0,0)$处，设学习率$\eta =0.1$、起始点$(x,y)=(3,2)$（这个点是随机定的）开始进行计算：

$i$	$x_i$	$y_i$	$\frac{\partial z}{\partial x}=2x$	$\frac{\partial z}{\partial y}=2y$	$\varDelta x$	$\varDelta y$	$z$	$\varDelta x\frac{\partial z}{\partial x}+\varDelta y\frac{\partial z}{\partial y}$	$f\left( x+\varDelta x, y+\varDelta y \right) -f\left( x \right) $
0	3.00	2.00	6.00	4.00	-0.60	-0.40	13.00	-5.20	-4.68	-0.52
1	2.40	1.60	4.80	3.20	-0.48	-0.32	8.32	-3.33	-3.00	-0.33
2	1.92	1.28	3.84	2.56	-0.38	-0.26	5.32	-2.13	-1.92	-0.21
3	1.54	1.02	3.07	2.05	-0.31	-0.20	3.41	-1.36	-1.23	-0.14
4	1.23	0.82	2.46	1.64	-0.25	-0.16	2.18	-0.87	-0.79	-0.09
5	0.98	0.66	1.97	1.31	-0.20	-0.13	1.40	-0.56	-0.50	-0.06
6	0.79	0.52	1.57	1.05	-0.16	-0.10	0.89	-0.36	-0.32	-0.04
7	0.63	0.42	1.26	0.84	-0.13	-0.08	0.57	-0.23	-0.21	-0.02
8	0.50	0.34	1.01	0.67	-0.10	-0.07	0.37	-0.15	-0.13	-0.01
9	0.40	0.27	0.81	0.54	-0.08	-0.05	0.23	-0.09	-0.08	-0.01
10	0.32	0.21	0.64	0.43	-0.06	-0.04	0.15	-0.06	-0.05	-0.01
11	0.26	0.17	0.52	0.34	-0.05	-0.03	0.10	-0.04	-0.03	0.00
12	0.21	0.14	0.41	0.27	-0.04	-0.03	0.06	-0.02	-0.02	0.00
13	0.16	0.11	0.33	0.22	-0.03	-0.02	0.04	-0.02	-0.01	0.00
14	0.13	0.09	0.26	0.18	-0.03	-0.02	0.03	-0.01	-0.01	0.00
15	0.11	0.07	0.21	0.14	-0.02	-0.01	0.02	-0.01	-0.01	0.00
16	0.08	0.06	0.17	0.11	-0.02	-0.01	0.01	0.00	0.00	0.00
17	0.07	0.05	0.14	0.09	-0.01	-0.01	0.01	0.00	0.00	0.00
18	0.05	0.04	0.11	0.07	-0.01	-0.01	0.00	0.00	0.00	0.00
19	0.04	0.03	0.09	0.06	-0.01	-0.01	0.00	0.00	0.00	0.00
20	0.03	0.02	0.07	0.05	-0.01	0.00	0.00	0.00	0.00	0.00
21	0.03	0.02	0.06	0.04	-0.01	0.00	0.00	0.00	0.00	0.00
22	0.02	0.01	0.04	0.03	0.00	0.00	0.00	0.00	0.00	0.00
23	0.02	0.01	0.04	0.02	0.00	0.00	0.00	0.00	0.00	0.00
24	0.01	0.01	0.03	0.02	0.00	0.00	0.00	0.00	0.00	0.00
25	0.01	0.01	0.02	0.02	0.00	0.00	0.00	0.00	0.00	0.00
26	0.01	0.01	0.02	0.01	0.00	0.00	0.00	0.00	0.00	0.00
27	0.01	0.00	0.01	0.01	0.00	0.00	0.00	0.00	0.00	0.00
28	0.01	0.00	0.01	0.01	0.00	0.00	0.00	0.00	0.00	0.00
29	0.00	0.00	0.01	0.01	0.00	0.00	0.00	0.00	0.00	0.00
30	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	0.00

以第一次计算为例：$(x,y)=(3,2)$，此时梯度为：

$$ \left[ \frac{\partial z}{\partial x}, \frac{\partial z}{\partial y} \right] =\left[ 2x, 2y \right] =\left[ 2\times 3, 2\times 2 \right] =\left[ 6, 4 \right] $$

所以说此时的“下山步长”为：

$$ -\eta \left[ \frac{\partial z}{\partial x}, \frac{\partial z}{\partial y} \right] =\left[ -0.6, -0.4 \right] $$

“下一步落脚点”的坐标为：

$$ \left[ x_i, y_i \right] -\eta \left[ \frac{\partial z}{\partial x}, \frac{\partial z}{\partial y} \right] =\left[ 3, 2 \right] -0.1\times \left[ 6, 4 \right] =\left[ 2.4, 1.6 \right] $$

由计算过程我们也可以看出：

$$ \varDelta x\frac{\partial z}{\partial x}+\varDelta y\frac{\partial z}{\partial y} \fallingdotseq f\left( x+\varDelta x, y+\varDelta y \right) -f\left( x \right) $$

两者之间并不是完全相等的，但是随着“逐渐靠近山下”，两者之间的差（最后一列）会越来越小，直到相等。

二、实现一个最简单的梯度下降示例

这里我们使用一个一元二次函数：

$$ f(x)=(x-2.5)^2-1 $$

来可视化梯度下降求其最小值的过程。

定义函数：

from sympy.abc import x
from sympy import lambdify, diff

# Define function
f_x = (x - 2.5) ** 2 - 1
# Calculate f(x)
f_x_fcn = lambdify(x, f_x)

# Calculate f'(x)
f_x_1_diff = diff(f_x, x)
f_x_1_diff_fcn = lambdify(x, f_x_1_diff)

x_arr = np.linspace(-1, 6, 200)
y_arr = f_x_fcn(x_arr)

这里我们使用了sympy.lambdify()创建了函数$f(x)$，然后使用sympy.diff()求出其关于x的导数。

由二元一次函数的性质可以求出，$f(x)$的最低点应该是出现在$x=2.5$处，而$f(2.5)=-1$ 。计算最值点：

lowest_point_x = 2.5
lowest_point_y = f_x_fcn(lowest_point_x)
(lowest_point_x, lowest_point_y)

绘图观察其图像：

# Plot image and data visualization
plot_x_min, plot_x_max, plot_y_min, plot_y_max = x_arr.min() - 1, x_arr.max() + 1, y_arr.min() - 1, y_arr.max() + 1
plt.plot(x_arr, y_arr)
plt.hlines(lowest_point_y, plot_x_min, lowest_point_x, 'r', '--')
plt.vlines(lowest_point_x, plot_y_min, lowest_point_y, 'r', '--')
plt.xlim(plot_x_min, plot_x_max)
plt.ylim(plot_y_min, plot_y_max)
plt.text(lowest_point_x, lowest_point_y+2, f'$x={lowest_point_x},f({lowest_point_x})={lowest_point_y}$', horizontalalignment ='center')
plt.show()

梯度下降求最小值：

参数	含义
`learn_rate`	学习率，调节步长
`_theta`	初始位置，开始“下山”的地方
`epsilon`	迭代终止条件：当上一步与这一步相比，走的距离小于这个值时，就说明已经到达山底了
`theta_history`	存放“下山”过程每一步的位置
`i`	计数器，看看走几步才可以下山

learn_rate, _theta, epsilon = 0.1, 0, 1e-8  
theta_history = []
i = 0  # Counter
while True:
    i += 1
    gradient = f_x_diff_fcn(_theta)
    last_theta = _theta
    theta_history.append(last_theta)
    _theta = _theta - learn_rate * gradient

    if abs(f_x_fcn(last_theta) - f_x_fcn(_theta)) < epsilon:
        hues.success(f'theta is {_theta}, f(x) = {f_x_fcn(_theta)}, the counter is {i}.')
        break

# Plot image and data visualization
plt.plot(x_arr, y_arr)
plt.scatter(theta_history, f_x_fcn(np.array(theta_history)), s=100, color='r', marker='+')
plt.show()

10:36:09 - SUCCESS - theta is 2.499891109642585, f(x) = -0.99999998814289, the counter is 45.

由上面的迭代过程可以看出，当设置学习率learn_rate=0.1、初始点为_theta=0、迭代终止条件epsilon=1e-8时，经过45次迭代，找到了最小值点2.499891109642585，其对应的函数值为-0.99999998814289。

当我们分别设置学习率为eta_list = [0.01, 0.1, 0.8, 1.1]时，对应的迭代次数分别为：

学习率	迭代次数
0.01	423
0.1	45
0.8	21
1.1	OverflowError

由此可见，设置一个正确的学习率，对于梯度下降过程的执行效率有着至关重要的作用。

三、二元函数的梯度下降

看懂了上面简单的一元函数梯度下降，再来一个复杂点的二元函数：

$$ z=f(x, y) =x - y + 2 x^2 + 2xy + y^2 $$

来看看。与一元函数相比，二元函数使用梯度下降求解极小值的过程才更像“走最陡峭的路下山”过程。

首先来定义下函数，并创造出数据：

from sympy import lambdify, diff
from sympy.abc import x, y
import numpy as np
from matplotlib import pyplot as plt
import hues

num = 400;  # number of mesh grids
x_array = np.linspace(-4, 4, num)
y_array = np.linspace(-4, 4, num)
xx, yy = np.meshgrid(x_array, y_array)

# 定义函数
f_xy = x - y + 2 * x * x + 2 * x * y + y * y
f_xy_fcn = lambdify([x, y], f_xy)
# 计算网格
f_xy_zz = f_xy_fcn(xx, yy)

然后绘制出图像来看下，看看曲面是什么样子的：

plt.figure()
ax = plt.axes(projection="3d")
ax.patch.set_facecolor("white")  #设置 axes 背景颜色
ax.plot_surface(xx, yy, f_xy_zz, alpha=0.9, cmap=plt.cm.jet)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("$z=f(x, y) =x - y + 2 x^2 + 2xy + y^2$")
ax.view_init(elev=30, azim=70)
plt.show()

由图中可以看出来，$z=f(x, y)$近似在$(x, y)=(4, 4)$以及$(x, y)=(-4, -4)$位置取得极大值（“山”的最高点），而在$(x, y)=(4, -4)$到$(x, y)=(-4, 4)$的“山谷”中取得极小值，即这里是山底。既然这个曲线肯定是有极小值的。也就是说，使用梯度下降的方法一定是可以找到最优解的。

接下来就可以愉快地进行梯度下降吧~

首先求出$z=f(x, y)$关于$x$和$y$的偏导数，对于一个一元函数，梯度就是其导数，而对于一个二元或者多元函数，梯度就是其各个自变量的偏导数。

$$ gradf(x, y)=\nabla f(x,y) = \{\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\} $$

# partial derivative with respect to x
df_dx = f_xy.diff(x)
df_dx_fcn = lambdify([x, y], df_dx)

# partial derivative with respect to y
df_dy = f_xy.diff(y)
df_dy_fcn = lambdify([x, y], df_dy)

定义了偏导数之后，就可以开始下山过程——梯度下降了🥰：

#梯度下降
learn_rate, epsilon = 0.05, 7e-9  # 学习率与迭代终止条件
start_x, start_y = -4, -4  # 开始迭代点
descent_point = [(start_x, start_y, f_xy_fcn(start_x, start_y))]  # 记录下山的每一步

i = 1 # Counter
while True:
    x, y = descent_point[-1][0], descent_point[-1][1]  # 取出“这一步”的x，y
    new_x = x - learn_rate * df_dx_fcn(x, y)  # 根据偏导数计算“下一步”
    new_y = y - learn_rate * df_dy_fcn(x, y)
    
    # 记录“下一步”的数据
    descent_point.append((new_x, new_y, f_xy_fcn(new_x, new_y)))

    if f_xy_fcn(x, y) - f_xy_fcn(new_x, new_y) < epsilon:  # 迭代终止条件
        hues.success(f'The counter is {i}, the last descent_point is {descent_point[-1]}.')
        break

    i += 1

输出：13:38:48 - SUCCESS - The counter is 226, the last descent_point is (-0.9997546279984584, 1.4996029797616182, -1.2499999167953932).
说明下山一共走了226步，山的最低点为(-0.9997546279984584, 1.4996029797616182, -1.2499999167953932)。

取出下山过程中，每一步所处位置的x、y、z坐标值：

descent_point_x = [i[0] for i in descent_point]
descent_point_y = [i[1] for i in descent_point]
descent_point_z = [i[2] for i in descent_point]

绘制出下山每一步的点：

plt.figure()
ax = plt.axes(projection="3d")
ax.patch.set_facecolor("white")  #设置 axes 背景颜色
ax.plot_surface(xx, yy, f_xy_zz, alpha=0.3, cmap=plt.cm.jet)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.view_init(elev=30, azim=70)
ax.plot(descent_point_x, descent_point_y, descent_point_z, 'r.')
plt.show()

也可以在二维平面上观察下：

# 创建画布
fig, ax = plt.subplots()

# 绘制函数f(x, y)的热图
colorbar = ax.contourf(xx, yy, f_xy_zz, 20, cmap='RdYlBu_r')
fig.colorbar(colorbar, fraction=0.046, pad=0.17, label=f'$z=f(x, y) =x - y + 2 x^2 + 2xy + y^2$',
             orientation='horizontal')

# 绘制出函数f(x, y)的最低点+0.1的区域
ax.contour(xx, yy, f_xy_zz, levels=[np.min(f_xy_zz) + 0.1],
           colors='red',
           linestyles='-')

# 设置x、y轴限制
ax.set_xlim(xx.min(), xx.max())
ax.set_ylim(yy.min(), yy.max())
# 轴名称
ax.set_xlabel('x')
ax.set_ylabel('y')
# x、y轴等间距
plt.gca().set_aspect('equal', adjustable='box')
# 绘制图形
ax.plot(descent_point_x, descent_point_y, descent_point_z, 'r.')
plt.show()

四、使用梯度下降求解线性回归问题

首先来创建一组数据：

import numpy as np
from matplotlib import pyplot as plt

np.random.seed(666)
x = 2 * np.random.random(size=100)
X = x.reshape(-1, 1)

y = x * 3 + 4 + np.random.normal(size=100)

X.shape, y.shape

X为一个100×1的特征矩阵，y为标签值，对应有100个。且两者之间“大概”存在为：

$$ y=3x+4 $$

的关系。

绘制出图形观察下：

plt.scatter(X, y)

线性回归问题，本质上是求损失函数：

$$ \begin{equation*} \begin{split} J &= \frac{1}{2m}\sum_{i=1}^{m} (y^{(i)} - \hat y ^{(i)})^2 \\ &= \frac{1}{2m}\sum_{i=1}^{m}\left(y^{(i)}-\theta_{0}-\theta_{1} X_{1}^{(i)}-\theta_{2} X_{2}^{(i)}-\ldots-\theta_{n} X_{n}^{(i)}\right)^{2} \end{split} \end{equation*} $$

令$X_0$为全为1的列向量：

特征矩阵$X$为$m×n$，则$y$为$m×1$，$X_0$为$m×1$。将$X_0$放置到$X$最左侧“组装”成$X_b$为$m×(n+1)$。$\theta$为$n×1$，将$\theta_{0}$补充到$\theta$的第一位，那么新的$\theta$为$(n+1)×1$。

$$ \begin{equation*} \begin{split} J &= \frac{1}{2m}\sum_{i=1}^{m}\left(y^{(i)}-\theta_{0}X_0^{(i)}-\theta_{1} X_{1}^{(i)}-\theta_{2} X_{2}^{(i)}-\ldots-\theta_{n} X_{n}^{(i)}\right)^{2} \\ &= \frac{1}{2m}\sum_{i=1}^{m}\left(y^{(i)}-\theta_{0}X_0^{(i)}-\theta_{1} X_{1}^{(i)}-\theta_{2} X_{2}^{(i)}-\ldots-\theta_{n} X_{n}^{(i)}\right)^{2} \\ &= \frac{1}{2m}(y-X_b\theta)^2 \end{split} \end{equation*} $$

达到最小值时所对应的最优解。这个时候我们就可以来定义损失函数了：

def loss_fun(theta, X_b, y):
    try:
        return np.sum((y - X_b.dot(theta) ** 2)) / len(X_b)
    except:
        hues.error('The return value is to large!!!')
        return np.inf

对于一个多元函数，这个时候的梯度变成了：

$$ \nabla J = \begin{bmatrix} \frac{\partial J}{\partial \theta_0}\\ \frac{\partial J}{\partial \theta_1}\\ \dots\\ \frac{\partial J}{\partial \theta_n} \end{bmatrix} = \frac{2}{m} \begin{bmatrix} \sum_{i=1}^{m} (X_b^{(i)}\theta - \hat y ^{(i)})\\ \sum_{i=1}^{m} (X_b^{(i)}\theta - \hat y ^{(i)})X_1^{(i)}\\ \dots\\ \sum_{i=1}^{m} (X_b^{(i)}\theta - \hat y ^{(i)})X_n^{(i)} \end{bmatrix} $$

所以求导函数可以定义为：

def get_der(theta, X_b, y):
    res = np.empty(len(theta))
    res[0] = np.sum(X_b.dot(theta) - y)
    for i in range(1, len(theta)):
        res[i] = (X_b.dot(theta) - y).dot(X_b[:, i])

    return res * 2 / len(X_b)

然后就可以封装一个求梯度的函数：

theta_history = []
def gradient_descent(X_b, y, ini_theta, eta, n_iters=1e4, eps=1e-8):
    theta = ini_theta
    i = 0
    theta_history.append(ini_theta)

    while i < n_iters:
        gradient = get_der(theta, X_b, y)
        last_theta = theta
        theta = theta - eta * gradient
        theta_history.append(theta)

        if abs(loss_fun(theta, X_b, y) - loss_fun(last_theta, X_b, y)) < eps:
            break

        i += 1

    return theta

Try it!

%%time
X_b = np.hstack([np.ones([len(X), 1]), X])
initial_theta = np.zeros(X_b.shape[1])
eta = 0.001

b, a = gradient_descent(X_b, y, initial_theta, eta)  # y = 3x + 4
b, a  # (3.9946613793913204, 3.029663997668306)

对于上面的$\nabla J$，我们可以进一步简化：

$$ \begin{equation*} \begin{split} \nabla J &= \frac{2}{m} \begin{bmatrix} X_b^{(1)}\theta -y^{(1)},& X_b^{(2)}\theta -y^{(2)},& \dots,& X_b^{(m)}\theta -y^{(m)} \end{bmatrix} \begin{bmatrix} X_0^{(1)}& X_1^{(1)}& \dots& X_n^{(1)} \\ X_0^{(2)}& X_1^{(2)}& \dots& X_n^{(2)} \\ \dots & \dots & \ddots & \vdots \\ X_0^{(m)}& X_1^{(m)}& \dots& X_n^{(m)} \\ \end{bmatrix} \\ &=\frac{2}{m} \cdot (X_b\theta - y)^T \cdot X_b \\ &=\frac{2}{m} \cdot X_b^T \cdot (X_b\theta - y) \end{split} \end{equation*} $$

其中$X_0$为全为1的列向量。

那么我们的求导函数就可以进一步优化：

def get_der(theta, X_b, y):
    return X_b.T.dot(X_b.dot(theta) - y) * 2 / len(X_b)

五、归一化一定要记得

使用波士顿房价数据集来验证我们上面实现的梯度下降：

加载数据：

data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep="\s+", skiprows=22, header=None)
data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])
target = raw_df.values[1::2, 2]

过滤数据：

X, y = data, target
X = X[y < 50]
y = y[y < 50]
X.shape, y.shape

为特征数据矩阵加上一列：

X_b = np.hstack([np.ones([len(X), 1]), X])

划分数据集：

# Split the Feature Dataset and Label Dataset
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test, = train_test_split(X_b, y, test_size=0.3, random_state=666)
X_train.shape, X_test.shape, y_train.shape, y_test.shape

梯度下降：

%%time
initial_theta = np.zeros(X_train.shape[1])
eta = 0.001

gradient_descent(X_train, y_train, initial_theta, eta)

这个时候我们运行代码，会发现提示Overflow溢出错误：

这是因为，当数据集中存在多列数据，且部分列数据的数量级非常大。eta如果设置的稍微一点点大，就会导致在计算部分列的导数时不收敛，进而导致溢出错误。

解决上面的问题有两种办法，第一就是我们手动设置一个极小的eta=1e-6，这样再次运行就不错报错了：

calc_theta = gradient_descent(X_train, y_train, initial_theta, eta=1e-6)
y_predict = X_test.dot(calc_theta)

from sklearn.metrics import r2_score

r2_score(y_test, y_predict) # 0.30594186189467965

虽然不再报错，但是我们发现，梯度下降求得的结果，R方只有0.306，这是完全不可用的。

接下来我们尝试增大迭代次数：

%%time
calc_theta = gradient_descent(X_train, y_train, initial_theta, eta=1e-6, n_iters=1e6)
y_predict = X_test.dot(calc_theta)

from sklearn.metrics import r2_score

r2_score(y_test, y_predict)  # 0.7440227019523431

当迭代次数上升至n_iters=1e6，R方有了显著的提高，但是程序耗时也大大增加。那么有没有什么方法可以既让耗时减小，同时又提升R方呢？

首先来分析出现该问题的根源——数据量过大，且数据表中有些列的数据量级高，进而导致Overflow溢出错误。这个问题，我们通过对数据进行归一化处理即可解决：

from sklearn.preprocessing import StandardScaler

scale_scaler = StandardScaler()
scale_scaler.fit(X_train)

X_train_std = scale_scaler.transform(X_train)
X_test_std = scale_scaler.transform(X_test)
X_train_std = np.hstack([np.ones([len(X_train), 1]), X_train_std])
X_test_std = np.hstack([np.ones([len(X_test), 1]), X_test_std])

归一化后，再次进行梯度下降：

%%time
initial_theta = np.zeros(X_train_std.shape[1])
calc_theta = gradient_descent(X_train_std, y_train, initial_theta, eta=0.001)
y_predict = X_test_std.dot(calc_theta)

from sklearn.metrics import r2_score

r2_score(y_test, y_predict)  # 0.7983573625036819

可以看到，虽然我们这次仍然设置了eta=0.001，但是却没有提示Overflow溢出错误，并且R方也得到了显著的提高，显著高于之前的0.306。说明对特征数据矩阵进行归一化处理可以解决梯度下降过程中的这个问题。

六、随机梯度下降法

上面我们实现的过程中，每一次“下山”的过程中，都会对各个方向的路径进行求导，以便找出最陡峭的下山路径，这个过程被称为批量梯度下降法（Batch Gradient Descent）。然而实际情况中，由于数据集的量级都十分大（特征多代表下山的路径多），因此我们可以采用随机梯度下降法（Stochastic Gradient Descent）来进行优化。

所谓随机梯度下降法，就是在当前位置，我们只随机挑选几条路，然后看看这几条路的陡峭程度（求导），然后从这几条中选出最陡峭的进行“下山”。

有了想法，实现起来也就很简单了：

def get_der(theta, X_b, y):
    rand_index = np.random.randint(len(X_b))
    X_b_i, y_i = X_b[rand_index], y[rand_index]
    # X_b_i, y_i
    return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2

上面的求导函数，我们每次只从一个维度（一个特征，即一条下山路线）计算其导数。

对于求梯度的函数，也需要更新一下：

theta_history = []
def gradient_descent(X_b, y, ini_theta, eta, n_iters=1e4, eps=1e-8):
    theta = ini_theta
    i = 0
    theta_history.append(ini_theta)
    t_0, t_1 = 5, 50

    def learning_rate(t):
        return t_0 / (t + t_1)

    while i < n_iters:
        gradient = get_der(theta, X_b, y)
        theta = theta - learning_rate(i) * gradient
        theta_history.append(theta)
        i += 1

    return theta

这里我们对学习率进行了动态计算，当刚开始“下山”的时候，“步子”迈大一点，而快到山底时，每次“只走一小步”。这是因为，采用随机梯度下降的时候，由于每次只计算了一个方向的梯度，可能在下降到山底的时候，所计算下一步方向的梯度非常大，即沿这个方向可能不再是下山，而是上山了，并不能保证这个方向真是最优的。

而且，采用随机梯度下降的过程因为计算量大大减小，因此我们可以只管“下山”，而不用像以前一样每“下去”一步，就要计算这一步和上一步相比，走了多远，如果距离足够小，就说明已经到了山底。

由上图中可以看出来，随机梯度下降并不像批量梯度下降，每一步的下一步都指向梯度最大（下降最快）的方向，但是其整体还是慢慢下山的。

七、在scikit-learn中使用梯度下降

from sklearn.linear_model import SGDRegressor

sgd_reg = SGDRegressor()
sgd_reg.fit(X_train_std, y_train)
sgd_reg.score(X_test_std, y_test)

注意 / WARNING
在scikit-learn中，随机梯度下降类被封装到linear_model包中，这说明其只可以用来计算线性相关的问题

sgd_reg = SGDRegressor(max_iter=100)
sgd_reg.fit(X_train_std, y_train)
sgd_reg.score(X_test_std, y_test)

当前页面是本站的「Google AMP」版。查看和发表评论请点击：完整版 »