动手实现CNN

时间:2022-09-28 12:55:26

动手实现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+2PK) // 2 + 1Wout=(Win+2PK) // 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×K1)×(Cout×Hout×Wout)

im2col

怎样实现卷积?

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

动手实现CNN

如果按照滑窗卷积进行计算,在计算机操作时,需要将其存入内存当中再操作(按照“行先序”):

动手实现CNN

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

动手实现CNN

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

这时候就要用到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 等)中都已经经过了极致的优化,有很高的性能。

动手实现CNN

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 中的填充模式右 SAMEVALID 两种。如果填充模式设置为 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 ,相较于输入尺寸变小。

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

动手实现CNN

填充预处理实现

只有当填充模式为 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+2PK) + 1P=(K1) / 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

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