吴恩达机器学习笔记（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 个单元（独热编码）。

先进行预处理，包括数据预处理、随机初始化参数：

import numpy as np
import scipy.io as scio
from sklearn.metrics import classification_report# load data
data = scio.loadmat('ex4data1.mat')
X = data['X'] 			# (5000, 400), 这里不转置
y = data['y'].flatten()	# (5000, )
y[y==10] = 0# parameter
(m, n) = (5000, 401)
L = 3 				 # layer of network
s = [0, 400, 25, 10] # size of each layer
lmd = 1				 # for regularization
alpha = 0.1
num_iters = 5000
#J_history = []# One-hot encoder
Y = np.zeros((m, s[L]))
for i in range(m):Y[i][y[i]] = 1# random init Theta, (25, 401) (10, 26)
Theta = [None] * L
for l in range(1, L):Theta[l] = np.random.rand(s[l+1], s[l]+1)eps = np.sqrt(6) / np.sqrt(s[l+1] + s[l])Theta[l] = Theta[l] * 2 * eps - eps

再实现前向传播与代价函数计算，代价函数这里用于绘图、梯度检验：

def sigmoid(z):return 1 / (1 + np.exp(-z))def forwardProp():a = [None] * (L+1)# input layera[1] = np.c_[np.ones(m), X] 				# (5000, 401)# hidden layerfor l in range(2, L):a[l] = sigmoid(a[l-1] @ Theta[l-1].T)	# (5000, 25)a[l] = np.c_[np.ones(m), a[l]]			# (5000, 26)# output layera[L] = sigmoid(a[L-1] @ Theta[L-1].T)		# (5000, 10)return adef J(a_L):'''a_L : output layer animation, (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):res += lmd / (2 * m) * np.sum(np.power(Theta[l][:, 1:], 2))return res

实现反向传播算法，注意这里偏置单元的 $delta$ 误差没有意义，算出的值要忽略：

def backProp():# init Delta -> 0, should be (s_{l+1}, s_l + 1)Delta = [None] * Lfor l in range(1, L):Delta[l] = np.zeros((s[l+1], s[l]+1))# init delta, should be (s_l, )delta = [None] * (L+1)# step 1a = forwardProp()#J_history.append(J(a[L]))# step 2for i in range(m):delta[L] = a[L][i, :] - Y[i, :]  # 降维成 (s_L, )for l in range(L-1, 1, -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] += delta[l+1][:,np.newaxis] @ a[l][i:i+1, :]  # 升维成 (s_l, 1)# step 3D = [None] * L  # init D, same as Deltafor l in range(1, L):D[l] = (1 / m) * (Delta[l] + lmd * np.c_[np.zeros(s[l+1]), Theta[l][:, 1:]])return D

实现梯度下降算法并检验结果，大约耗时 20 分钟：

# Gradient Descent
for i in range(0, num_iters):print(i)D = backProp()for l in range(1, L):Theta[l] -= alpha * D[l]# Prediction，调库生成混淆矩阵
y_pred = np.argmax(forwardProp()[L], axis=1)	# (5000, )
print(classification_report(y, y_pred, digits=3))

如果要进行梯度检验，只需要将下述函数插入到 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.0001J1 = J(forwardProp()[L])Theta[l][i, j] += 0.0002J2 = J(forwardProp()[L])Theta[l][i, j] -= 0.0001print(D[l][i, j], (J2 - J1) / 0.0002)

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

              precision    recall  f1-score   support0      0.950     0.984     0.967       5001      0.946     0.978     0.962       5002      0.944     0.914     0.929       5003      0.921     0.914     0.918       5004      0.948     0.940     0.944       5005      0.919     0.904     0.911       5006      0.949     0.964     0.956       5007      0.950     0.942     0.946       5008      0.935     0.924     0.930       5009      0.926     0.924     0.925       500accuracy                          0.939      5000macro avg      0.939     0.939     0.939      5000
weighted avg      0.939     0.939     0.939      5000