[人工智能]动手实现CNN

动手实现CNN

本文将使用 Python 从头到尾实现一个卷积层。当然，既然是自己动手实现，肯定不能用到 Pytorch、TensorFlow 这些现成的库。本文中仅会用到 numpy。文末有 github 仓库地址方便读者测试，仓库中会测试我们自己实现的 CNN 在 MNIST 手写数字识别上的准确率。

卷积的基本计算

先来回顾一下卷积基本计算，分别计算输出特征图的尺寸、卷积的计算量、卷积的参数量。这一小节就直接给出公式并简单解释一下，详细的计算细节和代码统计参数量/计算量的方法可参考：图像卷积及其计算（特征图尺寸、参数量、计算量）

记输入特征图尺寸 $H_{in}\times W_{in}\times C_{in}$ ，卷积核尺寸为 $K\times K$ ，卷积层填充为 $P$，步长为 $S$， 输出特征图通道数为：$C_{out}$ 。

输出特征图尺寸

按照卷积的滑窗方法计算即可。
$$\begin{gathered} H_{out}=(H_{in}+2P?K)//2+1 \\ {W}_{out}=(W_{in}+2P?K)//2+1 \end{gathered}$$
参数量

每个卷积核的参数量是 $C_{in}\times K\times K$ ，共 $C_{out}$ 个卷积核。
$$C_{in}\times C_{out}\times K\times K$$
计算量

先算出输出特征图每个位置的计算量，然后再乘以输出特征图的尺寸即可。前面括号中就是每个位置的计算量，其中前一项是乘法计算量，后一项是加法计算量，注意这里没有考虑偏置，如果有偏置需要再加一。后面括号是特征图的尺寸。
$$(C_{in}\times K\times K+C_{in}\times K\times K?1)\times (C_{out}\times H_{out}\times W_{out})$$

im2col

怎样实现卷积？

在初学卷积原理时，我们知道，CV 中的卷积就是拿一个卷积核，以滑动窗口的形式去特征图上进行计算，每个位置的计算结果就是输出特征图的一个元素（如下面动图所示）。这样理解当然是没问题的，但是在实际实现时，并不会按照这种方式去滑窗、取数、计算。为什么呢？访存开销太大。我们知道，现代计算机中 CPU 的频率比内存快得多，因此，在程序中，访存的开销要比实际计算的开销大得多。保证规律的、高效的访存，能大幅提升程序的性能。

如果按照滑窗卷积进行计算，在计算机操作时，需要将其存入内存当中再操作（按照“行先序”）：

这样一来，特征图y11,y12,y21,y22的计算如下所示：

按照卷积计算，实线标注出卷积计算中的访存过程（对应数据相乘），我们可以看到这一过程是非常散乱和混乱的。直接用卷积的计算方式是比较愚蠢的。

这时候就要用到im2col操作。

im2col原理

im2col 是卷积实现中的关键一步。一句话来介绍im2col操作的话，就是通过牺牲空间的手段（约扩增K×K倍），将特征图转换成庞大的矩阵来进行卷积计算。

其实思路非常简单：

把每一次循环所需要的数据都排列成列向量，然后逐一堆叠起来形成矩阵（按通道顺序在列方向上拼接矩阵）。

比如 $C_{in}\times H_{in}\times W_{in}$大小的输入特征图，$K\times K$ 大小的卷积核，输出大小为 $C_{out}\times W_{out}\times H_{out}$，输入特征图将按需求被转换成 $(K\times K\times C_{in})\times (W_{out}\times H_{out})$的矩阵，卷积核将被转换成 $C_{out}\times (K\times K\times C_{in})$ 的矩阵，这样二者直接进行矩阵乘法计算就实现了卷积。

如下图中：

• 对于输入特征图：将每个将要与卷积核进行计算的滑动窗口抽出来，作为一列向量，组成 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\times 4$ ，卷积核尺寸为 $3\times 3$ ：

• SAME 填充，在周围填充一圈 0，得到填充后的特征图尺寸为 $6\times 6$，再经过卷积，输出特征图尺寸为 $4\times 4$ ，与输入尺寸一致；
• VALID 填充，不作任何填充，经过卷积，输出特征图尺寸为 $2\times 2$ ，相较于输入尺寸变小。

注：图中数据手动计算得到，下面卷积实现完成后，可比对验证。

填充预处理实现

只有当填充模式为 SAME 的时候需要进行填充，此时填充会保证输入输出特征图的空间尺寸一致，VALID 模式不需要进行填充。

以特征图高为例，根据 SAME 模式中 $H_{out}=H_{in}$ ，计算填充的值 $P$ ：
$$\begin{gathered} H_{out}=(H_{in}+2P?K)+1 \\ P=(K?1)/2 \end{gathered}$$
代码实现：

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

1. 图像卷积及其计算（特征图尺寸、参数量、计算量）
2. 从头搭建一个深度学习框架
3. 动手实现一个带自动微分的深度学习框架
4. https://github.com/borgwang/tinynn
5. Pytorch autograd.grad与autograd.backward详解