4.2　反向传播神经网络两个基础算法详解

在正式介绍BP神经网络之前，需要先介绍两个非常重要的算法，即最小二乘法随机梯度和下降算法。

最小二乘法是统计分析中常用的逼近计算的一种算法，其交替计算结果使得最终结果尽可能地逼近真实结果。而随机梯度下降算法充分利用了深度学习的运算特性的迭代和高效性，通过不停地判断和选择当前目标下的最优路径，使得能够在最短路径下达到最优的结果，从而提高大数据的计算效率。

4.2.1　最小二乘法详解

最小二乘（Least Square, LS）法是一种数学优化技术，也是一种机器学习常用的算法。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

由于最小二乘法不是本章的重点内容，因此笔者只通过一个图示演示一下最小二乘法LS的原理。最小二乘法的原理如图4-5所示。

图4-5　最小二乘法的原理

从图4-5可以看到，若干个点依次分布在向量空间中，如果希望找出一条直线和这些点达到最佳匹配，那么最简单的方法就是希望这些点到直线的值最小，即下面的最小二乘法实现公式的值最小。

这里直接引用真实值与计算值之间的差的平方和，具体而言，这种差值有一个专门的名称——残差。基于此，表达残差的方式有以下3种：

·　∞范数：残差绝对值的最大值为 ，即所有数据点中残差距离的最大值。

·　L1范数：绝对残差和 ，即所有数据点残差距离之和。

·　L2范数：残差平方和 。

可以看到，所谓的最小二乘法也就是L2范数的一个具体应用。通俗地说，就是看模型计算出来的结果与真实值之间的相似性。

因此，最小二乘法的定义如下：

对于给定的数据( x _i , y _i )( i =1,…, m )，在取定的假设空间 H 中，求解 f ( x )∈ H ，使得残差 的L2范数最小。

看到这里，可能有读者会提出疑问，这里的 f ( x )该如何表示？

实际上，函数 f ( x )是一条多项式函数曲线：

由上面的公式可以知道，所谓的最小二乘法，就是找到一组权重 w ，使得 最小。那么问题又来了，如何能使最小二乘法的值最小？

对于求出的最小二乘法的结果，可以使用数学上的微积分处理方法，这是一个求极值的问题，只需要对权值依次求偏导数，最后令偏导数为0，即可求出极值点。

具体实现最小二乘法的代码如下（注意，为了简化起见，本示例使用一元一次方程组来演示拟合）。

【程序4-1】

     import numpy as np
     from matplotlib import pyplot as plt
     A = np.array([[5],[4]])
     C = np.array([[4],[6]])
     B = A.T.dot(C)
     AA = np.linalg.inv(A.T.dot(A))
     l=AA.dot(B)
     P=A.dot(l)
     x=np.linspace(-2,2,10)
     x.shape=(1,10)
     xx=A.dot(x)
     fig = plt.figure()
     ax= fig.add_subplot(111)
     ax.plot(xx[0,:],xx[1,:])
     ax.plot(A[0],A[1],'ko')
     ax.plot([C[0],P[0]],[C[1],P[1]],'r-o')
     ax.plot([0,C[0]],[0,C[1]],'m-o')
     ax.axvline(x=0,color='black')
     ax.axhline(y=0,color='black')
     margin=0.1
     ax.text(A[0]+margin, A[1]+margin, r"A",fontsize=20)
     ax.text(C[0]+margin, C[1]+margin, r"C",fontsize=20)
     ax.text(P[0]+margin, P[1]+margin, r"P",fontsize=20)
     ax.text(0+margin,0+margin,r"O",fontsize=20)
     ax.text(0+margin,4+margin, r"y",fontsize=20)
     ax.text(4+margin,0+margin, r"x",fontsize=20)
     plt.xticks(np.arange(-2,3))
     plt.yticks(np.arange(-2,3))
     ax.axis('equal')
     plt.show()

最终结果如图4-6所示。

图4-6　最小二乘法拟合曲线

4.2.2　道士下山的故事——梯度下降算法

在介绍随机梯度下降算法之前，这里先讲一个道士下山的故事。请读者参考图4-7。

图4-7　模拟随机梯度下降算法的演示图

为了便于理解，我们将其比喻成道士想要出去游玩的一座山。

设想道士有一天和道友一起到一座不太熟悉的山上去玩，在兴趣盎然中很快登上了山顶。但是天有不测，下起了雨。如果这时需要道士和其同来的道友用最快的速度下山，那么怎么办呢？

如果想以最快的速度下山，那么最快的办法就是顺着坡度最陡峭的地方走下去。但是由于不熟悉路，道士在下山的过程中，每走一段路程就需要停下来观望，从而选择最陡峭的下山路。这样一路走下来的话，可以在最短时间内走到山脚。

根据图4-7可以近似地表示为：

① → ② → ③ → ④ → ⑤ → ⑥ → ⑦

每个数字代表每次停顿的地点，这样只需要在每个停顿的地点选择最陡峭的下山路即可。

这就是道士下山的故事，随机梯度下降算法和这个类似。如果想要使用最迅捷的下山方法，那么最简单的方法就是在下降一个梯度的阶层后，寻找一个当前获得的最大坡度继续下降。这就是随机梯度算法的原理。

从上面的例子可以看到，随机梯度下降算法就是不停地寻找某个节点中下降幅度最大的那个趋势进行迭代计算，直到将数据收缩到符合要求的范围为止。通过数学公式表达的方式计算的话，公式如下：

在4.2.1节讲最小二乘法的时候，我们通过最小二乘法说明了直接求解最优化变量的方法，也介绍了求解的前提条件是计算值与实际值的偏差的平方最小。

但是在随机梯度下降算法中，对于系数需要不停地求解出当前位置下最优化的数据。使用数学方式表达的话，就是不停地对系数 θ 求偏导数。公式如下：

在公式中， θ 会向着梯度下降最快的方向减小，从而推断出 θ 的最优解。

因此，随机梯度下降算法最终被归结为：通过迭代计算特征值，从而求出最合适的值。求解 θ 的公式如下：

在公式中，α是下降系数。用较为通俗的话表示，就是用来计算每次下降的幅度大小。系数越大，每次计算中的差值就较大；系数越小，差值就越小，但是计算时间也会相应延长。

随机梯度下降算法的迭代过程如图4-8所示。

图4-8　随机梯度下降算法的迭代过程

从图4-8中可以看到，实现随机梯度下降算法的关键是拟合算法的实现。而本例的拟合算法的实现较为简单，通过不停地修正数据值从而达到数据的最优值。

随机梯度下降算法在神经网络特别是机器学习中应用较广，但是由于其天生的缺陷，噪音较大，使得在计算过程中并不是都向着整体最优解的方向优化，往往只是得到局部最优解。因此，为了克服这些困难，最好的办法就是增大数据量，在不停地使用数据进行迭代处理的时候，能够确保整体的方向是全局最优解，或者最优结果在全局最优解附近。

【程序4-2】

    x = [(2, 0, 3), (1, 0, 3), (1, 1, 3), (1,4, 2), (1, 2, 4)]
    y = [5, 6, 8, 10, 11]
    epsilon = 0.002
    alpha = 0.02
    diff = [0, 0]
    max_itor = 1000
    error0 = 0
    error1 = 0
    cnt = 0
    m = len(x)
    theta0 = 0
    theta1 = 0
    theta2 = 0
    while True:
        cnt += 1
        for i in range(m):
            diff[0] = (theta0 * x[i][0] + theta1 * x[i][1] + theta2 * x[i][2]) - y[i]
            theta0 -= alpha * diff[0] * x[i][0]
            theta1 -= alpha * diff[0] * x[i][1]
            theta2 -= alpha * diff[0] * x[i][2]
        error1 = 0
        for lp in range(len(x)):
            error1 += (y[lp] - (theta0 + theta1 * x[lp][1] + theta2 * x[lp][2])) **
2 / 2
        if abs(error1 - error0) < epsilon:
            break
        else:
            error0 = error1
    print('theta0 : %f, theta1 : %f, theta2 : %f, error1 : %f' % (theta0, theta1,
theta2, error1))
    print('Done: theta0 : %f, theta1 : %f, theta2 : %f' % (theta0, theta1, theta2))
    print('迭代次数: %d' % cnt)

最终结果打印如下：

     theta0 : 0.100684, theta1 : 1.564907, theta2 : 1.920652, error1 : 0.569459
     Done: theta0 : 0.100684, theta1 : 1.564907, theta2 : 1.920652
     迭代次数: 24

从结果来看，这里迭代24次即可获得最优解。

4.2.3　最小二乘法的梯度下降算法以及Python实现

从前面的介绍可以得知，任何一个需要进行梯度下降的函数都可以被比作一座山，而梯度下降的目标就是找到这座山的底部，也就是函数的最小值。根据之前道士下山的场景，最快的下山方式就是找到最为陡峭的山路，然后沿着这条山路走下去，直到下一个观望点。之后在下一个观望点重复这个过程，寻找最为陡峭的山路，直到山脚。

下面带领读者实现这个过程，求解最小二乘法的最小值，但是在开始之前，先介绍读者需要掌握的数学原理。

1．微分

高等数学中对函数微分的解释有很多，主要的有两种：

·　函数曲线上某点切线的斜率。

·　函数的变化率。

因此，对于一个二元微分的计算如下：

2．梯度

所谓的梯度，就是微分的一般形式，多元微分则是各个变量的变化率的总和，例子如下：

可以看到，求解的梯度值是分别对每个变量进行微分计算，之后用逗号隔开。这里用中括号“[]”将每个变量的微分值包裹在一起形成一个三维向量，因此可以认为微分计算后的梯度是一个向量。

综上所述，得出梯度的定义：在多元函数中，梯度是一个向量，而向量具有方向性，梯度的方向指出了函数在给定点上上升最快的方向。

这个与前面道士下山的过程联系在一起，如果需要到达山地，则需要在每一个观察点寻找梯度最陡峭的地方。梯度计算的值是在当前点上升最快的方向，反方向则是给定点下降最快的方向。梯度计算就是得出这个值的具体向量值，如图4-9所示。

图4-9　梯度计算

3．梯度下降的数学计算

前面已经给出了梯度下降的公式，此时对其进行变形：

此公式中的参数的含义如下：

J 是关于参数 θ 的函数，假设当前点为 θ ，如果需要找到这个函数的最小值，也就是山底的话，那么首先需要确定行进的方向，也就是梯度计算的反方向，之后走 α 的步长，走完这个步长之后，就到了下一个观察点。

α 的意义在4.2.2节已经介绍，是学习率或者步长，使用 α 来控制每一步走的距离。 α 过小会造成拟合时间过长，而 α 过大会造成下降幅度太大错过最低点。如图4-10所示为学习率太小（左）与学习率太大（右）的对比。

图4-10　学习率太小与学习率太大的对比

这里还要注意的是，梯度下降公式中的 ∇J(θ) 求出的是斜率最大值，也就是梯度上升最大的方向，而这里所需要的是梯度下降最大的方向，因此在 ∇J(θ) 前加一个负号。下面用一个例子演示梯度下降算法的计算。

假设这里的公式为：

J ( θ )= θ ²

此时的微分公式为：

∇J ( θ )=2 θ

设第一个值 θ ⁰ =1，α=0.3，则根据梯度下降公式：

θ ¹ = θ ⁰ − α ×2 θ ⁰ =1− α ×2×1=1−0.6=0.4

θ ² = θ ¹ − α ×2 θ ¹ =0.4− α ×2×0.4=0.4−0.24=0.16

θ ³ = θ ² − α ×2 θ ² =0.16− α ×2×0.16=0.16−0.096=0.064

这样依次经过运算，即可到 J(θ) 的最小值，也就是“山底”，如图4-11所示。

图4-11 J(θ) 的最小值

实现程序如下：

     import numpy as np
     x = 1
     def chain(x,gama = 0.1):
         x = x - gama * 2 * x
         return x
     for _ in range(4):
         x = chain(x)
         print(x)

多变量的梯度下降算法和前文所述的多元微分求导类似。例如，一个二元函数的形式如下：

此时对其的梯度微分为：

∇J ( θ )=2 θ ₁ +2 θ ₂

此时将设置：

J(θ ⁰ )=(2,5)， α =0.3

则依次计算的结果如下：

剩下的计算请读者自行完成。

如果把二元函数采用图像的方式展示出来，可以很明显地看到梯度下降的每个“观察点”坐标，如图4-12所示。

图4-12　梯度下降的每个“观察点”坐标

4．使用梯度下降算法求解最小二乘法

下面是本节的实战部分，使用梯度下降算法计算最小二乘法。假设最小二乘法的公式如下：

其中参数解释如下：

· m 是数据点总数。

· 是一个常量，这样在求梯度的时候，二次方微分后的结果就与 抵消了，自然就没有多余的常数系数，方便了后续的计算，同时对结果不会有影响。

· y 是数据集中每个点的真实 y 坐标的值。

其中 h _θ (x) 为预测函数，形式如下：

h _θ (x) = θ ₀ + θ ₁ x

根据每个输入 x ，都有一个经过参数计算后的预测值输出。

h _θ (x) 的Python实现如下：

     h_pred = np.dot(x,theta)

其中 x 是输入的、维度为[-1,2]的二维向量，-1的意思是维度不定。这里使用了一个技巧，即将 h _θ (x) 的公式转换成矩阵相乘的形式，而theta是一个[2,1]维度的二维向量。

依照最小二乘法实现的Python代码如下：

     def error_function(theta,x,y):
         h_pred = np.dot(x,theta)
         j_theta = (1./2*m) * np.dot(np.transpose(h_pred), h_pred)
         return j_theta

这里j_theta的实现同样是将原始公式转换成矩阵计算，即：

(h _θ (x)−y) ² =( h _θ (x)−y) ^T ×(h _θ (x)−y)

下面分析一下最小二乘法公式 J(θ) ，此时如果求 J(θ) 的梯度，则需要对其中涉及的两个参数 θ ₀ 和 θ ₁ 进行微分：

下面分别对两个参数的求导公式进行求导：

将分开求导的参数合并，可得新的公式如下：

公式最右边的常数1可以被去掉，公式变为：

使用矩阵相乘表示的公式为：

这里 (x) ^T ×(h _θ (x)−y) 已经转换为矩阵相乘的表示形式。使用Python表示如下：

     def gradient_function(theta, X, y):
         h_pred = np.dot(X, theta) - y
         return (1./m) * np.dot(np.transpose(X), h_pred)

如果读者对np.dot(np.transpose(X), h_pred)理解有难度，可以将公式使用逐个 x 值的形式列出来，这里就不罗列了。

最后是梯度下降的Python实现，代码如下：

     def gradient_descent(X, y, alpha):
         theta = np.array([1, 1]).reshape(2, 1)#[1,1]是theta的初始化参数，后面会修改
         gradient = gradient_function(theta,X,y)
         for i in range(17):
             theta = theta - alpha * gradient
             gradient = gradient_function(theta, X, y)
         return theta

或者使用如下代码：

    def gradient_descent(X, y, alpha):
        theta = np.array([1, 1]).reshape(2, 1) #[1,1]是theta的初始化参数，后面会修改
        gradient = gradient_function(theta,X,y)
        while not np.all(np.absolute(gradient) <= 1e-4):              #采用abs是因为gradient
计算的是负梯度
            theta = theta - alpha * gradient
            gradient = gradient_function(theta, X, y)
            print(theta)
        return theta

这两个代码段的区别在于：第一个代码段是固定循环次数，可能会造成欠下降或者过下降；而第二个代码段使用的是数值判定，可以设定阈值或者停止条件。

全部代码如下：

    import numpy as np
    m = 20
    # 生成数据集x，此时的数据集x是一个二维矩阵
    x0 = np.ones((m, 1))
    x1 = np.arange(1, m+1).reshape(m, 1)
    x = np.hstack((x0, x1)) #[20,2]
    y = np.array([
        3, 4, 5, 5, 2, 4, 7, 8, 11, 8, 12,
        11, 13, 13, 16, 17, 18, 17, 19, 21
    ]).reshape(m, 1)
    alpha = 0.01
    #这里的theta是一个[2,1]大小的矩阵，用来与输入x进行计算，以获得计算的预测值y_pred，而
y_pred用于与y计算误差
    def error_function(theta,x,y):
        h_pred = np.dot(x,theta)
        j_theta = (1./2*m) * np.dot(np.transpose(h_pred), h_pred)
        return j_theta
    def gradient_function(theta, X, y):
        h_pred = np.dot(X, theta) - y
        return (1./m) * np.dot(np.transpose(X), h_pred)
    def gradient_descent(X, y, alpha):
        theta = np.array([1, 1]).reshape(2, 1)  #[2,1]  这里的theta是参数
        gradient = gradient_function(theta,X,y)
    while not np.all(np.absolute(gradient) <= 1e-6):
        theta = theta - alpha * gradient
        gradient = gradient_function(theta, X, y)
    return theta
theta = gradient_descent(x, y, alpha)
print('optimal:', theta)
print('error function:', error_function(theta, x, y)[0,0])

打印结果和拟合曲线请读者自行完成。

现在回到前面的道士下山这个问题，这个下山的道士实际上就代表了反向传播算法，而要寻找的下山路径其实就代表着算法中一直在寻找的参数 θ ，山上当前点的最陡峭的方向实际上就是代价函数在这一点的梯度方向，在场景中观察最陡峭的方向所用的工具就是微分。 hBaPitHOhR+RFEq4Rs1sWMae/pPgn3KvbYqVjv94WUbKQ28rioTG1Z5wuBuyTFid