4. 线性回归的梯度下降求解过程

在掌握梯度下降基本原理即基本性质之后，我们可以尝试在手动创建的线性回归数据集上进行梯度下降方法的参数求解。
围绕 Lesson 3.3 中的数据集进行线性回归参数求解，梯度下降求解过程仍然按照上述过程执行：
（1）确定数据集和模型

仍然还是创建扰动项不大、基本满足y=2x1−x2+1规律的数据集

# 设置随机数种子
np.random.seed(24)   

# 扰动项取值为0.01
features, labels = arrayGenReg(delta=0.01)

2）设置初始参数

np.random.seed(24)  
w = np.random.randn(3, 1)
w
#array([[ 1.32921217],
#       [-0.77003345],
#       [-0.31628036]])

（3）构建损失函数与梯度表达式

# 计算w取值时SSE
SSELoss(features, w, labels)
#array([[2093.52940481]])
# 计算w取值时MSE
MSELoss(features, w, labels)
#array([[2.0935294]])
# 计算w取值时梯度
lr_gd(features, w, labels)
#array([[-1.15082596],
#       [ 0.3808195 ],
#       [-2.52877477]])

（3）构建损失函数与梯度表达式

# 计算w取值时SSE
SSELoss(features, w, labels)
#array([[2093.52940481]])
# 计算w取值时MSE
MSELoss(features, w, labels)
#array([[2.0935294]])
# 计算w取值时梯度
lr_gd(features, w, labels)
#array([[-1.15082596],
#       [ 0.3808195 ],
#       [-2.52877477]])

（4）执行梯度下降计算
接下来，借助此前所定义的参数更新函数进行梯度下降参数求解：

w = w_cal(features, w, labels, gd_cal = lr_gd, lr = 0.1, itera_times = 100)
w
#array([[ 1.99961892],
#       [-0.99985281],
#       [ 0.99970541]])
# 计算w取值时SSE
SSELoss(features, w, labels)
#array([[0.09300731]])
# 计算w取值时MSE
MSELoss(features, w, labels)
#array([[9.30073138e-05]])

至此，我们就完成了梯度下降算法的手动实现过程。

三、梯度下降算法评价

1. 梯度下降算法优势

以线性回归的损失函数构建了梯度表达式来进行梯度下降的计算，尽管线性回归的损失函数可以用最小二乘法直接进行求解，但梯度下降算法在某些场景下仍然具备一定优势。

这些优势主要体现在两个方面，其一是相比大规模数值矩阵运算，梯度下降所遵循的迭代求解效率更高（尽管大规模矩阵运算也可以通过分块矩阵的划分来减少每一次计算的数据量），其二则是对于某些最小二乘法无法计算全域唯一最优解的情况，梯度下降仍然能够有效进行最小值点（或者解空间）的寻找。

对于通过严格数学公式计算得出的结果，类似最小二乘法计算结果，也被称为解析解。

由某种计算方法或者计算流程得出的结果，例如梯度下降算得的结果，也被称为数值解。

2. 梯度下降算法局限

但是，对于上述梯度下降算法还是存在一定的局限性。当损失函数不是凸函数时，也就是有可能存在局部最小值或者鞍点时，梯度下降并不一定能够准确找到全域最小值。

2.1 局部最小值点（local minimun）陷阱

所谓局部最小值，指的是该点左右两端取值都大于该点，但是该点不是全域最小值点。例如函数：

f(x)=x⋅cos(πx)

该函数函数图像如下：

x = np.arange(-1, 2, 0.1)
y = x * np.cos(np.pi * x)
fig = plt.plot(x, y)[0]
fig.axes.annotate('local minimun', xytext=(-0.77, -1), 
                  arrowprops=dict(arrowstyle='->'), xy=(-0.3, -0.25))
fig.axes.annotate('global minimun', xytext=(0.6, 0.8), 
                  arrowprops=dict(arrowstyle='->'), xy=(1.1, -0.95))
plt.xlabel('x')
plt.ylabel('x·cos(pi·x)')

def f_1(x):
    return (x*np.cos(np.pi*x))
f_1(-1)
#1.0
def f_gd_1(x):
    return (np.cos(np.pi*x)-x*np.pi*(np.sin(np.pi*x)))
f_gd_1(-1)
#-1.0000000000000004

以及相应的梯度更新函数。

def gd_1(lr = 0.02, itera_times = 20, w = -1):
    """
    梯度下降计算函数
    :param lr: 学习率
    :param itera_times：迭代次数
    :param w：参数初始取值
    :return results：每一轮迭代的参数计算结果列表
    """                              
    results = [w]
    for i in range(itera_times):
        w -= lr * f_gd_1(w)              # 梯度计算公式
        results.append(w)
    return results

当我们在以 0.02 作为学习率，参数初始值为 -1，迭代 5000 轮时参数每一轮迭代后参数的计算结果如下：

res = gd_1(itera_times = 5000)
res[-1]
#-0.2738526868008511

我们发现，参数最终在停留在 -0.27 附近。当然，对于上述一维梯度下降，我们也可以绘制对应参数变化的轨迹图来进行观察。

def show_trace_1(res):
    """
    梯度下降轨迹绘制函数
    """
    f_line = np.arange(-1, 2, 0.01)
    plt.plot(f_line, [f_1(x) for x in f_line])
    plt.plot(res, [f_1(x) for x in res], '-o')
    plt.xlabel('x')
    plt.ylabel('Loss(x)')

show_trace_1(res)

我们发现，参数最终停留在局部最小值附近，并且无法跨越该局部最小值点。而根据此前所介绍的梯度下降移动距离衰减理论，最终迭代 5000 轮之后该点梯度也应该是一个趋于 0 的值。当然，这也符合局部最小值的特点，导数为 0。

f_gd_1(res[-1])
#-1.2212453270876722e-15

当然，对于梯度下降的局部最小值陷阱，根本原因还是在于梯度下降过程中参数移动的距离和梯度直接挂钩，而梯度下降的该特点不仅导致了局部最小值陷阱，还有另外一类更加常见的陷阱——鞍点陷阱。

2.2 鞍点（saddle point）陷阱

鞍点，简单来理解就是那些不是极值点但梯度为 0 的点。
极值，指的是那些连续函数上导数为 0、并且所有两边单调性相反的点，极值包括局部最小值、最小值点、局部最大值和最大值点四类。
鞍点和极值点的区别在于导数为 0 点左右两边单调性相同，例如：

f(x)=x3

我们通过绘制该函数的函数图像来观察鞍点。

x = np.arange(-2, 2, 0.1)
y = np.power(x, 3)
fig = plt.plot(x, y)[0]
fig.axes.annotate('saddle point', xytext=(-0.52, -5), 
                  arrowprops=dict(arrowstyle='->'), xy=(0, -0.2))
plt.xlabel('x')
plt.ylabel('x**3')

进一步的，我们来看下鞍点是如何影响梯度下降过程的。首先是由损失函数到梯度计算表达式的过程，由于损失函数只有一个参数，所以梯度表达式为：

f′(x)=3x2

我们可以定义该函数、导函数的函数及轨迹绘制函数。

# x**3函数
def f_2(x):
    return np.power(x, 3)
f_2(-1)
#-1
# x**3导函数
def f_gd_2(x):
    return (3*np.power(x, 2))
f_gd_2(-1)
#3

# 梯度更新函数
def gd_2(lr = 0.05, itera_times = 200, w = 1):
    """
    梯度下降计算函数
    :param lr: 学习率
    :param itera_times：迭代次数
    :param w：参数初始取值
    :return results：每一轮迭代的参数计算结果列表
    """                              
    results = [w]
    for i in range(itera_times):
        w -= lr * f_gd_2(w)              # 梯度计算公式
        results.append(w)
    return results

当我们在以 0.05 作为学习率，参数初始值为 1，迭代 5000 轮时参数每一轮迭代后参数的计算结果如下：

res = gd_2(itera_times=5000)
res[-1]
#0.0013297766246039373

def show_trace_2(res):
    """
    梯度下降轨迹绘制函数
    """
    f_line = np.arange(-1, 2, 0.01)
    plt.plot(f_line, [f_2(x) for x in f_line])
    plt.plot(res, [f_2(x) for x in res], '-o')
    plt.xlabel('x')
    plt.ylabel('Loss(x)')

show_trace_2(res)

同样，梯度下降无法跨越鞍点抵达更小值的点，并且我们判断 5000 次迭代后参数点的梯度已经非常小了：

f_gd_2(res[-1])
#5.304917614029122e-06

我们可以绘制一个包含了鞍点的二维图像进行观察：

x, y = np.mgrid[-1:1:31j,-1:1:31j]
z = x**2- y**2
ax = plt.axes(projection='3d')
ax.plot_wireframe(x, y, z,**{'rstride':2,'cstride':2})
ax.plot([0],[0],[0],'rx')
ticks =[-1,0,1]
plt.xticks(ticks)
plt.yticks(ticks)
ax.set_zticks(ticks)
plt.xlabel('x')
plt.ylabel('y')

值得注意的是，在实际情况中，鞍点出现的频率高于局部最小值点。

3. 梯度下降算法本质与改善方法

根据局部最小值和鞍点的讨论，我们不难发现梯度下降的本质作用其实是让参数点移动到梯度为 0 的点，当损失函数是严格意义的凸函数时，梯度为 0 的点就是全域最小值点，但如果损失函数不是凸函数，那么梯度为 0 的点就有可能是局部最小值点或者鞍点。

此时受到局部最小值点或者鞍点梯度为 0 的影响，梯度下降无法从该点移出。

尽管大多数线性模型的损失函数都是凸函数，但很多复杂机器学习模型所构建的损失函数都不一定是严格凸函数，要避免局部最小值点或者鞍点陷阱，我们就必须在梯度下降算法基础上进行改进。- 有一种最基础也是最通用的改进办法，就是每次在构建损失函数的时候代入一小部分数据，从而让参数有机会跳出陷阱，这就是所谓的随机梯度下降和小批量梯度下降。