深度学习代码实践（五）- 从0搭建一个神经网络：从多元方程到梯度下降反向求导

回到多元方程的求解方法

对于构建一个神经网络来说，需要求出每一个神经元的参数。

每一层都有线性变换，加上非线性变换组成。 神经网络的求解变成一个多元方程的求解问题。

图：要求解的线性函数变换（单层）

        y1 ​= x1​w11 ​ + x2​w21​+ ... + xnwn1 + b1

        y2 ​= x1​w12 ​ + x2​w22​ + ... + xnwn2 + b2​

        y3​= x1​w13 ​+ x2​w23 ​+ ... + xnwn3 + b3​

假设要解的这一层输入变量 x 有两个， x1, x2，输出有3个值 y1, y2, y3， 即要解方程

 方程式简写为代数式：

如何理解方程式

理想情况下，如果我们能够找到 一组参数 [[w11, w12, w13],[w21,w22,w23]], 以及一组阈值（称之为偏置） b1, b2, b3 使得对于所有输入的 [x1, x2] 经过上述运算， 能够得到期望的值 y。

那么 [[w11, w12, w13],[w21,w22,w23]], b1, b2, b3 就是我们要找的答案。

简写一下前面这个方程组：

Y1 = X1*W + B

Y2 = X2*W + B

Y3 = X3*W + B

...

Ym = Xm*W + B

总共有 m个样本， 其中 Yi = [yi1, yi2, yi3] ， Xi = [xi1, xi2]，

每一个样本的输入有 xi1, xi2 两个值，输出有 3个维度值。

实际多变量多样本的情况

实际上，在有很多维输入变量，或者有大量样本时，很难计算出来一组准确的 W，B 使得方程式严格相等。往往实际求出来的参数，计算得到的 Y， 跟预期的 Y 会有差异，就是说计算结果存在误差。

退一步的预期：误差最小化（数学优化问题）

假设 y 就是对训练集样本经过计算的结果标签（值）, t 是对于训练样本结果的预期标签（值），理想情况下就是这两个值的误差最小。

误差如何衡量？

回顾一下问题：

输入：28x28 维像素数组，即 784 个变量的输入 [x1, x2, x3, ... x784]，

输出: 对应10个数字概率的向量 y = [z0, z1, z2, z3, z4, ...., z9] ，

zi 是表示这个输入的图像是数字 i 的概率。 比如:

y = [ 0.1, 0.04, 0.2, 0.6, 0, 0, 0, 0.01, 0.02, 0.03]

这个输出里面， z3 = 0.6， 数字3的概率是 0.6，认为识别的图像对应数字3.

用 t 来表示正确解，也用10个元素的向量表示， 如以下 t 表示数字3：

t = [0, 0, 0, 1, 0, 0, 0, 0, 0, 0,]

这里将正确解标记为1， 其他标签标记为0的表示方法，称为 one-hot 表示法。

误差的常用计算方法

均方差误差(Mean Squared Error， MSE)

 

 

yi是神经网络节点计算得到的输出值，ti是样本对应于节点的目标值（监督数据的实际值）。 Ed表示 样本 d 的均方误差。

交叉熵误差(Cross Entropy Loss Function) - 单个样本

 

i表示对应位置的数字，

t = [t1, t2, t3, ... t10] t为实际的分类结果， ti 的值为0或1. 可以理解成概率。

y = [y1, y2, y3, ... y10]  为预测的结果，yi 的值范围 [0,1]， 表示一个概率

def cross_entropy_error(y, t):

delta = le -7

return -np.sum( t * np.log(y+delta) )

为了避免得到 log(0) 负无穷大的值， 函数中加上一个微小值 delta 来做防护。

所有样本的加和交叉熵误差/损失函数

对于 m 个样本，总体的交叉熵误差（m个交叉熵之和，再求平均）

计算示例:

例：需要根据图片动物的轮廓、颜色等特征，来预测动物的类别，有三种可预测类别：猫、狗、猪。假设模型通过sigmoid/softmax的方式得到对于每个预测结果的概率值， 3个输入样本的预测值，真实值如下：

预测(y)	真实(t)	是否正确
0.3 0.3 0.4	0 0 1 (猪)	正确
0.3 0.4 0.3	0 1 0 (狗)	正确
0.1 0.2 0.7	1 0 0 (猫)	错误

对于样本1和样本2以非常微弱的优势判断正确，对于样本3的判断则彻底错误。

现在我们计算前面三个样本的交叉熵损失函数值：

样本三是判断错误的样本，最后算出来的损失最大。

对所有样本的loss求平均：

 

 

以另外1个样本为例：

仿射变换得到 y = [2.0，1.0，0.1]，

经过 softmax 转换得到预测样本是猫、狗、猪的概率值 p = [0.7, 0.2, 01]，

交叉熵的计算过程：

Code - 交叉熵函数 

def cross_entropy_error(y, t):
    if y.ndim == 1:
        t = t.reshape(1, t.size)
        y = y.reshape(1, y.size)

    # 监督数据是one-hot-vector的情况下，转换为正确解标签的索引
    if t.size == y.size:
        t = t.argmax(axis=1)

    bat_size = y.shape[0] # batch size
    return -np.sum(np.log(y[np.arange(bat_size), t] + 1e-7)) / bat_size

问题的转化：如何求解使得误差最小（接近0）的 W, B - 求解函数极值

如何求解得到 W, B， 使得计算得到的 Yr， 跟预期的 Ye的误差最小。

输入的样本数很大时，大量的样本，得到各自的多个Y值，误差有多个，从全局来看，问题是求 W, B，使得所有样本计算结果Y 跟 预期的 Ye 的累积误差最小。

新问题的定义（数学化定义）

同时，我们可以使用估算法/假设法来逐步逼近这个函数， 使得计算出来的值跟预期值误差最小。

单层（一个隐藏层）的问题定义，对于一组输入的 x1, x2, b：

使用交叉熵误差： 

 

 

求解的问题：

对于输入的 m 组样本，每组有已知的 X = [x1, x2]， 以及打过的正确标签 t = [y1, y2, y3] ,

求使得总体的交叉熵误差最小的 W, B，

就是说求 W， B，使得经过函数 X*W + B 计算得到的 y 跟 正确标签 t 的差异最小：

L 为训练集总体的交叉熵（损失函数）， 其中训练集中总共有m个样本， 其中 Ek 为第k个样本的交叉熵。

这样的话，问题变成了一个求函数最小值的问题。

估算法/函数极小值/导数

凸函数（Convex Function）

凸子集C中任意两个向量 x1, x2, 有下式成立：

 

凸函数在某段区间只有一个低谷的函数。对于凸函数，求极小值是一个求导数的过程。

如下二维函数， 函数的极小值出现在切线为0处（即导数 df(y)/dx = 0 处）

 对于三维的函数：

 

名词定义：导数 Derivative

当函数y=f（x）的自变量x在一点x0上产生一个增量Δx时，

函数输出值的增量Δy与自变量增量Δx的比值在Δx趋于0时的极限，a如果存在，a即为在x0处的导数，记作f'(x0）或df(x0)/dx。

一个函数在某一点的导数描述了这个函数在这一点附近的变化率。

如果函数的自变量和取值都是实数的话，函数在某一点的导数即是这一点的斜率。

导数的本质是通过极限的概念对函数进行局部的线性逼近。导数表示某个瞬间的变化率。

名词定义：微分 Differential

y是x的函数 y=f(x) 。

自变数x有微小的变化量时(d/dx)，函数的值 y 也会跟着变动的量(dy)。

自变量的微分记为 dx, 函数值的微分记为 dy

名词定义：偏导数

有多个变量的函数的导数称为偏导数。用数学式表示的话，可以写成 

 例如函数：

 

偏导数的计算方法：

计算偏导数时，将多个变量中的某一个变量定位目标变量（要求偏导数的变量），并将

其他变量固定位某个值。

求 x0=2，x1=3 时，关于 x0 的偏导数 

>>> numerical_diff(function_a, 12)
5.099999999949034
>>> def function_tmp1(x0):
...     return x0*x0 + 3*3
...
>>> numerical_diff(function_tmp1, 2.0)
3.9999999999995595

求 x0=2，x1=3 时，关于 x1 的偏导数  

>>> def function_tmp2(x1):
...     return 2*2 + x1*x1
...

>>> numerical_diff(function_tmp2, 3.0)
6.000000000012662

名词定义：梯度

像 这样的由全部变量的偏导数汇总而成的向量，称之为梯度（Gradient）。

梯度表示的是各点处的函数值减小最多的方向。

随机梯度下降 - 求损失函数的最小值

定义：从任意一个随机点开始，沿着当前梯度的父方向， 迭代更新权重参数，找到目标函数最小值的方法。

对于三维空间的直观理解：

从一座大山上的某处位置（起点）开始下山（目标是到达山底），由于我们不知道怎么下山，于是决定走一步算一步。在每走到一个位置的时候，求解当前位置的梯度，沿着梯度的负方向，也就是当前最陡峭的位置向下走一步。然后继续求解当前位置梯度，向这一步所在位置沿着最陡峭最易下山的位置走一步。这样一步步的走下去，一直走到觉得我们已经到了山脚。

随机梯度下降计算过程

 

1. 随机取一个自变量的值 x0，如图中点 A；
2. 对应该自变量算出对应点的函数值：f(x0)；
3. 计算 f(x0) 处目标函数 f(x) 的导数。 对于一元的函数，图中A点的导数，大约为 -1.4 (斜率）, d(f(x))/dx = -1.4 （x 增加 0.01 时，y 大约增加 -0.014）
4. 从 f(x0) 开始，沿着该处目标函数导数的反方向，按一个指定的步长 delta，向前“走一步”到B点，B点位置对应的自变量 x1 = x0 - delta * d(f(x))/dx
5. 继续重复2-4，直至退出迭代（达到指定迭代次数，或 f(x) 近似收敛到最优解，使得 d(f(x))/dx = 0）。

这里的 f(x) 不仅仅适用于一元函数，也适用于我们使用前面分析过的交叉熵函数。 因此可以用梯度下降的方法找到使得误差最下的参数组 W, B。

这里5个步骤的过程中，最重要的是求导数。 如何求取导数，有不同的方法。 比如可以用数值微分法（死算），也可以用微积分的链式求导法则（快速计算），使用反向传播的方法来求导。 我们先来总结、定义一下要求导的函数。

对于下山的过程来说，由于是在一个三维空间，在任何一点都有相对于三个维度的梯度值。梯度是一个向量（不是标量）。

问题进一步聚焦 - 线性方程求导数

单层（一个隐藏层）的问题定义，对于一组输入的 x1, x2, b：

求解的问题：

对于输入的 m 组样本，每组有已知的 X = [x1, x2]， 以及打过的正确标签 t = [y1, y2, y3] ,

求使得总体的交叉熵误差（L）最小的 W, B，

就是说求 W， B，使得经过函数 X*W + B 计算得到的 y 跟 正确标签 t 的差异最小：

其中训练集中总共有m个样本， 其中 Ek 为第k个样本的交叉熵。

因此求最小值的问题，变成了 L 对 W, B 求导。

注意：

1.对于从训练集来找最优的 W, B 这个优化问题来说， x 是已知的，函数其实变成了 y 对 w, b 的函数， L 对 w, b 的函数。 不再是对 x 的函数。

2.如下图，多维变量 W 是一个矩阵， 因此 L 对于 W 的梯度，也是一个矩阵。就是说 L 对于 w11, w21, w12, w22, w13, w23 分别有一个梯度。 同样对于 b1, b2, b3 也都有各自的梯度。

 

 

经过前面的分析，我们通过这些步骤来迭代，梯度下降找到最优的 W， B

1. 确认模型特征： 样本维度（X的维度），输出的特征维度（Y的维度），以及损失函数（如交叉熵），设定一个合适的步长 delta， 如 delta = 0.01 ；算法的终止距离 d，比如 d = 0.001（算法终止距离：迭代计算得到 损失函数值，与上一次的损失函数值之差）
2. 参数初始化：随机初始化分配 W, B 的值 （前面曲线中任意找到一点）；
3. 确定当前位置的损失函数的梯度：从这一个点(W,B)， 计算 L 对 W, L 对 B 的导数（梯度）
4. 用步长 delta 乘以损失函数的梯度，得到当前位置下降的距离；
5. 更新所有的 参数 W, B；
6. 检查 使用更新的 W，B 计算得到的损失函数值，梯度下降的距离是否都小于 d，或者超过了最大的迭代次数 。如果满足任一条件，则算法终止，当前所有的 参数 W，B 即为最优解的最终结果。否则进入步骤5，继续更新。

对于求导问题的进一步定义

输入：训练集多个样本的 X, Y, 训练集各个样本对应的标签 t， 以及当前所在点的 W, B。

要求的导数： L 对于 W, B 的导数。

 

求导（L对W, B的导数）的多种方法

1. 数值微分法（numerical differentiation）
2. 符号微分法（symbolic differentiation）：即解析性求导，如 对于 y = x^{2} , dy/dx = 2x
3. 自动微分法（automatic differentiation）：使用链式法则求导

数值微分法计算导数

 

函数y=f(x) 在点 x 处的导函数: 

h 也写作∆x, 希腊字母∆读作delta，表示趋近于0的微小变化。

这个函数称之为向前差商。还有向后差商，中心差商。

实际上，按照前面这个函数，求出来的导数，实际上是 PQ 两个点中间点的近似导数。 如下图 P 点，斜率 l 是更准确的 P 点的导数，这个导数可以用中心差商算出来。

 

 中心差商求导

 下面来实现这个函数

def numerical_diff(f, x):
    h = 1e-4 # 0.0001
    return (f(x+h) - f(x-h)) / (2*h)

代码比较简单，入参 f 是一个函数。 x 是输入的变量。

numerical_diff 表示数值微分的含义。

举例， 对于函数 y = 0.2x^2 + 0.3x, 求 x = 5, x=12 处的导数

>>> import numpy as np
>>> import matplotlib.pylab as plt
>>> def function_a(x):
...     return 0.2*x*x + 0.3*x

>>> x = np.arange(0.0, 30, 0.1) # 以0.1 为单位, 从0 到 30的数组x
>>> y = function_a(x)
>>> plt.xlabel("x")
Text(0.5, 0, 'x')
>>> plt.ylabel('f(x)')
Text(0, 0.5, 'f(x)')

>>> plt.plot(x,y)
[<matplotlib.lines.Line2D object at 0x111920ca0>]
>>> plt.show()

>>> def numerical_diff(f, x):
...     h = 1e-4 # 0.0001
...     return (f(x+h) - f(x-h)) / (2*h)
...
>>> numerical_diff(function_a, 5)
2.2999999999928633
>>> numerical_diff(function_a, 12)
5.099999999949034

数值微分：求梯度 （求多维变量x的梯度） 

def function_2(x):
    return x[0]**2 + x[1]**2

def numerical_gradient(f, x):
    h = 1e-4 # 0.0001
    grad = np.zeros_like(x) # 生成和x形状相同的数组
    for idx in range(x.size):
        tmp_val = x[idx]
        
        # f(x+h)的计算
        x[idx] = tmp_val + h
        fxh1 = f(x)
        
        # f(x-h)的计算
        x[idx] = tmp_val - h
        fxh2 = f(x)
        
        grad[idx] = (fxh1 - fxh2) / (2*h)
        x[idx] = tmp_val # 还原值
    return grad

求点 (3, 5)、(3, 0)、(4, 2) 处的梯度

>>> numerical_gradient(function_2, np.array([3, 5]))
array([25000, 45000])
>>> numerical_gradient(function_2, np.array([3, 0]))
array([25000,     0])
>>> numerical_gradient(function_2, np.array([4, 2]))
array([35000, 15000])

数值微分的梯度下降实现 （求解函数 y = x0^2 + x1^2 的最小值处的 x ）

def gradient_descent(f, init_x, lr=0.01, step_num=100):
    ## 任意设定的随机初始值 init_100x
    x = init_x
    for i in range(step_num):
        ## 在 x 处的数值梯度， 沿着梯度下降的方向改变x，执行 step_num 次
        grad = numerical_gradient(f, x)
        x -= lr * grad

    return x

>>> init_x = np.array([-3.0, 4.0])
>>> gradient_descent(function_2, init_x=init_x, lr=0.1, step_num=100)
array([-6.11110793e-10,  8.14814391e-10])

如上梯度下降求解到的 (-6.1e-10, 8.1e-10) 非常接近于0。

反向传播求导数/链式法则求导（计算图求导）

为什么使用反向传播求导（自动求导）？如果使用数值微分的方法来求导，因为有大量的样本，大量的参数，数值微分的方法，需要经过大量的运算。

为什么反向传播能够工作？ 实际上是利用了前向传播过程的计算结果，计算过程，避免了实质上是重复的计算过程。

链式法则（包含正向/反向传播)求导/自动微分法依赖的假设：

所有函数都是由一系列简单的基本操作组合而成。

仅针对基本操作进行符号微分，基于基本操作的微分使用微积分中的链式法则综合求解整个函数的导数，保存中间的数值结果，并不会导出一个封闭形式的公式。

复合函数是由多个函数构成的函数， 对于函数 

 

当求z对x的偏导数时，将其它变量(y)看作常数。 用换元得到下式：

计算图的反向传播

若y = f(x),则

反向传播的计算顺序是，将信号E乘以节点的局部导数  ，

然后将结果传递给下一个节点。这里所说的局部导数是指正向传播中y = f(x)的导数，

也就是y关于x的导数  。一次计算下去，可以计算出梯度。

 

加法节点的反向传播将上游的值原封不动地输出到下游 (z = x + y)。

 

乘法的反向传播会乘以输入信号的翻转值(z=x*y) 。

使用链式法则来求导

 

由于 y1 = x1x11 + x2w21 + b，

损失函数 L 对 w11 的导数只与 y1, z1 有关，跟 y2, y3, z2, z3 都没有关系。

手工推到 L 对 w11 的导数

除了使用链式法则来直接推导反向传播的求导， 也可以用计算图的方法来理解反向传播求导。

 

Code - 反向传播求导 

# coding: utf-8
import sys, os
sys.path.append(os.pardir)  # 为了导入父目录的文件而进行的设定
from common.functions import *
from common.gradient import numerical_gradient


class TwoLayerNet:

    def __init__(self, input_size, hidden_size, output_size, weight_init_std=0.01):
        # 初始化权重
        self.params = {}
        self.params['W1'] = weight_init_std * np.random.randn(input_size, hidden_size)
        self.params['b1'] = np.zeros(hidden_size)
        self.params['W2'] = weight_init_std * np.random.randn(hidden_size, output_size)
        self.params['b2'] = np.zeros(output_size)

    def predict(self, x):
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']

        a1 = np.dot(x, W1) + b1
        z1 = sigmoid(a1)
        a2 = np.dot(z1, W2) + b2
        y = softmax(a2)

        return y

    # x:输入数据, t:监督数据
    def loss(self, x, t):
        y = self.predict(x)

        return cross_entropy_error(y, t)

    def accuracy(self, x, t):
        y = self.predict(x)
        y = np.argmax(y, axis=1)
        t = np.argmax(t, axis=1)

        accuracy = np.sum(y == t) / float(x.shape[0])
        return accuracy

    # 数值微分梯度下降， x:输入数据, t:监督数据
    def numerical_gradient(self, x, t):
        loss_W = lambda W: self.loss(x, t)

        grads = {}
        grads['W1'] = numerical_gradient(loss_W, self.params['W1'])
        grads['b1'] = numerical_gradient(loss_W, self.params['b1'])
        grads['W2'] = numerical_gradient(loss_W, self.params['W2'])
        grads['b2'] = numerical_gradient(loss_W, self.params['b2'])

        return grads

    # 反向传播梯度下降
    def gradient(self, x, t):
        W1, W2 = self.params['W1'], self.params['W2']
        b1, b2 = self.params['b1'], self.params['b2']
        grads = {}

        batch_num = x.shape[0]

        # forward
        a1 = np.dot(x, W1) + b1
        z1 = sigmoid(a1)
        a2 = np.dot(z1, W2) + b2
        y = softmax(a2)

        # backward
        dy = (y - t) / batch_num
        grads['W2'] = np.dot(z1.T, dy)
        grads['b2'] = np.sum(dy, axis=0)

        da1 = np.dot(dy, W2.T)
        dz1 = sigmoid_grad(a1) * da1
        grads['W1'] = np.dot(x.T, dz1)
        grads['b1'] = np.sum(dz1, axis=0)

        return grads

训练神经网络，分别计算训练集，测试集的准确率

 

# coding: utf-8
 import sys, os
 sys.path.append(os.pardir)  # 为了导入父目录的文件而进行的设定
 import numpy as np
 import matplotlib.pyplot as plt
 from dataset.mnist import load_mnist
 from two_layer_net import TwoLayerNet

 # 读入数据
 (x_train, t_train), (x_test, t_test) = load_mnist(normalize=True, one_hot_label=True)

 #network = TwoLayerNet(input_size=784, hidden_size=50, output_size=10)
 network = TwoLayerNet(input_size=784, hidden_size=20, output_size=10)

 iters_num = 10000  # 适当设定循环的次数
 train_size = x_train.shape[0]
 batch_size = 100
 learning_rate = 0.1

 train_loss_list = []
 train_acc_list = []
 test_acc_list = []

 iter_per_epoch = max(train_size / batch_size, 1)

 for i in range(iters_num):
     batch_mask = np.random.choice(train_size, batch_size)
     x_batch = x_train[batch_mask]
     t_batch = t_train[batch_mask]

     # 计算梯度
     grad = network.numerical_gradient(x_batch, t_batch)
     #grad = network.gradient(x_batch, t_batch)

     # 更新参数
     for key in ('W1', 'b1', 'W2', 'b2'):
         network.params[key] -= learning_rate * grad[key]

     loss = network.loss(x_batch, t_batch)
     train_loss_list.append(loss)

     if i % iter_per_epoch == 0:
         train_acc = network.accuracy(x_train, t_train)
         test_acc = network.accuracy(x_test, t_test)
         train_acc_list.append(train_acc)
         test_acc_list.append(test_acc)
         print("train acc, test acc | " + str(train_acc) + ", " + str(test_acc))

# 绘制图形
 markers = {'train': 'o', 'test': 's'}
 x = np.arange(len(train_acc_list))
 plt.plot(x, train_acc_list, label='train acc')
 plt.plot(x, test_acc_list, label='test acc', linestyle='--')
 plt.xlabel("epochs")
 plt.ylabel("accuracy")
 plt.ylim(0, 1.0)
 plt.legend(loc='lower right')
 plt.show()

分别运行数值微分梯度下降，以及反向传播梯度下降的性能。

经过多次迭代，准确率越来越高，在测试集，训练集上面都达到了 94% 的准确率。

 

训练集、测试集上面的准确度趋势图

 

到这里为止， 从头到位实现了一个神经网络，并用于手写数字体的识别。 完整代码：

GitHub - davideuler/beauty-of-math-in-deep-learning: Beauty of math in deep learning

名词翻译

error function: 误差函数，对于单个样本预测值和实际值的偏差

loss function/cost function：损失函数/代价函数，指的是同一个含义。

weight/paramater: 深度学习中一般用权重 weight, 机器学习中一般用参数parameter， 同一含义。

active functiontransfer function/: 激活函数active function和转移函数transfer function， 同一个含义，都是叠加的非线性函数的说法。

Perceptron: 感知机/感知器 

参考

1.《深度学习入门 - 基于 Python 的理论与实现》

2.Learning representations by back-propagating errors-[BP系列]

[BP系列]-Learning representations by back-propagating errors_Dream__Zh的博客-CSDN博客

3.What is Deep Learning and How does it work?

https://towardsdatascience.com/what-is-deep-learning-and-how-does-it-work-2ce44bb692ac

4.零基础入门深度学习(3) - 神经网络和反向传播算法

零基础入门深度学习(3) - 神经网络和反向传播算法 - 作业部落 Cmd Markdown 编辑阅读器

5.神经网络浅讲：从神经元到深度学习

神经网络浅讲：从神经元到深度学习 - 计算机的潜意识 - 博客园

6.误差反向

深度学习入门之4--误差反向传播法_代码哥19950715的博客-CSDN博客

7.损失函数｜交叉熵损失函数

损失函数｜交叉熵损失函数 - 知乎

8.交叉熵误差

交叉熵误差（cross entropy error)_布鲁克林有一棵树-CSDN博客_交叉熵误差