卷积神经网络 (convolution neural network) 是当下图像识别最为热门而有效的手段。这种网络一般有若干个卷积模型 (convolution model) 加上若干个线性模型 (linear model) 组合而成。卷积模型部分相当于提取图像的特征,而线性模型部分是根据这些特征做出预测。卷积模型的想法是通过卷积运算提取图像的局部信息。我们在本笔记中首先介绍卷积运行,然后按照 CS231n Convolutional Neural Networks for Visual Recognition 这篇参考资料的思路实现卷积模型。我们会在下一篇笔记中使用卷积模型建立一个简单的卷积神经网络来做 Fashion MNIST 识别。

import numpy as np
import matplotlib.pyplot as plt

%run import_npnet.py
npnet = import_npnet(5)

卷积运算

给定两个矩阵 \(A = (A_{ij})\)\(B = (B_{kl})\),两者的卷积 (convolution) \(C = B \ast A\) 定义为

$$C_{ij} = \sum_{k,l} B_{kl}A_{i+k-1, j+l-1}.$$


其中 \(B\) 称为卷积核 (kernel)。

比如我们有两个矩阵

$$A = \begin{pmatrix} 1 & 3 & 5 \\ 3 & 2 & 1 \\ 0 & 1 & 0 \end{pmatrix}, \hspace{5mm} B = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}.$$


\(B \ast A\) 的结果是一个 \(2 \times 2\) 的矩阵

$$C = B \ast A = \begin{pmatrix} 0 & 3 \\ 2 & 0 \end{pmatrix}.$$


\(C\) 中的第一个 0 是 \(A\) 中第一个 \(2 \times 2\) 子矩阵 \(

\begin{pmatrix} 1 & 3\\ 3 & 2\end{pmatrix}
\) 与 \(B\) 按位乘法得到数字之和。更细致地展开来说

$$ 1 \times 0 + 3 \times 1 + 3 \times (-1) + 2 \times 0 = 0.$$

接着我们向右移动 \(A\)\(2 \times 2\) 子矩阵的位置,得到 \(A\) 的下一个 \(2 \times 2\) 子矩阵 \(

\begin{pmatrix} 3 & 5 \\ 2 & 1\end{pmatrix}
\),其跟 \(B\) 的按位乘法得到数字之和为3。当我们无法向右移子矩阵的位置时,我们返回最左边然后向下一移一行,得到 \(A\) 的子矩阵 \(
\begin{pmatrix} 3 & 2 \\ 0 & 1 \end{pmatrix}
\)
。这个子矩阵跟 \(B\) 的卷积结果时2。最后一个子矩阵时 \(
\begin{pmatrix} 2 & 1 \\ 1 & 0 \end{pmatrix}
\)
,跟 \(B\) 的卷积结果是0。我们不难发现,上述的卷积运算可以化作矩阵乘法。首先我们将所有子矩阵变成行向量得到矩阵 \(A^\prime\),然后将 \(B\) 变成列向量 \(B^\prime\)

$$A' = \begin{pmatrix} 1 & 3 & 3 & 2 \\ 3 & 5 & 2 & 1 \\ 3 & 2 & 0 & 1 \\ 2 & 1 & 1 & 0 \end{pmatrix}, \hspace{5mm} B' = \begin{pmatrix} 0 \\ 1 \\ -1 \\ 0 \end{pmatrix}.$$

那么将矩阵乘法 \(A^\prime \cdot B^\prime\) 得到的结果化成 \(2 \times 2\) 矩阵就能得到 \(C = B \ast A\)。我们按照这个思路,首先搞清楚如何将矩阵的所有子矩阵合在一起变成一个大的矩阵,以及如何将矩阵乘法的结果变回卷积的结果。在这之前我们介绍几个卷积运算中重要的概念。在卷积运算中每次向左或者向下移动的数量叫做步长 (stride),我们记向左的步长为 \(sW\),向下的步长为 \(sH\)。在卷积运算之前我们常对矩阵 \(A\) 进行向四周填充0 (zero padding),我们将向上下填充的行数叫做 \(pH\),向左右填充的列数叫做 \(pW\)。比如我们对上述 \(A\) 进行 \(pH = pW = 1\) 的0填充能得到

$$\begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 3 & 5 & 0 \\ 0 & 3 & 2 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix}.$$

im2row

在卷积模型中,输入是一个四维 \(N \times C \times H \times W\) 的矩阵。其中 \(N\) 表示图片的个数,\(C\) 表示图片的频道 (channel) 数目,\(H\) 代表图片的高度,而 \(W\) 代表图片的宽度。比如1张普通的 RGB \(1920 \times 1080\) 的照片是一个 \(1 \times 3 \times 1920 \times 1080\) 的矩阵。假设卷积模型中有 \(D\) 个卷积核,每个卷积核的宽为 \(kW\) 且高为 \(kH\)。那么卷积得到的结果是一个 \(N \times D \times H^\prime \times W^\prime\) 的思维矩阵,其中

$$\begin{align*} & H^\prime = \frac{H + 2 * pH - kH}{sH} + 1, \\ & W^\prime = \frac{W + 2 * pW - kW}{sW} + 1. \end{align*}$$

卷积模型中的卷积运算是对前一节卷积运算的推广。上一节中的卷积运算相当于 \(N = C = 1\) 的特殊情形。现在我们的图片是有 \(C\) 层频道的,所以每个卷积核是一个 \(C \times kH \times kW\) 的三维矩阵。如果我们希望将卷积运算化作矩阵乘法,我们需要将输入和卷积核变成二维矩阵。将 \(D\) 个这样的卷积核变成一个二维矩阵是不难的,这个二维矩阵是一个规模为 \((C*kH*kW) \times D\) 的矩阵,每一列代表一个卷积核。困难的是如何将 \(N \times C \times H \times W\) 的四维矩阵变成一个二维矩阵。对输入规模和输出规模来分析,我们需要将这个四维矩阵变成一个规模为 \((N*H^\prime*W^\prime) \times (C*kH*kW)\) 的二维矩阵。假设我们有一个简单的 \(2 \times 2 \times 2 \times 3\) 输入:

img = [
    [
        [[1, 2, 3],
         [4, 5, 6]],

        [[6, 5, 4],
         [3, 2, 1]]
    ],
    [
        [[0, 1, 0],
         [2, 2, 3]],

        [[2, 3, 3],
         [6, 6, 6]]
    ]
]

我们假设 \(kH = kW = 2, sH = sW = 1\)\(pH = pW = 0\)。我们不难计算出 \(H^\prime = 1, W^\prime = 2\)。我们希望得到的对应的二维数组为

row = [
    [1, 2, 4, 5, 6, 5, 3, 2],
    [2, 3, 5, 6, 5, 4, 2, 1],
    [0, 1, 2, 2, 2, 3, 6, 6],
    [1, 0, 2, 3, 3, 3, 6, 6]
]

row 的第一行是 img[0] 所有层中第一个 \(2 \times 2\) 矩阵得到的数字所构成的行向量,第二行是 img[0] 所有层中第二个 \(2 \times 2\) 矩阵得到的数字所构成的行向量。第三第四行是从 img[1] 得来的。如果我们分析 row 中数字来自于 img 的索引 (index),我们可以得到

row_indices = [
    0[000, 001, 010, 011, 100, 101, 110, 111],
    0[001, 002, 011, 012, 101, 102, 111, 112],
    1[000, 001, 010, 011, 100, 101, 110, 111],
    1[001, 002, 011, 012, 101, 102, 111, 112]
]

n[...chw...] 此处的数字来自于 img[n, c, h, w]。我们发现这些索引是很有规律的,所以我们可以通过 numpy 的高级所以 (advanced indexing) 来将 img 变成 row

# import
def im2row_indices(img_shape, kernel, stride, padding):
    N, C, H, W = img_shape
    kH, kW = kernel
    sH, sW = stride
    pH, pW = padding
    assert (H + 2 * pH - kH) % sH == 0
    assert (W + 2 * pW - kW) % sW == 0
    oH = (H + 2 * pH - kH) // sH + 1
    oW = (W + 2 * pW - kW) // sW + 1

    # k 代表的是 chw 中 c 的变化规律
    k = np.repeat(np.arange(C), kH * kW)

    # j 代表的是 chw 中 h 的变化规律
    j0 = np.repeat(np.arange(kH), kW)
    j0 = np.tile(j0, C)
    j1 = sH * np.repeat(np.arange(oH), oW)
    j = j0 + j1.reshape(-1, 1)

    # i 代表的是 chw 中 w 的变化规律
    i0 = np.tile(np.arange(kW), kH * C)
    i1 = sW * np.tile(np.arange(oW), oH)
    i = i0 + i1.reshape(-1, 1)

    return i, j, k


def im2row(img, kernel, stride, padding):
    C = img.shape[1]
    kH, kW = kernel
    pH, pW = padding

    i, j, k = im2row_indices(img.shape, kernel, stride, padding)
    img = np.pad(img, ((0, 0), (0, 0), (pH, pH), (pW, pW)), mode='constant')
    rows = img[:, k, j, i]
    rows = rows.reshape(-1, kW * kH * C)
    return rows

我们验证一下之前给出的例子:

img = [
    [
        [[1, 2, 3],
         [4, 5, 6]],

        [[6, 5, 4],
         [3, 2, 1]]
    ],
    [
        [[0, 1, 0],
         [2, 2, 3]],

        [[2, 3, 3],
         [6, 6, 6]]
    ]
]
img = np.array(img)
row = im2row(img, (2, 2), (1, 1), (0, 0))
print(row)
[[1 2 4 5 6 5 3 2]
 [2 3 5 6 5 4 2 1]
 [0 1 2 2 2 3 6 6]
 [1 0 2 3 3 3 6 6]]

Nice! 现在一个四维的输入通过上述的 im2row 变成一个二维的数组。那么当我们计算向后传播的导数时得到的是一个关于这个二维数组的导数。我们需要将这个二维数组的导数还原成关于四维输入的导数。使用链式法则 (chain rule) 以及 im2row_indices,我们能得到下列这个将二维数组导数还原成四维输入导数的函数。

# import
def row2im(row, img_shape, kernel, stride, padding):
    N, C, H, W = img_shape
    kH, kW = kernel
    pH, pW = padding
    H += 2 * pH
    W += 2 * pW
    img = np.zeros((N, C, H, W), dtype=row.dtype)

    i, j, k = im2row_indices(img_shape, kernel, stride, padding)
    row = row.reshape(N, -1, kW * kH * C)
    np.add.at(img, (slice(None), k, j, i), row)

    return img[:, :, pH:H-pH, pW:W-pW]

需要注意的是,row2im 并不是将 row 转成原来的 img,所以它不是 im2row 的逆函数。

卷积模型

通过 im2rowrow2im,卷积模型其实就相当于一个线性模型来。使用之前的符号,卷积模型需要的参数有:输入图片的层数 \(C\) ,输出图片的层数 \(D\),每一层卷积核的高和宽 \(kH, kW\),步长的高和宽 \(sH, sW\),和0填充的高和宽 \(pH, pW\)。模型的参数是一个 \((C*kH*kW) \times D\) 的二维矩阵,每一列代表一个卷积核,还有大小为 \(D\) 的偏移量 (bias)。

# import
def to_tuple(x):
    if isinstance(x, tuple):
        ans = x
    elif isinstance(x, list):
        ans = tuple(x)
    else:
        ans = (x, x)
    assert len(ans) == 2
    return ans


class Conv(npnet.Model):

    def __init__(self, in_channels, out_channels,
                 kernel=3, stride=1, padding=1):
        super().__init__()
        self.in_channels = in_channels
        self.out_channels = out_channels
        self.kernel = to_tuple(kernel)
        self.stride = to_tuple(stride)
        self.padding = to_tuple(padding)

        ksize = in_channels * self.kernel[0] * self.kernel[1]
        self.parameters['W'] = np.random.randn(ksize, out_channels)
        self.parameters['W'] *= np.sqrt(2.0 / (in_channels + out_channels))
        self.parameters['b'] = np.random.randn(out_channels)

        self._init_gradients()

    def forward(self, input):
        self.input_shape = input.shape
        params = self.parameters
        row = im2row(input, self.kernel, self.stride, self.padding)
        self.row = row
        output = row @ params['W'] + params['b']
        output = output.transpose()

        # 正确地将二维数组转成四维数组作为卷积的结果
        N, D, H, W = self._get_output_shape(input.shape)
        output = output.reshape(D, N, H, W)
        output = output.transpose(1, 0, 2, 3)
        return output

    def backward(self, grad):
        grad = grad.transpose(1, 0, 2, 3)
        grad = grad.reshape(grad.shape[0], -1)
        grad = grad.transpose()

        self.gradients['W'] += self.row.T @ grad
        self.gradients['b'] += np.sum(grad, axis=0)
        grad = grad @ self.parameters['W'].T
        grad = row2im(grad,
                      self.input_shape,
                      self.kernel,
                      self.stride,
                      self.padding)
        return grad

    def _get_output_shape(self, input_shape):
        N, _, H, W = input_shape
        kH, kW = self.kernel
        sH, sW = self.stride
        pH, pW = self.padding
        assert (H + 2 * pH - kH) % sH == 0
        assert (W + 2 * pW - kW) % sW == 0
        oH = (H + 2 * pH - kH) // sH + 1
        oW = (W + 2 * pW - kW) // sW + 1
        return N, self.out_channels, oH, oW

跟往常一样,我们需要验证偏导数是否正确。

conv = Conv(4, 8)
npnet.GradientCheck(conv, input_shape=(20, 4, 28, 28), allowed_error=1e-5).check()
True

池化

池化 (pooling) 的操作跟卷积很像,也都是考虑图片的局部信息。不同于卷积层 (convolution layer) 的是,池化层 (pooling layer) 没有参数。常见的池化有最大池化 (max pooling) 以及平均池化 (average pooling)。最大池化是取局部的最大值,而平均池化是取局部的平均值。平均池化可以看作一种特殊的卷积。

池化有很多作用,比如可以降低图片的大小,增加不变性等等。用兴趣的可以参考知乎的一个问答:CNN网络的pooling层有什么用?。使用 im2rowrow2im,我们能很容易的封装一个抽象的池化类:

# import
class Pooling(npnet.Model):

    def __init__(self, pool_func, pool_der,
                 kernel=2, stride=2, padding=0):
        super().__init__()
        self.pool_func = pool_func
        self.pool_der = pool_der
        self.kernel = to_tuple(kernel)
        self.stride = to_tuple(stride)
        self.padding = to_tuple(padding)

    def forward(self, input):
        self.input_shape = input.shape
        N, C, H, W = input.shape
        input = input.reshape(N * C, 1, H, W)
        row = im2row(input, self.kernel, self.stride, self.padding)
        self.row = row
        output = self.pool_func(row)
        output = output.transpose()
        N, C, H, W = self._get_output_shape(self.input_shape)
        output = output.reshape(N, C, H, W)
        return output

    def backward(self, grad):
        grad = grad.reshape(1, -1)
        grad = grad.transpose()
        grad = self.pool_der(self.row, grad)
        N, C, H, W = self.input_shape
        shape = (N * C, 1, H, W)
        grad = row2im(grad,
                      shape,
                      self.kernel,
                      self.stride,
                      self.padding)
        grad = grad.reshape(N, C, H, W)
        return grad

    def _get_output_shape(self, input_shape):
        N, C, H, W = input_shape
        kH, kW = self.kernel
        sH, sW = self.stride
        pH, pW = self.padding
        assert (H + 2 * pH - kH) % sH == 0
        assert (W + 2 * pW - kW) % sW == 0
        oH = (H + 2 * pH - kH) // sH + 1
        oW = (W + 2 * pW - kW) // sW + 1
        return N, C, oH, oW

有了这个池化的基类 Pooling,我们只需要定制 pool_funcpool_der 就能定义特定的池化层。

最大池化

# import
class MaxPool(Pooling):

    def __init__(self, kernel=2, stride=2, padding=0):
        def max_func(row):
            out = np.max(row, axis=1)
            return out

        def max_der(row, grad):
            cache = np.argmax(row, axis=1)
            out = np.zeros_like(row)
            out[np.arange(cache.size), cache] = grad.ravel()
            return out

        super().__init__(pool_func=max_func,
                         pool_der=max_der,
                         kernel=kernel,
                         stride=stride,
                         padding=padding)
npnet.GradientCheck(MaxPool(), input_shape=(10, 3, 10, 10)).check()
True

平均池化

# import
class AvgPool(Pooling):

    def __init__(self, kernel=2, stride=2, padding=0):
        def avg_func(row):
            out = np.average(row, axis=1)
            return out

        def avg_der(row, grad):
            m = row.shape[1]
            grad = grad / m
            grad = np.repeat(grad, m, axis=1)
            return grad

        super().__init__(pool_func=avg_func,
                         pool_der=avg_der,
                         kernel=kernel,
                         stride=stride,
                         padding=padding)
npnet.GradientCheck(AvgPool(), input_shape=(10, 3, 10, 10)).check()
True

我们已经将需要的模型准备好了,下一篇笔记将使用这些模型搭建卷积神经网络来完成 Fashion MNIST 的识别任务。