ML学习笔记 #07 神经网络：反向传播

本文最后更新于：2023年1月5日 凌晨

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

  • \(x_i\) 表示输入层的第 \(i\) 个输入特征，其中 \(x_0\) 偏置项省略；
  • \(a^{(j)}_i\) 表示第 \(j\) 层的第 \(i\) 个激活单元，其中 \(a_0^{(j)}\) 偏置单元省略；
  • \(s_j\) 表示第 \(j\) 层的神经元数量（不包含偏置项），\(\hat 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\geqslant 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()=-_{i=1}m{}+{j=1}^n{{j}{2}}\

$$

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

  1. 我们使用了 \(\text{sigmoid}\) 函数作为激活函数引入非线性，其单增 \(S\) 形曲线使其在 \(y=0\) 和 \(1\) 附近导数较小；
  2. 通常我们会随机初始化参数后进行梯度下降，如果使用均方误差代价函数，对 \(\theta\) 求偏导后含有 \(\text{sigmoid}\) 的导数项，会使前期的梯度下降十分缓慢；
  3. 如果使用交叉熵代价函数，则求偏导后不含 \(\text{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)})_k \right) +(1-y_{k}^{(i)})\ln \left( 1-h_{\Theta}(x^{(i)})_k \right) \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)}\)）.

反向传播算法 | Back Propagation

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

为什么需要「反向传播」？我们知道梯度反映的是 \(J\) 在 \(\Theta_{ji}^{(l)}\) 方向的变化率，即：当参数发生微小增量时，代价函数的增量。对于前几层网络，参数的增量需要经过多层全连接传播，才会影响到输出（代价函数）。

如果我们从输入层开始计算每个参数的梯度，则每个参数都要依赖其后续路径的传播，每个参数的计算都要先算一遍后续路径的梯度。而如果从输出层开始逐层向前计算并存储了梯度，则解耦了依赖关系，前一层可以利用后一层的计算结果。有点类似「动态规划」中「递归求解子问题」的思想。

数学推导：问题转化

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

我们首先对第 \(l+1\) 层的第 \(j\) 个神经元（偏置单元除外）定义一个 \(\text{delta}\) 误差： \[ \delta_j^{(l+1)}=\frac{\partial J}{\partial z_j^{(l+1)}} \] 其中 \(z_j^{(l+1)}=\sum\limits_{k=0}^{s_{l}}\Theta^{(l)}_{jk}a^{(l)}_{k}\)，即 \(a_j^{(l+1)}\) 神经元取 \(\text{sigmoid}\) 激活前的输出值。换句话说，\(a^{(l+1)}_j=g\left(z_j^{(l+1)}\right)\)。

于是根据链式求导法则，我们有： \[ \frac{\partial J}{\partial \Theta^{(l)}_{ji}}=\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)} \] 写作矩阵形式： \[ \boxed{\Delta ^{\left( l \right)}=\delta ^{\left( l+1 \right)}\cdot \left( a^{\left( l \right)} \right) ^T}\tag{1} \] 显然，一旦求出当前层的 \(\text{delta}\) 误差，那么传入当前层的各参数的梯度也可求出。所以现在问题转化为求解 \(\delta^{(l+1)}\)。

数学推导：输出层误差

既然要递归求解子问题，首先就是要算输出层，即 \(\delta_j^{(L)}\)： \[ \begin{aligned} \delta _{j}^{\left( L \right)}&=\frac{\partial J}{\partial z_{j}^{\left( L \right)}}=\frac{\partial J}{\partial a_{j}^{\left( L \right)}}\cdot \frac{\partial a_{j}^{\left( L \right)}}{\partial z_{j}^{\left( L \right)}}\\ &=\frac{\partial J}{\partial a_{j}^{\left( L \right)}}\cdot g'\left( z_{j}^{\left( L \right)} \right)\\ &=-\left( \frac{y_j}{a_{j}^{\left( L \right)}}-\frac{1-y_j}{1-a_{j}^{\left( L \right)}} \right) \cdot \left( a_{j}^{\left( L \right)}\left( 1-a_{j}^{\left( L \right)} \right) \right)\\ &=a_{j}^{\left( L \right)}-y_j\\ \end{aligned} \]

写作矩阵形式： \[ \boxed{\delta^{(L)}=a^{(L)}-y}\tag{2} \]

此处推导过程的一些注释：

关于第二行，请注意：$a^{( L )}=g( z^{( L )} ) $；

关于第三行，请注意 \(\text{sigmoid}\) 函数的性质：\(g'(z)=g(z)(1-g(z))\)；

以及第三行偏导项的计算，请注意 \(a^{(L)}=h_\Theta(x)\)，所以 \(J(\Theta)\) 在现在的假设条件下可以写作： \[ J(\Theta)=-\left[\sum_{k=1}^Ky_k\ln(a^{(L)}_k)+(1-y_k)\ln(1-a_k^{(L)})\right] \]

数学推导：隐藏层误差

下面计算第 \(l\) 层（\(2\leqslant l<L\)）的 \(\delta_j^{(l)}\)： \[ \begin{aligned} \delta _{j}^{\left( l \right)}&=\frac{\partial J}{\partial z_{j}^{\left( l \right)}}\\ &=\frac{\partial J}{\partial z^{\left( l+1 \right)}}\cdot \frac{\partial z^{\left( l+1 \right)}}{\partial a_{j}^{\left( l \right)}}\cdot \frac{\partial a_{j}^{\left( l \right)}}{\partial z_{j}^{\left( l \right)}}\\ &={\delta ^{\left( l+1 \right)}}^T\cdot \Theta _{\cdot ,j}^{\left( l \right)}\cdot g'\left( z_{j}^{\left( l \right)} \right)\\ &={\delta ^{\left( l+1 \right)}}^T\cdot \Theta _{\cdot ,j}^{\left( l \right)}\cdot \left( a_{j}^{\left( l \right)} \left( 1-a_{j}^{\left( l \right)} \right) \right)\\ \end{aligned} \] 观察式子的前半部分，可以看出下一层所有结点的 \(\text{delta}\) 误差均会影响当前层的任一结点，而影响则是通过该结点链出的权重参数反向传播。

转置后，写作矩阵形式： \[ \boxed{\delta ^{\left( l \right)}=\left( {\Theta ^{\left( l \right)}}^T\delta ^{\left( l+1 \right)} \right) \odot \left( a^{\left( l \right)}\odot \left( 1-a^{\left( l \right)} \right) \right) }\tag{3} \] 其中 \(\odot\) 表示 Hadmard 积，即两个向量对应位置相乘。

步骤总结

以上是对一组数据的推导，我们得到了三个重要的结果 \((1)(2)(3)\)： \[ \boxed{ \begin{aligned} \Delta ^{\left( l \right)}&=\delta ^{\left( l+1 \right)}\cdot \left( a^{\left( l \right)} \right) ^T&&1\leqslant l<L\\ \delta ^{\left( L \right)}&=a^{\left( L \right)}-y\\ \delta ^{\left( l \right)}&=\left( {\Theta ^{\left( l \right)}}^T\delta ^{\left( l+1 \right)} \right) \odot \left( a^{\left( l \right)}\odot \left( 1-a^{\left( l \right)} \right) \right) &&2\leqslant l<L\\ \end{aligned} } \] 而 \(m\) 组数据只需要在一些地方进行累加即可。设数据集为 \(\left\{ \left( x^{\left( i \right)},y^{\left( i \right)} \right) \mid 1\leqslant i\leqslant m \right\}\)，则反向传播算法的步骤如下：

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

  2. 遍历数据集，设当前数据为 \((x^{(i)},y^{(i)})\)：

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

  3. 计算 \(D\) 矩阵，对偏置项以外的参数进行正则化： \[ D_{ij}^{(l)}:=\begin{cases} \frac{1}{m}\left( \Delta _{ij}^{(l)}+\lambda \Theta _{ij}^{(l)} \right)& \mathrm{if}\; j\ne 0\\ \frac{1}{m}\Delta _{ij}^{(l)}& \mathrm{if}\; j=0\\ \end{cases} \] 这就是最终的梯度矩阵：\(\begin{aligned}\frac{\partial J}{\partial \Theta_{ij}^{(l)}}=D^{(l)}_{ij}\end{aligned}\)。

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

代码实现

参数展开 | Unrolling Parameters

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

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

梯度检验 | Gradient Checking

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

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

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

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

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

参数随机初始化

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

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

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

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

练习

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def sigmoid(z):
	return 1 / (1 + np.exp(-z))

def forwardProp():
	a = [None] * (L+1)
	# input layer
	a[1] = np.c_[np.ones(m), X] 				# (5000, 401)
	# hidden layer
	for 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 layer
	a[L] = sigmoid(a[L-1] @ Theta[L-1].T)		# (5000, 10)
	return a

def 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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def backProp():
	# init Delta -> 0, should be (s_{l+1}, s_l + 1)
	Delta = [None] * L
	for 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 1
	a = forwardProp()
	#J_history.append(J(a[L]))
	# step 2
	for 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 3
	D = [None] * L  # init D, same as Delta
	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

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

1
2
3
4
5
6
7
8
9
10
# 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() 后即可，正式训练时删去：

1
2
3
4
5
6
7
8
9
10
11
# 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)

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
              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
学习笔记 > 机器学习
#ML
ML学习笔记 #07 神经网络：反向传播
https://hwcoder.top/ML-Note-7
作者
Wei He
发布于
2022年2月6日
许可协议