在《神经网络的梯度推导与代码验证》之CNN的前向传播和反向梯度推导 中，我们学习了CNN的前向传播和反向梯度求导，但知识仍停留在纸面。本篇章将基于深度学习框架tensorflow验证我们所得结论的准确性，以便将抽象的数学符号和实际数据结合起来，将知识固化。更多相关内容请见《神经网络的梯度推导与代码验证》系列介绍。
需要用到的库有tensorflow和numpy，其中tensorflow其实版本>=2.0.0就行
import tensorflow as tf
import numpy as np
np.random.seed(0)
然后是定义交叉熵损失函数：
def get_crossentropy(y_pred, y_true):
return -tf.reduce_sum(y_true * tf.math.log(y_pred))
接下来正式进入代码主体：
with tf.GradientTape(persistent=True) as t:
# -------input-----------
x = tf.constant(np.random.randn(1, 9, 9, 1).astype(np.float32))
y_true = np.array([[0.3, 0.5, 0.2]]).astype(np.float32)
t.watch(x)
# -----conv2d l1----------
l1 = tf.keras.layers.Conv2D(filters=1, kernel_size=(3, 3), strides=2) z_l1 = l1(x)
t.watch([z_l1])
a_l1 = tf.nn.relu(z_l1)
t.watch([a_l1])
# -----max pooling--------
l2 = tf.keras.layers.MaxPool2D(pool_size=(2, 2)) z_l2 = l2(a_l1)
t.watch([z_l2])
a_l2 = tf.keras.layers.Flatten()(z_l2)
t.watch([a_l2])
# ---------FNN------------
l3 = tf.keras.layers.Dense(3) z_l3 = l3(a_l2)
t.watch([z_l3])
a_l3 = tf.math.softmax(z_l3)
t.watch([a_l3])
# ---------loss----------
loss = get_crossentropy(y_pred=a_l3, y_true=y_true)
上面这段代码实现的是输入--卷积层--池化层--全连接层--输出的前向过程。因为是咱们的目的是做验证，所以没必要把网络加深，让这个麻雀五脏俱全即可。
input部分，输入x是形状为(1, 9, 9, 1)的张量，可以理解为一张单通道的9*9的图。标签y_true是3维的概率向量。在上面的代码中我们只考虑一条（x，y_true）样本。
--------前向传播验证---------
我们先验证下第7+9行代码是否在对x做stride=2的（默认valid模式）卷积操作。
先瞧一瞧这一层的卷积核和偏置：
np.squeeze(l1.kernel.numpy())
Out[4]:
array([[-0.3447126 , -0.23770776, -0.20545131],
[ 0.40415084, -0.56749415, 0.13746339],
[-0.5106965 , -0.36734173, -0.18053415]], dtype=float32)
l1.bias.numpy()
Out[8]: array([0.], dtype=float32)
然后再看看输入x和输出z_l1：
np.squeeze(x)
Out[6]:
array([[ 1.7640524 , 0.4001572 , 0.978738 , 2.2408931 , 1.867558 ,
-0.9772779 , 0.95008844, -0.1513572 , -0.10321885],
[ 0.41059852, 0.14404356, 1.4542735 , 0.7610377 , 0.12167501,
0.44386324, 0.33367434, 1.4940791 , -0.20515826],
[ 0.3130677 , -0.85409576, -2.5529897 , 0.6536186 , 0.8644362 ,
-0.742165 , 2.2697546 , -1.4543657 , 0.04575852],
[-0.18718386, 1.5327792 , 1.4693588 , 0.15494743, 0.37816253,
-0.88778573, -1.9807965 , -0.34791216, 0.15634897],
[ 1.2302907 , 1.2023798 , -0.3873268 , -0.30230275, -1.048553 ,
-1.420018 , -1.7062702 , 1.9507754 , -0.5096522 ],
[-0.4380743 , -1.2527953 , 0.7774904 , -1.6138978 , -0.21274029,
-0.89546657, 0.3869025 , -0.51080513, -1.1806322 ],
[-0.02818223, 0.42833188, 0.06651722, 0.3024719 , -0.6343221 ,
-0.36274117, -0.67246044, -0.35955316, -0.8131463 ],
[-1.7262826 , 0.17742614, -0.40178093, -1.6301984 , 0.46278226,
-0.9072984 , 0.0519454 , 0.7290906 , 0.12898292],
[ 1.1394007 , -1.2348258 , 0.40234163, -0.6848101 , -0.87079716,
-0.5788497 , -0.31155252, 0.05616534, -1.1651498 ]],
dtype=float32)
np.squeeze(z_l1)
Out[7]:
array([[-0.00542112, -0.17352474, -1.3421125 , -1.6447177 ],
[-1.1239526 , 1.6031268 , 1.1616374 , -0.78091574],
[-0.14451274, 1.5910958 , 2.1035302 , 1.1354219 ],
[-1.1602874 , 1.0651501 , 1.8656987 , 0.4581319 ]],
dtype=float32)
为提高效率，我们就验证一丢丢就好：
np.sum(np.squeeze(x)[0:3, 0:3] * np.squeeze(l1.kernel.numpy()))
Out[8]: -0.0054210722
np.sum(np.squeeze(x)[2:5, 0:3] * np.squeeze(l1.kernel.numpy()))
Out[9]: -1.1239524
np.sum(np.squeeze(x)[0:3, 2:5] * np.squeeze(l1.kernel.numpy()))
Out[10]: -0.17352472
对比红色标记的部分，发现我们定义的卷积层确实有乖乖做stride=2的valid卷积操作。
接下来再看看pool_size=(2, 2)池化层的操作：
np.squeeze(a_l1)
Out[11]:
array([[0. , 0. , 0. , 0. ],
[0. , 1.6031268, 1.1616374, 0. ],
[0. , 1.5910958, 2.1035302, 1.1354219],
[0. , 1.0651501, 1.8656987, 0.4581319]], dtype=float32)
np.squeeze(z_l2)
Out[12]:
array([[1.6031268, 1.1616374],
[1.5910958, 2.1035302]], dtype=float32)
也没有什么问题。
接下来的FNN的前向传播验证可参考FNN（DNN）前向和反向传播过程的代码验证，这里就略过了。
----------反向梯度计算的验证----------
老样子，我们从尾向头开始验证
先看$\frac{\partial l}{\partial\boldsymbol{z}\_ l3}$：
# -----dl_dz3------
dl_dz3 = t.gradient(loss, z_l3)
my_dl_dz3 = a_l3 - y_true
dl_dz3
Out[13]: <tf.Tensor: shape=(1, 3), dtype=float32, numpy=array([[ 0.20361423, -0.4315467 , 0.22793245]], dtype=float32)>
my_dl_dz3
Out[14]: <tf.Tensor: shape=(1, 3), dtype=float32, numpy=array([[ 0.20361423, -0.4315467 , 0.22793244]], dtype=float32)>
上面dl_dz3 = t.gradient(loss, z_l3)表示用tensorflow微分工具求得的$\frac{\partial l}{\partial\boldsymbol{z}\_ l3}$；而带my_前缀的则是根据先前推导得到的结论手动计算出来的结果。后续的符号同样沿用这样的命名规则。
$\frac{\partial l}{\partial\boldsymbol{z}\_ l3}$应当满足：
$\frac{\partial l}{\partial\boldsymbol{z}\_ l3} = \boldsymbol{a}\_ l3 - \boldsymbol{y}_{true}$
推导过程参考《神经网络的梯度推导与代码验证》之FNN（DNN）的前向传播和反向梯度推导 的2.2小节。
从结果看来并没有什么问题。
我们直接跳过FNN层的反向梯度验证（可参考FNN的代码验证），来到max pooling层：
（跳过FNN的反向梯度验证意味着此时我们是已经能够计算得到$\frac{\partial l}{\partial\boldsymbol{z}\_ l2}$）
于是我们直接验证$\frac{\partial l}{\partial\boldsymbol{z}\_ l1}$
下面这是tensorflow自动微分给出的$\frac{\partial l}{\partial\boldsymbol{z}\_ l1}$结果
inverse_mp = np.squeeze(t.gradient(loss, z_l1)) inverse_mp
Out[15]:
array([[ 0. , 0. , 0. , 0. ],
[ 0. , 0.08187968, -0.07518297, 0. ],
[ 0. , -0.2259186 , 0.5417712 , 0. ],
[ 0. , 0. , 0. , 0. ]],
dtype=float32)
接下来我们手动实现下CNN的前向传播和反向梯度推导中关于max pooling层的结论：
$\boldsymbol{\delta}_{k}^{l - 1} = upsample\left( \boldsymbol{\delta}_{k}^{l} \right) \odot \sigma^{'}\left( \boldsymbol{z}_{k}^{l - 1} \right)$
我们先看看卷积层出来的结果z_l1，因为z_l1是个高维张量不太好直接观察，所以用squeeze去掉了一些多余的axis，于是batch_size和channel这两个为1的维度被我去掉，但并不影响验证。
flat_z_l1 = np.squeeze(z_l1) flat_z_l1
Out[16]:
array([[-0.00542112, -0.17352474, -1.3421125 , -1.6447177 ],
[-1.1239526 , 1.6031268 , 1.1616374 , -0.78091574],
[-0.14451274, 1.5910958 , 2.1035302 , 1.1354219 ],
[-1.1602874 , 1.0651501 , 1.8656987 , 0.4581319 ]],
dtype=float32)
# z_l1经过激活函数relu之后变成了下面的a_l1
flat_a_l1 = np.squeeze(a_l1) flat_a_l1
Out[17]:
array([[0. , 0. , 0. , 0. ],
[0. , 1.6031268, 1.1616374, 0. ],
[0. , 1.5910958, 2.1035302, 1.1354219],
[0. , 1.0651501, 1.8656987, 0.4581319]], dtype=float32)
# a_l1池化后得到下面的结果，我们记住池化前后的元素位置用于后面up sampling
flat_z_l2 = np.squeeze(z_l2) flat_z_l2
Out[18]:
array([[1.6031268, 1.1616374],
[1.5910958, 2.1035302]], dtype=float32)
因为前面假设我们已知$\frac{\partial l}{\partial\boldsymbol{z}\_ l2}$，所以这里直接用tensorflow的自动微分工具给出结果：
# 下面是dl_dz2的结果
dl_dz2 = np.squeeze(t.gradient(loss, z_l2)) dl_dz2
Out[23]:
array([[ 0.08187968, -0.07518297],
[-0.2259186 , 0.5417712 ]], dtype=float32)
根据$\frac{\partial l}{\partial\boldsymbol{z}\_ l1} = upsample\left( \frac{\partial l}{\partial\boldsymbol{z}\_ l2} \right) \odot \sigma^{'}\left( {\boldsymbol{z}\_ l1} \right)$我们简单脑补一下下面这个东西就好，没有必要写代码实现这个upsampling（upsampling时，四个元素的位置要参考max pooling时的位置记录）。
按照公式给出的规则我们脑补出来的结果确实跟下面这个自动微分的结果是一致的。
dl_dz1 = t.gradient(loss, z_l1) np.squeeze(dl_dz1)
Out[26]:
array([[ 0. , 0. , 0. , 0. ],
[ 0. , 0.08187968, -0.07518297, 0. ],
[ 0. , -0.2259186 , 0.5417712 , 0. ],
[ 0. , 0. , 0. , 0. ]],
dtype=float32)
最后是略为麻烦的，验证递推公式：$\boldsymbol{\delta}^{l - 1} = \left( \frac{\partial\boldsymbol{z}^{l}}{\partial\boldsymbol{z}^{l - 1}} \right)^{T}\boldsymbol{\delta}^{l} = \boldsymbol{\delta}^{l}*rot180\left( \boldsymbol{W}^{l} \right) \odot \sigma^{'}\left( \boldsymbol{z}^{l - 1} \right)$， 其中$\boldsymbol{\delta}^{l} = \frac{\partial l}{\partial\boldsymbol{z}^{\boldsymbol{l}}}$。
在我们这个“麻雀”小例子中，我们只有一层CNN层，也就意味着，这里loss对z_l1 = l1(x) 的导数就是上述公式中的$\boldsymbol{\delta}^{l}$，而loss对x的导数则对映着公式中的$\boldsymbol{\delta}^{l-1}$，而且因为x位于输入层，所以是没有激活函数的，或者说激活函数是线性函数，即$\sigma\left( \boldsymbol{z}^{l - 1} \right)=\boldsymbol{z}^{l - 1}$
于是我们最终要验证的东西简化成了：
$\frac{\partial l}{\partial\boldsymbol{x}} = \frac{\partial l}{\partial\boldsymbol{z}\_ l1}*rot180\left( \boldsymbol{W}^{1} \right)$
虽说一般情况下，我们不会去求loss对输入的导数就是了，因为通常我们训练的对象是神经网络而非样本数据。但在生成问题下，例如GAN，这种玩法还是很常见的。
回到正题，我们现在已经明确要对比什么了。接下来开始码代码吧。。
现已知loss对z_l1的导数
# 已知dl_dz1
dl_dz1 = t.gradient(loss, z_l1)
我们要通过它来求loss对输入x的导数，也就是下面这个东西，这是tensorflow给出的答案，接下来我们要通过递推公式得到一个跟这一样的答案。
# -----dl_dx--------
dl_dx = np.squeeze(t.gradient(loss, x)) np.squeeze(dl_dx)
Out[]:
array([[ . , . , . , . , -0.02591492,
0.0461329 , -0.02298165, . , . ],
[ . , . , . , . , 0.00783426,
-0.0234191 , 0.00530995, . , . ],
[ . , . , . , . , 0.05378128,
-0.02789585, 0.07432017, . , . ],
[ . , . , . , . , . ,
. , . , . , . ],
[ . , . , . , . , . ,
. , 0.02716039, -0.04835004, 0.02408614],
[ . , . , . , . , . ,
. , -0.00821077, 0.02454463, -0.00556515],
[-0.02953471, 0.05257672, -0.02619171, . , . ,
. , -0.05636601, 0.02923652, -0.07789201],
[ 0.00892854, -0.02669027, 0.00605165, . , . ,
. , . , . , . ],
[ 0.06129343, -0.03179233, 0.08470119, . , . ,
. , . , . , . ]],
dtype=float32)
我们注意到 l1 = tf.keras.layers.Conv2D(filters=1, kernel_size=(3, 3), strides=2)，这是一个stride>1的情况。让我们回到CNN的前向传播和反向梯度推导的3.2.2中的那个stride>1的例子。当stride>1时，$\frac{\partial l}{\partial\boldsymbol{z}\_\boldsymbol{l}1}$矩阵的元素是要散开成一个较大的矩阵（空的位置补0，类似upsampling的感觉）然后跟旋转了180度的卷积核做valid模式的卷积就能完成跨越CNN层的梯度递推。
下面是关键点回顾：
- 前向计算有：
$conv2D\left( {\boldsymbol{a}^{l - 1},\boldsymbol{W}^{l},^{'}valid^{'},~stride = 2} \right) = \left\lbrack \begin{array}{lll} z_{11} & 0 & z_{12} \\ 0 & 0 & 0 \\ z_{21} & 0 & z_{22} \\ \end{array} \right\rbrack\overset{down~sampling}{\rightarrow}\boldsymbol{z}^{l} = \left\lbrack \begin{array}{ll} z_{11} & z_{12} \\ z_{21} & z_{22} \\ \end{array} \right\rbrack$
- 计算反向梯度时有：
loss $l$对矩阵$\boldsymbol{z}^{l}$的导数，即$\boldsymbol{\delta}^{l}$，它经过上采用后再跟$rot180\left( \boldsymbol{W}^{l} \right)$进行stride=1的full模式的卷积运算的结果就是上面例子2的最终结果，即：
$\boldsymbol{\delta}^{l} = \left\lbrack \begin{array}{ll} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \\ \end{array} \right\rbrack\overset{up~sampling}{\rightarrow}\left\lbrack \begin{array}{lll} \delta_{11} & 0 & \delta_{12} \\ 0 & 0 & 0 \\ \delta_{21} & 0 & \delta_{22} \\ \end{array} \right\rbrack$
$\frac{\partial l}{\partial\boldsymbol{a}^{l - 1}} = conv2D\left( {\left\lbrack \begin{array}{lll} \delta_{11} & 0 & \delta_{12} \\ 0 & 0 & 0 \\ \delta_{21} & 0 & \delta_{22} \\ \end{array} \right\rbrack,rot180\left( \boldsymbol{W}^{l} \right),'full'} \right)$
但是，就算是照着上面这个模式做，自己写起代码来也是有点麻烦的，而tensorflow恰好提供了用于完成上面这种操作的layer，Conv2DTranspose。
inverse_conv2d = tf.keras.layers.Conv2DTranspose(filters=, kernel_size=(, ), strides=)
inverse_conv2d.build(dl_dz1.shape)
上面两行代码完成了设置和初始化Conv2DTranspose层，它做的事情正好就是：先对输入做扩张比率为strides的upsampling，然后将它和旋转了180度的自己的卷积核做full模式的卷积操作。
但目前Conv2DTranspose还不能用，因为前传和反传期间卷积核只是被旋转，其元素的值是不变的，因此初始化的Conv2DTranspose层的卷积核和偏置还需要被赋值成Conv2D时所用的那些：
inverse_conv2d.kernel = l1.kernel
inverse_conv2d.bias = l1.bias
然后就是见证奇迹的时刻了：
my_dl_dx = np.squeeze(inverse_conv2d(dl_dz1)) my_dl_dx
Out[]:
array([[ . , . , . , . , -0.02591492,
0.0461329 , -0.02298165, . , . ],
[ . , . , . , . , 0.00783426,
-0.0234191 , 0.00530995, . , . ],
[ . , . , . , . , 0.05378128,
-0.02789585, 0.07432017, . , . ],
[ . , . , . , . , . ,
. , . , . , . ],
[ . , . , . , . , . ,
. , 0.02716039, -0.04835004, 0.02408614],
[ . , . , . , . , . ,
. , -0.00821077, 0.02454463, -0.00556515],
[-0.02953471, 0.05257672, -0.02619171, . , . ,
. , -0.05636601, 0.02923652, -0.07789201],
[ 0.00892854, -0.02669027, 0.00605165, . , . ,
. , . , . , . ],
[ 0.06129343, -0.03179233, 0.08470119, . , . ,
. , . , . , . ]],
dtype=float32)
跟前面自动微分工具给出的结果一致。
最后是loss对$W$，$b$的导数的验证。
根据结论，loss对$W$的导数满足：
$\frac{\partial l}{\partial\boldsymbol{W}^{l}} = {\sum\limits_{k}{\sum\limits_{l}\delta_{k,l}^{l}}}\sigma\left( z_{k + x,l + y}^{l - 1} \right) = conv2D\left( {\sigma\left( \boldsymbol{z}^{l - 1} \right),\boldsymbol{\delta}^{l}~,'valid'} \right)$
需要注意的是，对于stride>1的情况，如果要求$\frac{\partial\boldsymbol{l}}{\partial\boldsymbol{W}^{\boldsymbol{l}}}$和$\frac{\partial\boldsymbol{l}}{\partial\boldsymbol{b}^{\boldsymbol{l}}}$，我们需要将$\boldsymbol{\delta}^{\boldsymbol{l}}$做扩展比率=strides上采样先。
那我们先看看原本$\boldsymbol{\delta}^{\boldsymbol{l}}$的样子，为方便查看依然用np.squeeze去掉多余的axis。
np.squeeze(dl_dz1)
Out[]:
array([[ . , . , 0.13701844, . ],
[ . , . , . , . ],
[ . , . , . , -0.14360355],
[ 0.15615712, . , . , . ]],
dtype=float32)
接下来我们手动给它做up sampling：
new_kernel = np.squeeze(dl_dz1)
col = np.zeros()
new_kernel = np.column_stack((new_kernel[:, ], col, new_kernel[:, ], col, new_kernel[:, ], col, new_kernel[:, ]))
row = np.zeros()
new_kernel = np.row_stack((new_kernel[, :], row, new_kernel[, :], row, new_kernel[, :], row, new_kernel[, :])) np.squeeze(new_kernel)
Out[]:
array([[ . , . , . , . , 0.13701844,
. , . ],
[ . , . , . , . , . ,
. , . ],
[ . , . , . , . , . ,
. , . ],
[ . , . , . , . , . ,
. , . ],
[ . , . , . , . , . ,
. , -0.14360355],
[ . , . , . , . , . ,
. , . ],
[ 0.15615712, . , . , . , . ,
. , . ]], dtype=float32)
从原来4*4扩展成了7*7，即4个元素间，放了3个0。
完成了对$\boldsymbol{\delta}^{\boldsymbol{l}}$的上采样后接下来就是计算$\sigma\left( \boldsymbol{z}^{l - 1} \right)$，它是卷积层的输入，在本例中就是输入x，所以$\sigma\left( \boldsymbol{z}^{l - 1} \right)=x$。
剩下的工作就是做卷积的。既然有卷积层那就不用我们自己动手写卷积操作了。我们定义好CNN层，然让将其卷积核替换成我们事先准备好的那个上采样后的$\boldsymbol{\delta}^{\boldsymbol{l}}$。
conv2d = tf.keras.layers.Conv2D(filters=, kernel_size=new_kernel.shape)
conv2d.build(inverse_conv2d(dl_dz1).shape) new_kernel = new_kernel[:, :, np.newaxis, np.newaxis].astype(np.float32)
conv2d.kernel = new_kernel
再就是对比结果了：
dl_dW = np.squeeze(t.gradient(loss, l1.kernel)) dl_dW
Out[]:
array([[-0.10748424, -0.00609609, 0.364539 ],
[ 0.1810095 , -0.13556898, -0.19335817],
[-0.135235 , -0.18566257, -0.23072638]], dtype=float32)
my_dl_dW = np.squeeze(conv2d(x)) my_dl_dW
Out[]:
array([[-0.10748424, -0.00609609, 0.364539 ],
[ 0.1810095 , -0.13556898, -0.19335817],
[-0.135235 , -0.18566257, -0.23072638]], dtype=float32)
验证无误。
剩下就是$\frac{\partial\boldsymbol{l}}{\partial\boldsymbol{b}^{\boldsymbol{l}}}$了。
根据公式，它满足$\frac{\partial l}{\partial b^{l}} = \frac{\partial l}{\partial b_{x,y}^{l}} = {\sum\limits_{k}{\sum\limits_{l}\delta_{k,l}^{l}}}$，这个简单，就只是对$\boldsymbol{\delta}^{\boldsymbol{l}}$的所有元素做累加和：
dl_db = np.squeeze(t.gradient(loss, l1.bias))
dl_db
Out[]: array(0.15956555, dtype=float32) my_dl_db = np.sum(new_kernel.astype(np.float32))
my_dl_db
Out[]: 0.15956555
也没有问题。
（欢迎转载，转载请注明出处。欢迎留言或沟通交流： lxwalyw@gmail.com）