动手实现CNN
本文将使用 Python 从头到尾实现一个卷积层。当然,既然是自己动手实现,肯定不能用到 Pytorch、TensorFlow 这些现成的库。本文中仅会用到 numpy。文末有 github 仓库地址方便读者测试,仓库中会测试我们自己实现的 CNN 在 MNIST 手写数字识别上的准确率。
卷积的基本计算
先来回顾一下卷积基本计算,分别计算输出特征图的尺寸、卷积的计算量、卷积的参数量。这一小节就直接给出公式并简单解释一下,详细的计算细节和代码统计参数量/计算量的方法可参考:图像卷积及其计算(特征图尺寸、参数量、计算量)
记输入特征图尺寸 H i n × W i n × C i n H_{in}\times W_{in}\times C_{in} Hin×Win×Cin ,卷积核尺寸为 K × K K\times K K×K ,卷积层填充为 P P P,步长为 S S S, 输出特征图通道数为: C o u t C_{out} Cout 。
输出特征图尺寸
按照卷积的滑窗方法计算即可。
H
o
u
t
=
(
H
i
n
+
2
P
−
K
)
/
/
2
+
1
W
o
u
t
=
(
W
i
n
+
2
P
−
K
)
/
/
2
+
1
H_{out}=(H_{in}+2P-K)\ //\ 2\ +\ 1\\ W_{out}=(W_{in}+2P-K)\ //\ 2\ +\ 1
Hout=(Hin+2P−K) // 2 + 1Wout=(Win+2P−K) // 2 + 1
参数量
每个卷积核的参数量是
C
i
n
×
K
×
K
C_{in}\times K\times K
Cin×K×K ,共
C
o
u
t
C_{out}
Cout 个卷积核。
C
i
n
×
C
o
u
t
×
K
×
K
C_{in}\times C_{out}\times K\times K
Cin×Cout×K×K
计算量
先算出输出特征图每个位置的计算量,然后再乘以输出特征图的尺寸即可。前面括号中就是每个位置的计算量,其中前一项是乘法计算量,后一项是加法计算量,注意这里没有考虑偏置,如果有偏置需要再加一。后面括号是特征图的尺寸。
(
C
i
n
×
K
×
K
+
C
i
n
×
K
×
K
−
1
)
×
(
C
o
u
t
×
H
o
u
t
×
W
o
u
t
)
(C_{in}\times K\times K+C_{in}\times K\times K-1)\times (C_{out}\times H_{out}\times W_{out})
(Cin×K×K+Cin×K×K−1)×(Cout×Hout×Wout)
im2col
怎样实现卷积?
在初学卷积原理时,我们知道,CV 中的卷积就是拿一个卷积核,以滑动窗口的形式去特征图上进行计算,每个位置的计算结果就是输出特征图的一个元素(如下面动图所示)。这样理解当然是没问题的,但是在实际实现时,并不会按照这种方式去滑窗、取数、计算。为什么呢?访存开销太大。我们知道,现代计算机中 CPU 的频率比内存快得多,因此,在程序中,访存的开销要比实际计算的开销大得多。保证规律的、高效的访存,能大幅提升程序的性能。
如果按照滑窗卷积进行计算,在计算机操作时,需要将其存入内存当中再操作(按照“行先序”):
这样一来,特征图y11,y12,y21,y22的计算如下所示:
按照卷积计算,实线标注出卷积计算中的访存过程(对应数据相乘),我们可以看到这一过程是非常散乱和混乱的。直接用卷积的计算方式是比较愚蠢的。
这时候就要用到im2col操作。
im2col原理
im2col 是卷积实现中的关键一步。一句话来介绍im2col操作的话,就是通过牺牲空间的手段(约扩增K×K倍),将特征图转换成庞大的矩阵来进行卷积计算。
其实思路非常简单:
把每一次循环所需要的数据都排列成列向量,然后逐一堆叠起来形成矩阵(按通道顺序在列方向上拼接矩阵)。
比如 C i n × H i n × W i n C_{in}\times H_{in}\times W_{in} Cin×Hin×Win大小的输入特征图, K × K K\times K K×K 大小的卷积核,输出大小为 C o u t × W o u t × H o u t C_{out}\times W_{out}\times H_{out} Cout×Wout×Hout,输入特征图将按需求被转换成 ( K × K × C i n ) × ( W o u t × H o u t ) (K\times K\times C_{in})\times (W_{out}\times H_{out}) (K×K×Cin)×(Wout×Hout)的矩阵,卷积核将被转换成 C o u t × ( K × K × C i n ) C_{out}\times (K\times K\times C_{in}) Cout×(K×K×Cin) 的矩阵,这样二者直接进行矩阵乘法计算就实现了卷积。
如下图中:
- 对于输入特征图:将每个将要与卷积核进行计算的滑动窗口抽出来,作为一列向量,组成 col 数组
- 对于卷积核:也是平铺开来
这样,col 数组直接与平铺开的卷积核进行矩阵乘法,就实现了卷积操作。而矩阵乘法的访存明显比滑窗卷积要规律的多,并且,矩阵乘法在各个高性能计算库(如 OpenBlas 等)中都已经经过了极致的优化,有很高的性能。
im2col实现
im2col 的实现思路非常简单:先计算出根据输入特征图尺寸(padding之后的)、卷积核的宽高和步长,先计算出输出特征图的尺寸。然后遍历输出特征图的每个位置,从输入特征图中提取出计算该位置的窗口,放到 col 数组中即可。
下面是一个 im2col 函数的 python 实现,代码量很少,去掉注释和空行一共也就十五六行。笔者已经在代码中添加了详尽的注释,供读者参考。
import numpy as np
def im2col(img, k_h, k_w, s_h, s_w):
# 输入参数
# img: 输入特征图, bs, h, w, c
# k_h, k_w: 卷积核的高, 卷积核的宽
# s_h, s_w: 高方向的步长, 宽方向的步长
# 返回值
# col: 输入特征图转换后的结果 形状为: bs * out_h * out_w, k_h * k_w * in_c
batch_size, in_h, in_w, in_c = img.shape
# 计算输出特征图的尺寸
out_h = (in_h - k_h) // s_h + 1
out_w = (in_w - k_w) // s_w + 1
# 初始化一个col数组
col = np.empty( (batch_size * out_h * out_w, k_h * k_w * in_c), dtype=np.float32)
# 批次跨度, 接下来会作为填充col时的步长
# 其值为输出特征图的尺寸, 用它作为步长, 才能正确的填充批次中每个样本的对应位置
batch_span = out_h * out_w
# 遍历输出特征图的每个位置, 填充col数组
for r in range(out_h):
r_start = r * s_h # 输入特征图中的行开始位置
col_r = r * out_h # col数组中的开始位置
for c in range(out_w):
c_start = c * s_w # 输入特征图中的列开始位置
patch = img[ :, r_start: r_start + k_h, c_start: c_start + k_w, : ] # 取到当前输出位置对应输入中计算卷积的窗口
patch = patch.reshape(batch_size, -1) # 拉直为向量
col[col_r:: batch_span, :] = patch # 填充到col数组中
return col
测试
我们在测试中使用的例子图示相同:
- 输入特征图:尺寸为 4x4,数值为 0-15;
- 卷积核宽高均为 3,步长均为 1。
if __name__ == '__main__':
data16 = [[ 0, 1, 2, 3],
[4, 5, 6, 7],
[8, 8, 10, 11],
[12, 13, 14, 15]]
x = np.array(data16).reshape(1, 4, 4, 1)
k_h, k_w = (3, 3)
s_h, s_w = (1, 1)
y = im2col(x, k_h, k_w, s_h, s_w)
print(y.shape)
以上测试代码的输出为:
[[ 0. 1. 2. 4. 5. 6. 8. 9. 10.]
[ 1. 2. 3. 5. 6. 7. 9. 10. 11.]
[ 4. 5. 6. 8. 9. 10. 12. 13. 14.]
[ 5. 6. 7. 9. 10. 11. 13. 14. 15.]]
可以与上面 im2col 的示意图进行对比,完全一致。
填充
如果卷积层设定有填充的话,需要在将特征图送入 im2col 之前先预处理,完成填充。
填充模式简介
TensorFlow 中的填充模式右 SAME
和 VALID
两种。如果填充模式设置为 SAME
,则保证输入特征图大小和输出特征图大小是一致的,TensorFlow 会据此计算出高/宽两维度上填充的行/列数。如果是 VALID
则图片经过卷积后会变小,不需要进行任何填充。
以下图为例,输入单通道特征图尺寸为 4 × 4 4\times 4 4×4 ,卷积核尺寸为 3 × 3 3\times 3 3×3 :
-
SAME
填充,在周围填充一圈 0,得到填充后的特征图尺寸为 6 × 6 6\times 6 6×6,再经过卷积,输出特征图尺寸为 4 × 4 4\times 4 4×4 ,与输入尺寸一致; -
VALID
填充,不作任何填充,经过卷积,输出特征图尺寸为 2 × 2 2\times 2 2×2 ,相较于输入尺寸变小。
注:图中数据手动计算得到,下面卷积实现完成后,可比对验证。
填充预处理实现
只有当填充模式为 SAME
的时候需要进行填充,此时填充会保证输入输出特征图的空间尺寸一致,VALID
模式不需要进行填充。
以特征图高为例,根据 SAME
模式中
H
o
u
t
=
H
i
n
H_{out}=H_{in}
Hout=Hin ,计算填充的值
P
P
P :
H
o
u
t
=
(
H
i
n
+
2
P
−
K
)
+
1
P
=
(
K
−
1
)
/
2
H_{out}=(H_{in}+2P-K)\ +\ 1\\ P=(K-1)\ /\ 2
Hout=(Hin+2P−K) + 1P=(K−1) / 2
代码实现:
get_padding
函数获取高/宽两维度的填充行/列数,然后通过 np.pad
完成填充。
import numpy as np
def padding_preprocess(img, kernel_size, mode):
_, in_h, in_w, _ = img.shape
k_h, k_w = kernel_size
padding = get_padding_2d((in_h, in_w), (k_h, k_w), mode)
return np.pad( img, pad_width=padding, mode='constant' )
def get_padding_2d(in_shape, k_shape, mode):
def get_padding_1d(w, k):
if mode == 'SAME':
pads = (w - 1) + k - w
half = pads // 2
padding = (half, half) if pads % 2 == 0 else (half, half + 1)
else:
padding = (0, 0)
return padding
h_pad = get_padding_1d(in_shape[0], k_shape[0])
w_pad = get_padding_1d(in_shape[1], k_shape[1])
return (0, 0), h_pad, w_pad, (0, 0)
测试
同样的 4x4 的输入特征图,卷积核尺寸为 3 ,填充模式为 SAME
:
if __name__ == '__main__':
data16 = [[ 0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11],
[12, 13, 14, 15]]
x = np.array(data16 ).reshape(1, 4, 4, 1)
res = padding_preprocess(x, kernel_size=(3, 3), mode='SAME')
print(np.squeeze(res))
测试的结果:
[[ 0 0 0 0 0 0]
[ 0 0 1 2 3 0]
[ 0 4 5 6 7 0]
[ 0 8 9 10 11 0]
[ 0 12 13 14 15 0]
[ 0 0 0 0 0 0]]
可以看到,与预期是一致的,特征图周围填充了一圈 0。填充后尺寸为 6x6,经过卷积核为 3 的卷积计算,输出特征图的尺寸为 4x4,与输入一致,符合 SAME
模式要求。读者有兴趣的话可以尝试一下更大的卷积核,然后计算一下上述 SAME
模式填充的实现,能否保证卷积输出的特征图是否与输入尺寸一致。
Conv2d前向计算
实现
至此,准备工作都做的差不多了,接下来就可以动手实现卷积的前向计算过程。
具体的实现过程实际上已经在前面介绍组件时介绍过:就是先将输入特征图张量转换为 col 数组,将卷积核也平铺开来,然后二者直接进行矩阵乘法计算。
import numpy as np
def forward(inputs, kernel, stride, padding_mode):
k_h, k_w, _, out_c = kernel['w'].shape
s_h, s_w = stride
# 填充预处理
X = padding_preprocess(inputs, (k_h, k_w), padding_mode)
# im2col
col = im2col(X, k_h, k_w, s_h, s_w)
W = kernel['w'].reshape(-1, out_c)
# 矩阵乘法
Z = col @ W
# 计算并调整尺寸
batch_sz, in_h, in_w, _ = X.shape
out_h = (in_h - k_h) // s_h + 1
out_w = (in_w - k_w) // s_w + 1
Z = Z.reshape(batch_sz, Z.shape[0] // batch_sz, out_c)
Z = Z.reshape(batch_sz, out_h, out_w, out_c)
# 计算偏置
Z += kernel['b']
return Z
测试
还是一直用的例子,设置不再赘述,卷积核我们就简单地设置权重全 1,偏置全 0,方便与上方图中手动计算的结果比较
if __name__ == '__main__':
data16 = [[ 0, 1, 2, 3],
[4, 5, 6, 7],
[8, 9, 10, 11],
[12, 13, 14, 15]]
kernel_size = 3
in_c = 1
out_c = 1
batch_size = 1
stride = 1
weight = np.ones( (kernel_size, kernel_size, in_c, out_c) )
bias = np.zeros( (out_c) )
kernel = dict()
kernel['w'] = weight
kernel['b'] = bias
x = np.array(data16 ).reshape(batch_size, 4, 4, in_c)
k_h, k_w = (3, 3)
s_h, s_w = (1, 1)
y = forward(x, kernel, (stride, stride), 'SAME')
print(np.squeeze(y))
测试结果:
[[10. 18. 24. 18.]
[27. 45. 54. 39.]
[51. 81. 90. 63.]
[42. 66. 72. 50.]]
可以看到,得到的结果与上方图中手动计算的结果一致。
Conv2D反向传播
实现
在反向传播阶段,我们要完成两件事情:
- 计算后一层梯度对于本层权重的梯度,用来更新本层参数
- 计算后一层梯度对于本层输入的梯度,用来将梯度继续传播到前一层
注意:最终实现出的 Conv2D 肯定是一个类(class),其中卷积的权重/偏置参数、步长/填充设置都是以成员变量的形式维护。在前面,为了简明地介绍卷积前向计算的原理和代码实现,所以都是用的函数来讲解和测试。但是在反向传播阶段,由于必须用到一些前向时的上下文信息(如输入特征尺寸等),因此需要这里的 ctx 成员变量时避不开的。但是这不妨碍读者理解反向传播的计算方式,卷积的反向传播函数实现如下。
def backward(self, grad):
# 获取尺寸信息, 由于卷积处理的输入尺寸可变, 因此需要一个在forward和backward时维护一个ctx变量
k_h, k_w, in_c, out_c = self.kernel_shape
s_h, s_w = self.stride
batch_sz, in_h, in_w, in_c = self.ctx['X_shape']
pad_h, pad_w = self.padding[1: 3]
# 计算关于本层权重/偏置参数的梯度
flat_grad = grad.reshape(-1, out_c)
d_W = self.ctx['col'].T @ flat_grad
self.grads['w'] = d_W.reshape(self.kernel_shape)
self.grads['b'] = np.sum(flat_grad, axis=0)
# 计算关于输入的梯度, 使得梯度能够继续反传
d_X = grad @ self.ctx['W'].T
d_in = np.zeros(shape=self.ctx['X_shape'], dtype=np.float32)
for i, r in enumerate(range(0, in_h - k_h + 1, s_h)):
for j, c in enumerate(range(0, in_w - k_w + 1, s_w)):
patch = d_X[:, i, j, :]
patch = patch.reshape(batch_sz, k_h, k_w, in_c)
d_in[:, r: r+k_h, c: c+k_w, :] += patch
# 将关于填充(padding)的梯度裁掉
d_in = d_in[:, pad_h[0]: in_h - pad_h[1], pad_w[0]: in_w - pad_w[1], :]
return d_in
有读者可能会好奇,我们在 Pytorch 中实现自定义层的时候,只需要实现前向 forward
函数就好了呀,为什么这里动手实现 CNN 还需要实现反向传播函数呢?实际上,对于常见算子(如加减乘除等),Pytorch 实现了自动微分机制,我们只要实现自定义层的前向过程,Pytorch 会自动实现它的反向传播过程。关于 Pytroch 的自动微分机制可参考:Pytorch autograd.grad与autograd.backward详解。
有探索精神的读者又会问:那为什么不实现顺便把自动微分实现了呢?是可以的!动手实现自动微分,可参考这篇博客:动手实现一个带自动微分的深度学习框架。文末也给出了该博客对应的自动微分实现方式的最简单的实践仓库,方便读者更好地理解自动微分的机制,并能够参考实现一个简单的自动微分框架。
Conv2D类
在实现了前向计算和反向传播之后,我们的卷积算子就基本完成了。接下来要做的事情就是将上述功能封装成一个 class,并将卷积的权重/偏置参数、梯度、步长/填充设置、上下文信息以成员变量的形式维护。
最终代码如下:
class Conv2D(Layer):
def __init__(self,
kernel,
stride=(1, 1),
padding="SAME",
w_init=XavierUniformInit(),
b_init=ZerosInit()):
super().__init__('Conv2D')
self.kernel_shape = kernel
self.stride = stride
self.initializers = {'w': w_init, 'b': b_init}
self.shapes = {'w': self.kernel_shape, 'b': self.kernel_shape[-1]}
self.padding_mode = padding
self.padding = None
self.is_init = False
def _init_parameters(self):
self.params['w'] = self.initializers['w'](shape=self.shapes['w'])
self.params['b'] = self.initializers['b'](shape=self.shapes['b'])
self.is_init = True
def _inputs_preprocess(self, inputs):
_, in_h, in_w, _ = inputs.shape
k_h, k_w, _, _ = self.kernel_shape
# padding calculation
if self.padding is None:
self.padding = get_padding_2d(
(in_h, in_w), (k_h, k_w), self.padding_mode)
return np.pad(inputs, pad_width=self.padding, mode='constant')
def forward(self, inputs):
if not self.is_init:
self._init_parameters()
k_h, k_w, _, out_c = self.kernel_shape
s_h, s_w = self.stride
X = self._inputs_preprocess(inputs)
col = im2col(X, k_h, k_w, s_h, s_w)
W = self.params['w'].reshape(-1, out_c)
Z = col @ W
batch_sz, in_h, in_w, _ = X.shape
Z = Z.reshape(batch_sz, Z.shape[0] // batch_sz, out_c)
out_h = (in_h - k_h) // s_h + 1
out_w = (in_w - k_w) // s_w + 1
Z = Z.reshape(batch_sz, out_h, out_w, out_c)
Z += self.params['b']
self.ctx = {'X_shape': X.shape, 'col': col, 'W': W}
return Z
def backward(self, grad):
# read size
k_h, k_w, in_c, out_c = self.kernel_shape
s_h, s_w = self.stride
batch_sz, in_h, in_w, in_c = self.ctx['X_shape']
pad_h, pad_w = self.padding[1: 3]
# grads wrt parameters
flat_grad = grad.reshape(-1, out_c)
d_W = self.ctx['col'].T @ flat_grad
self.grads['w'] = d_W.reshape(self.kernel_shape)
self.grads['b'] = np.sum(flat_grad, axis=0)
# grads wrt inputs
d_X = grad @ self.ctx['W'].T
d_in = np.zeros(shape=self.ctx['X_shape'], dtype=np.float32)
for i, r in enumerate(range(0, in_h - k_h + 1, s_h)):
for j, c in enumerate(range(0, in_w - k_w + 1, s_w)):
patch = d_X[:, i, j, :]
patch = patch.reshape(batch_sz, k_h, k_w, in_c)
d_in[:, r: r+k_h, c: c+k_w, :] += patch
# cut off gradients of padding
d_in = d_in[:, pad_h[0]: in_h - pad_h[1], pad_w[0]: in_w - pad_w[1], :]
return d_in
测试代码地址
测试代码
github 仓库地址:https://github.com/Adenialzz/MiniNN/blob/simplestnn/simplestnn/conv_layer.py
测试方式:
git clone git@github.com:Adenialzz/MiniNN.git -b simplestnn
cd simplestnn/
python main.py
自动微分
github 仓库地址:https://github.com/Adenialzz/MiniNN/tree/simplestnn-ad
测试方式:
git clone git@github.com:Adenialzz/MiniNN.git -b simplestnn-ad
cd simplestnn-ad/
python main.py
Ref
- 图像卷积及其计算(特征图尺寸、参数量、计算量)
- 从头搭建一个深度学习框架
- 动手实现一个带自动微分的深度学习框架
- https://github.com/borgwang/tinynn
- Pytorch autograd.grad与autograd.backward详解