吴恩达机器学习笔记（8）—神经网络：反向传播算法（附代码）

上一节介绍了神经网络模型可以利用前向传播做预测，并引入了一系列符号来描述神经网络模型：

• $x_i$ 表示输入层的第 $i$ 个输入特征，其中 $x_0$ 偏置项省略；
• $a_i^{(j)}$ 表示第 $j$ 层的第 $i$ 个激活单元，其中 $a_0^{(j)}$ 偏置单元省略；
• $s_j$ 表示第 $j$ 层的神经元数量（不包含偏置项），$a^{(j)}\in {\mathbb{R}}^{s_j+1}$ 表示加上偏置项后的 $a^{(j)}$；
• ${\Theta }^{(j)}\in {\mathbb{R}}^{s_{j+1}\times (s_j+1)}$ 表示第 $j$ 层到第 $j+1$ 层的权重矩阵，$\Theta$ 表示所有权重矩阵的集合。

其中的权重参数，与前面的回归问题类似，需要通过最小化代价函数来逼近。

一、代价函数

为了描述代价函数，我们引入新的符号：

• $L$ 表示神经网络的层数，包含输入层、隐藏层、输出层，因此共有 $L-1$ 次传播；
• $K⩾3$ 表示多分类问题的类别数，相应的有 $s_L=K$，$h_{\theta }(x)\in {\mathbb{R}}^K$；
• $h_{\theta }(x^{(i)})\in {\mathbb{R}}^K$，表示从第 $i$ 个样本经过神经网络得到的输出结果，$h_{\theta }(x^{(i)}{)}_k$ 表示其第 $k$ 维。

回顾正则化逻辑回归中的交叉熵代价函数：
$$J(\theta )=-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}\log h_{\theta }(x^{(i)})+(1-y^{(i)})\log (1-h_{\theta }(x^{(i)}))\right]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$

在神经网络中，我们沿用交叉熵代价函数而不是均方误差，原因如下：

1. 我们使用了 sigmoid 函数作为激活函数引入非线性，其单增 $S$ 形曲线使其在 $y=0$ 和 $1$ 附近导数较小；
2. 通常我们会随机初始化参数后进行梯度下降，如果使用均方误差代价函数，对 $\theta$ 求偏导后含有 sigmoid 的导数项，会使前期的梯度下降十分缓慢；
3. 如果使用交叉熵代价函数，则求偏导后不含 sigmoid 的导数项，前期收敛更快。

具体地，对于一个 $K$ 分类的神经网络而言，可以看作是 $K$ 个二分类的交叉熵代价函数之和：
$$J(\Theta )=-\frac{1}{m}\left[\sum _{i=1}^m\sum _{k=1}^K{y}_k^{(i)}\ln {\left(h_{\Theta }(x^{(i)})\right)}_k+(1-y_k^{(i)})\ln {\left(1-h_{\Theta }(x^{(i)})\right)}_k\right]+\frac{\lambda }{2m}\sum _{l=1}^{L-1}\sum _{i=1}^{s_l}\sum _{j=1}^{s_{l+1}}{\left({\Theta }_{ji}^{(l)}\right)}^2$$

注意这里有一个符号混用的地方：上标 $(i)$ 在 $x^{(i)}$ 和 $y^{(i)}$ 中表示第 $i$ 个数据的输入 $(\in {\mathbb{R}}^{n+1})$ 和输出 $(\in {\mathbb{R}}^K)$ ，但是在 ${\Theta }^{(l)}$ 中表示第 $l$ 层的 $\Theta (\in {\mathbb{R}}^{s_{l+1}\times (s_l+1)})$ 。

二、反向传播算法

现在我们的目标就是最小化代价函数 $J(\Theta )$ ，并找到使之最小的参数 $\Theta$ ，我们仍使用梯度下降法解决最小化问题（尽管 $J(\Theta )$ 可能为非凸函数）。因此需要计算梯度 $\frac{\partial J}{\partial {\Theta }_{ji}^{(l)}}$ ，而计算的关键就是「反向传播」：从输出层开始逐层向前计算。

• 为什么需要「反向传播」？我们知道梯度反映的是 $J$ 在 ${\Theta }_{ji}^{(l)}$ 方向的变化率，即：当参数发生微小增量时，代价函数的增量。对于前几层网络，参数的增量需要经过多层全连接传播，才会影响到输出（代价函数）。
• 如果我们从输入层开始计算每个参数的梯度，则每个参数都要依赖其后续路径的传播，每个参数的计算都要先算一遍后续路径的梯度。而如果从输出层开始逐层向前计算并存储了梯度，则解耦了依赖关系，前一层可以利用后一层的计算结果。有点类似「动态规划」中「递归求解子问题」的思想。

1、数学推导：问题转化

为了推导方便，先假设只有一组数据 $(x,y)$ ，这样可以省去上标和求和的麻烦。并且令 $\lambda =0$ 省略正则化项，等到最后的偏导结果再加上。

我们首先对第 $l+1$ 层的第 $j$ 个神经元（偏置单元除外）定义一个 delta 误差：
$${\delta }_j^{(l+1)}=\frac{\partial J}{\partial z_j^{(l+1)}}$$

其中 $z_j^{(l+1)}={\sum }_{k=0}^{s_l}{\Theta }_{jk}^{(l)}{a}_k^{(l)}$ ，即 $a_j^{(l+1)}$ 神经元取 sigmoid 激活前的输出值。换句话说，$a_j^{(l+1)}=g\left(z_j^{(l+1)}\right)$ 。

于是根据链式求导法则，我们有：
$$\frac{\partial J}{\partial {\Theta }_{ji}^{(l)}}=\frac{\partial J}{\partial z_j^{(l+1)}}\cdot \frac{\partial z_j^{(l+1)}}{\partial {\Theta }_{ji}^{(l)}}={\delta }_j^{(l+1)}{a}_i^{(l)}$$

写作矩阵形式：
$${\Delta }^{(l)}={\delta }^{(l+1)}\cdot {\left(a^{(l)}\right)}^T$$

显然，一旦求出当前层的 $delta$ 误差，那么传入当前层的各参数的梯度也可求出。所以现在问题转化为求解 ${\delta }^{(l+1)}$ 。

2、数学推导：输出层误差

既然要递归求解子问题，首先就是要算输出层，即 ${\delta }_j^{(L)}$ ：
$$\begin{array}{rl}{\delta }_j^{(L)}& =\frac{\partial J}{\partial z_j^{(L)}}=\frac{\partial J}{\partial a_j^{(L)}}\cdot \frac{\partial a_j^{(L)}}{\partial z_j^{(L)}}\\ & =\frac{\partial J}{\partial a_j^{(L)}}\cdot g^{\prime }\left(z_j^{(L)}\right)\\ & =-\left(\frac{y_j}{a_j^{(L)}}-\frac{1-y_j}{1-a_j^{(L)}}\right)\cdot \left(a_j^{(L)}\left(1-a_j^{(L)}\right)\right)\\ & =a_j^{(L)}-y_j\end{array}$$

写作矩阵形式：
$${\delta }^{(L)}=a^{(L)}-y$$

• 此处推导过程的一些注释：
• 关于第二行，请注意：$a^{(L)}=g(z^{(L)})$ ；
• 关于第三行，请注意 sigmoid 函数的性质：$g^{\prime }(z)=g(z)(1-g(z))$ ；
• 以及第三行偏导项的计算，请注意 $a^{(L)}=h_{\Theta }(x)$ ，所以 $J(\Theta )$ 在现在的假设条件下可以写作：
  $$J(\Theta )=-\left[\sum _{k=1}^K{y}_k\ln \left(a_k^{(L)}\right)+(1-y_k)\ln \left(1-a_k^{(L)}\right)\right]$$

3、数学推导：隐藏层误差

下面计算第 $l$ 层 $(2⩽l<L)$ 的 ${\delta }_j^{(l)}$：
$$\begin{array}{rl}{\delta }_j^{(l)}& =\frac{\partial J}{\partial z_j^{(l)}}\\ & =\frac{\partial J}{\partial z_j^{(l+1)}}\cdot \frac{\partial z_j^{(l+1)}}{\partial a_j^{(l)}}\cdot \frac{\partial a_j^{(l)}}{\partial z_j^{(l)}}\\ & ={\delta }^{(l+1)T}\cdot {\Theta }_j^{(l)}\cdot g^{\prime }\left(z_j^{(l)}\right)\\ & ={\delta }^{(l+1)T}\cdot {\Theta }_j^{(l)}\cdot \left(a_j^{(l)}\left(1-a_j^{(l)}\right)\right)\end{array}$$

观察式子的前半部分，可以看出下一层所有结点的 $delta$ 误差均会影响当前层的任一结点，而影响则是通过该结点链出的权重参数反向传播。

转置后，写作矩阵形式：
$${\delta }^{(l)}=\left({\Theta }^{(l)T}{\delta }^{(l+1)}\right)\odot \left(a^{(l)}\odot \left(1-a^{(l)}\right)\right)$$

其中 $\odot$ 表示 Hadamard 积，即两个向量对应位置相乘。

4、步骤总结

以上是对一组数据的推导，我们得到了三个重要的结果 (1)(2)(3)：
$$\{\begin{array}{ll}{\Delta }^{(l)}={\delta }^{(l+1)}\cdot {\left(a^{(l)}\right)}^T& 1⩽l<L\\ {\delta }^{(L)}=a^{(L)}-y& \\ {\delta }^{(l)}=\left({\Theta }^{(l)T}{\delta }^{(l+1)}\right)\odot \left(a^{(l)}\odot \left(1-a^{(l)}\right)\right)& 2⩽l<L\end{array}$$

而 $m$ 组数据只需要在一些地方进行累加即可。设数据集为 $\{(x^{(i)},y^{(i)})\mid 1⩽i⩽m\}$ ，则反向传播算法的步骤如下：

1. 所有 ${\Delta }^{(l)}$ 置零；

2. 遍历数据集，设当前数据为 $(x^{(i)},y^{(i)})$ ：
   i. 以 $x^{(i)}$ 为输入做前向传播，得到输出 $a^{(L)}$ ；
   ii. 公式 (2) 计算输出层误差，公式 (3) 计算隐藏层误差；
   iii. 公式 (1) 更新各层的 ${\Delta }^{(l)}$ 矩阵，进行累加：${\Delta }^{(l)}:={\Delta }^{(l)}+{\delta }^{(l+1)}\cdot {\left(a^{(l)}\right)}^T$；

3. 计算 $D$ 矩阵，对偏置项以外的参数进行正则化：
   $$D_{ij}^{(l)}:=\{\begin{array}{ll}\frac{1}{m}\left({\Delta }_{ij}^{(l)}+\lambda {\Theta }_{ij}^{(l)}\right)& \text{if }j\ne 0\\ \frac{1}{m}{\Delta }_{ij}^{(l)}& \text{if }j=0\end{array}$$

这就是最终的梯度矩阵：$\frac{\partial J}{\partial {\Theta }_{ij}^{(l)}}=D_{ij}^{(l)}$ 。

现在，我们可以用 $D_{ij}^{(l)}$ 做一次梯度下降了，整个步骤称为一代（epoch）。注意，我们的参数 ${\Theta }^{(l)}$ 应该初始化为 $[-ϵ,ϵ]$ 中的随机值。

三、参数展开

实现过程中，由于每一层都有一个 ${\Theta }^{(l)}$ 和 $D^{(l)}$ 矩阵，总共有 $2\times L$ 个矩阵，为了方便存储和计算，我们最好都转换为向量，并将各层的向量拼接成一个长向量。这样做可以方便调用封装好的高级优化函数接口，如 fminunc 或 scipy.optimize.minimize。

同样，函数输出的结果向量也可以通过 reshape 转为矩阵，再应用于前向传播。

四、梯度检验

当我们为一个较为复杂的模型（例如神经网络）实现梯度下降算法时，可能会出现一些不易察觉的错误，导致虽然代价看上去在不断减小，但最终的结果可能并不是最优解，误差可能高出一个量级。

为了避免这样的问题，可采取梯度的数值检验方法：通过数值计算得到梯度的近似值，再和反向传播的结果比对，检验是否接近。

在任意一代计算后，得到一组 $\theta$ 向量及 $D$ 矩阵，通过差商近似 $J(\theta )$ 的各偏导：
$$\frac{\partial }{\partial {\theta }_k}J(\theta )\approx \frac{J({\theta }_1,\cdots ,{\theta }_k+ϵ,\cdots ,{\theta }_n)-J({\theta }_1,\cdots ,{\theta }_k-ϵ,\cdots ,{\theta }_n)}{2ϵ}$$

其中，$ϵ$ 是非常小的常数，为了避免计算误差，通常选取 $10^{-4}$。现在我们可以比较这些偏导估计值与对应位置 $D_{ij}^{(l)}$ 的值，它们应该非常接近。

注意：梯度检验耗时巨大，复杂度远大于反向传播，一旦验证了神经网络反向传播的代码正确后，不应进行梯度检验（删掉 or 注释掉）。

五、参数随机初始化

前文我们提到参数 ${\Theta }^{(l)}$ 应该初始化为 $[-ϵ,ϵ]$ 中的随机值，而不是简单初始化为全零。如果将 ${\Theta }^{(l)}$ 初始化为全零，则下一层的 $z^{(l+1)}$ 也为全零，则 $a^{(l+1)}$ 就全为 $0.5$，且所有隐藏层都会得到这个结果；与此同时，反向传播公式 (3) 中，由输出层到隐藏层时 $delta$ 误差也会清零。这将导致网络无法传播！

那么，如果将参数 ${\Theta }^{(l)}$ 全都初始化为同一个非零常量呢？由于全连接，同一层的 $a^{(l+1)}$ 将为相同值；在反向传播公式 (3) 中，由输出层到隐藏层时 $delta$ 误差也会变成相同值，则所有的 $D^{(l)}$ 矩阵也为常数矩阵，参数 ${\Theta }^{(l)}$ 将再次被更新为新的相同常量。这样下去所有的激活单元都成了摆设，因为每一层都在重复计算着相同特征。

1. 上述问题也被称为对称权重问题，解决方法就是随机初始化，实现时可以将 $[0,1]$ 的随机数映射到 $[-ϵ,ϵ]$ ，这里的 $ϵ$ 与梯度检验中的无关。

2. 实际应用中 $ϵ$ 的选取，通常与 ${\Theta }^{(l)}$ 前后两层的神经元个数有关，一种常用的有效的取值是：
   $${ϵ}^{(l)}=\frac{\sqrt{6}}{\sqrt{s^{(l)}+s^{(l+1)}}}$$

六、代码实现

下面以 Coursera 上的多分类数据集 ex4data1.mat 为例，这是一个手写数字的数据集，与上一节的数据集类似。共 5000 组数据，每组数据输入是一个由 $20\times 20$ 灰度矩阵压缩而来的 400 维向量，输出是 0 到 9 之间的整数。

题目还提供一组训练好的权重参数 ex4weight.mat，用以检验前向传播和代价函数的正确性，此处跳过。采用题目推荐的网络层数 $L=3$，则输入层有 401 个单元（含偏置项），隐藏层有 26 个单元（含偏置项），输出层有 10 个单元（独热编码）。

实现反向传播算法并检验结果，大约耗时 20 分钟：

# 导入必要的库
import numpy as np                  # 用于数值计算和矩阵操作
import scipy.io as scio             # 用于读取MATLAB格式的数据文件
from sklearn.metrics import classification_report  # 用于生成分类结果的详细报告

# 加载数据集
# 从MAT文件中读取数据，该文件包含手写数字识别的训练样本
data = scio.loadmat('ex4data1.mat')
X = data['X']                       # 提取特征数据，形状为(5000, 400)，每一行是一个20×20像素的手写数字展开成的向量
y = data['y'].flatten()             # 提取标签数据并将其展平为一维数组(5000, )
y[y==10] = 0                        # 原数据中数字0用10表示，这里将其转换为0以便统一处理

# 神经网络参数设置
(m, n) = (5000, 401)                # m: 样本数量，n: 特征数量+1(偏置项)
L = 3                               # 神经网络的层数（输入层+隐藏层+输出层）
s = [0, 400, 25, 10]                # 各层神经元数量，s[1]为输入层(400个特征)，s[2]为隐藏层(25个神经元)，s[3]为输出层(10个类别)
lmd = 1                             # 正则化参数，用于防止过拟合
alpha = 0.1                         # 梯度下降的学习率
num_iters = 5000                    # 梯度下降的迭代次数
#J_history = []                     # 用于记录每次迭代的损失值（当前注释掉）

# 对标签进行One-hot编码
# 将形状为(m, )的标签转换为(m, 10)的矩阵，每一行只有一个1，表示对应类别
Y = np.zeros((m, s[L]))
for i in range(m):
    Y[i][y[i]] = 1

# 随机初始化权重参数Theta
# Theta[l]表示从第l层到第l+1层的权重矩阵
Theta = [None] * L                  # 索引0未使用，符合层编号习惯
for l in range(1, L):
    # 生成(s[l+1], s[l]+1)形状的随机矩阵，+1是因为要包含偏置项的权重
    Theta[l] = np.random.rand(s[l+1], s[l]+1)
    # 使用Xavier初始化方法调整权重范围，有助于训练收敛
    eps = np.sqrt(6) / np.sqrt(s[l+1] + s[l])
    Theta[l] = Theta[l] * 2 * eps - eps  # 将权重调整到[-eps, eps]范围内

def sigmoid(z):
    """sigmoid激活函数，将输入映射到(0,1)区间"""
    return 1 / (1 + np.exp(-z))

def forwardProp():
    """前向传播计算
    返回各层的激活值a，其中a[l]表示第l层的输出
    """
    a = [None] * (L+1)  # a[0]未使用，索引从1开始对应各层
    
    # 输入层处理：添加偏置项(全为1的列)
    a[1] = np.c_[np.ones(m), X]  # 形状为(5000, 401)
    
    # 隐藏层计算
    for l in range(2, L):
        # 计算加权输入并应用sigmoid激活函数
        a[l] = sigmoid(a[l-1] @ Theta[l-1].T)  # 形状为(5000, 25)
        a[l] = np.c_[np.ones(m), a[l]]         # 添加偏置项，形状变为(5000, 26)
    
    # 输出层计算：最后一层不需要再添加偏置项
    a[L] = sigmoid(a[L-1] @ Theta[L-1].T)  # 形状为(5000, 10)
    
    return a

def J(a_L):
    """计算损失函数值
    参数:
        a_L: 输出层的激活值，形状为(m, s[L])
    返回:
        包含正则化项的交叉熵损失
    """
    # 计算交叉熵损失（不包含正则化项）
    res = - (1 / m) * np.sum(Y * np.log(a_L) + (1-Y) * np.log(1-a_L))
    
    # 添加正则化项（注意排除偏置项对应的权重）
    for l in range(1, L):
        # Theta[l][:, 1:]表示排除每一行的第一个元素（偏置项权重）
        res += lmd / (2 * m) * np.sum(np.power(Theta[l][:, 1:], 2))
    
    return res

def backProp():
    """反向传播算法，计算各权重矩阵的梯度
    返回:
        D: 与Theta结构相同的梯度矩阵
    """
    # 初始化Delta矩阵，用于累积梯度（与Theta形状相同）
    Delta = [None] * L
    for l in range(1, L):
        Delta[l] = np.zeros((s[l+1], s[l]+1))
    
    # 初始化delta，用于存储各层的误差项
    delta = [None] * (L+1)  # delta[0]和delta[1]未使用（输入层无误差）
    
    # 步骤1：执行前向传播，获取各层激活值
    a = forwardProp()
    #J_history.append(J(a[L]))  # 记录当前迭代的损失值（当前注释掉）
    
    # 步骤2：计算各层误差并累积梯度
    for i in range(m):  # 对每个样本计算
        # 输出层误差：预测值与真实值的差
        delta[L] = a[L][i, :] - Y[i, :]  # 形状为(s[L], )
        
        # 反向计算隐藏层误差（从倒数第二层开始到第二层）
        for l in range(L-1, 1, -1):
            # 误差计算：权重矩阵转置与上层误差的乘积，再乘以sigmoid函数的导数
            # [1:]表示排除偏置项对应的误差（因为偏置项没有输入来源）
            delta[l] = ((Theta[l].T @ delta[l+1]) * (a[l][i, :] * (1 - a[l][i, :])))[1:]
        
        # 累积梯度：误差项与激活值的外积
        for l in range(1, L):
            # delta[l+1][:,np.newaxis]将行向量转为列向量，a[l][i:i+1, :]保持为行向量
            # 外积结果为(s[l+1], s[l]+1)，与Theta[l]形状一致
            Delta[l] += delta[l+1][:,np.newaxis] @ a[l][i:i+1, :]
    
    # 步骤3：计算平均梯度并添加正则化项
    D = [None] * L  # D的结构与Theta相同
    for l in range(1, L):
        # 偏置项的权重不参与正则化（第一列保持原样）
        D[l] = (1 / m) * (Delta[l] + lmd * np.c_[np.zeros(s[l+1]), Theta[l][:, 1:]])
    
    return D

# 梯度下降算法：迭代更新权重参数
for i in range(0, num_iters):
    print(i)  # 打印当前迭代次数，便于观察训练进度
    D = backProp()  # 计算梯度
    # 更新各层权重：Theta = Theta - alpha * 梯度
    for l in range(1, L):
        Theta[l] -= alpha * D[l]

# 模型预测与评估
# 执行前向传播获取输出层结果，取最大值所在索引作为预测类别
y_pred = np.argmax(forwardProp()[L], axis=1)  # 形状为(5000, )

# 生成并打印分类报告，包括精确率、召回率、F1分数等指标
print(classification_report(y, y_pred, digits=3))

5000 轮次后，最终的预测结果如下：

              precision    recall  f1-score   support

           0      0.950     0.984     0.967       500
           1      0.946     0.978     0.962       500
           2      0.944     0.914     0.929       500
           3      0.921     0.914     0.918       500
           4      0.948     0.940     0.944       500
           5      0.919     0.904     0.911       500
           6      0.949     0.964     0.956       500
           7      0.950     0.942     0.946       500
           8      0.935     0.924     0.930       500
           9      0.926     0.924     0.925       500

    accuracy                          0.939      5000
   macro avg      0.939     0.939     0.939      5000
weighted avg      0.939     0.939     0.939      5000

如果要进行梯度检验，只需要将下述函数插入到 D = backProp() 后即可，正式训练时删去：

# Gradient Check
def gradCheck():
	for l in range(1, L):
		for i in range(s[l+1]):
			for j in range(s[l]+1):
				Theta[l][i, j] -= 0.0001
				J1 = J(forwardProp()[L])
				Theta[l][i, j] += 0.0002
				J2 = J(forwardProp()[L])
				Theta[l][i, j] -= 0.0001
				print(D[l][i, j], (J2 - J1) / 0.0002)