【455】Python 徒手實現 卷積神經網絡 CNN

參考：CNNs, Part 1: An Introduction to Convolutional Neural Networks

參考：CNNs, Part 2: Training a Convolutional Neural Network

---

目錄

1. 動機（Motivation）
2. 數據集（Dataset）
3. 卷積（Convolutions）
4. 池化（Pooling）
5. Softmax
6. 訓練概述（Training Overview）
7. 反向傳播：Softmax（Backprop: Softmax）
8. 反向傳播：池化層（Backprop: Max Pooling）
9. 反向傳播：卷積層（Backprop: Conv）
10. 訓練 CNN（Training a CNN）
11. Keras 實現

---

1. 動機（Motivation）

　　通過普通的神經網絡可以實現，但是現在圖片越來越大，如果通過 NN 來實現，訓練的參數太多。例如 224 x 224 x 3 = 150,528，隱藏層設置為 1024 就需要訓練參數 150,528 x 1024 = 1.5 億 個，這還是第一層，因此會導致我們的網絡很龐大。

　　另一個問題就是特征位置在不同的圖片中會發生變化。例如小貓的臉在不同圖片中可能位於左上角或者右下角，因此小貓的臉不會激活同一個神經元。

2. 數據集（Dataset）

　　我們使用手寫數字數據集 MNIST 。

　　

　　每個數據集都以一個 28x28 像素的數字。

　　普通的神經網絡也可以處理這個數據集，因為圖片較小，另外數字都集中在中間位置，但是現實世界中的圖片分類問題可就沒有這么簡單了，這里只是拋磚引玉哈。

3. 卷積（Convolutions）

　　CNN 相較於 NN 來說主要是增加了基於 convolution 的卷積層。卷基層包含一組 filter，每一個 filter 都是一個 2 維的矩陣。以下為 3x3 filter：

　　

　　我們可以通過輸入的圖片和上面的 filter 來做卷積運算，然后輸出一個新的圖片。包含以下步驟：

• • 將 filter 疊加在圖片的頂部，一般是左上角
  • 然后執行對應元素的相乘
  • 將相乘的結果進行求和，得到輸出圖片的目標像素值
  • 重復以上操作在所有位置上

　　執行效果如下所示：

　　

　　3.1 有用嗎？

　　通過卷積可以提取圖片中的特定線條，垂直線條或者水平線條，以下為 vertical Sobel filter and horizontal Sobel filter 的結果：

　　  

　　卷積可以幫助我們查找一些圖片特征（例如邊緣）。

　　3.2 Padding（填充）

　　可以通過在周圍補 0 實現輸出前后圖像大小一致，如下所示：

　　

　　這叫做 "same padding"，不過一般不用 padding，叫做 "valid" padding。

　　3.3 卷基層

　　CNN 包含卷基層，卷基層通過一組 filter 將輸入的圖片轉為輸出的圖片。卷基層的主要參數是 filter 的個數。

　　對於 MNIST CNN，我使用一個含有 8 個 filter 的卷基層，意味着它將 28x28 的輸入圖片轉為 26x26x8 的輸出集：

　　

　　卷基層的 8 個 filter 分別產生 26x26 的輸出，只有 3 x 3 (filter size) x 8 (nb_filters) = 72 權重值。

　　3.4 卷積層代碼實現

　　簡單起見，我們使用 3x3 的filter，首先實現一個 卷基層的類：

import numpy as np

class Conv3x3:
    # A Convolution layer using 3x3 filters.

    def __init__(self, num_filters):
        self.num_filters = num_filters

        # filters is a 3d array with dimensions (num_filters, 3, 3)
        # We divide by 9 to reduce the variance of our initial values
        self.filters = np.random.randn(num_filters, 3, 3) / 9

　　Conv3x3 類只需要一個參數：filter 個數。通過 NumPy 的 randn() 方法實現。之所以在初始化的時候除以 9 是因為對於初始化的值不能太大也不能太小，參考：Xavier Initialization。 

　　接下來，具體實現卷基層：

class Conv3x3:
    # ...

    def iterate_regions(self, image):
        '''
        Generates all possible 3x3 image regions using valid padding.
        - image is a 2d numpy array
        '''
        h, w = image.shape

        for i in range(h - 2):
            for j in range(w - 2):
                im_region = image[i:(i + 3), j:(j + 3)]
                yield im_region, i, j
        # 將 im_region, i, j 以 tuple 形式存儲到迭代器中
        # 以便后面遍歷使用

    def forward(self, input):
        '''
        Performs a forward pass of the conv layer using the given input.
        Returns a 3d numpy array with dimensions (h, w, num_filters).
        - input is a 2d numpy array
        '''
        # input 為 image，即輸入數據
        # output 為輸出框架，默認都為 0，都為 1 也可以，反正后面會覆蓋
        # input: 28x28
        # output: 26x26x8
        h, w = input.shape
        output = np.zeros((h - 2, w - 2, self.num_filters))

        for im_region, i, j in self.iterate_regions(input):
            # 卷積運算，點乘再相加，ouput[i, j] 為向量，8 層
            output[i, j] = np.sum(im_region * self.filters, axis=(1, 2))
        # 最后將輸出數據返回，便於下一層的輸入使用
        return output

4. 池化（Pooling） 　　

　　圖片的相鄰像素具有相似的值，因此卷基層中很多信息是冗余的。通過池化來減少這個影響，包含 max, min or average，下圖為基於 2x2 的 Max Pooling：

　　

　　與卷積計算類似，只是這個更容易，只是計算最大值並賦值。池化層將會把 26x26x8 的輸入轉為 13x13x8 的輸出：

　　

　　4.1 池化層代碼實現

import numpy as np

class MaxPool2:
    # A Max Pooling layer using a pool size of 2.

    def iterate_regions(self, image):
        '''
        Generates non-overlapping 2x2 image regions to pool over.
        - image is a 2d numpy array
        '''
        # image: 26x26x8
        h, w, _ = image.shape
        new_h = h // 2
        new_w = w // 2

        for i in range(new_h):
            for j in range(new_w):
                im_region = image[(i * 2):(i * 2 + 2), (j * 2):(j * 2 + 2)]
                yield im_region, i, j

    def forward(self, input):
        '''
        Performs a forward pass of the maxpool layer using the given input.
        Returns a 3d numpy array with dimensions (h / 2, w / 2, num_filters).
        - input is a 3d numpy array with dimensions (h, w, num_filters)
        '''
        # input: 卷基層的輸出，池化層的輸入
        h, w, num_filters = input.shape
        output = np.zeros((h // 2, w // 2, num_filters))

        for im_region, i, j in self.iterate_regions(input):
            output[i, j] = np.amax(im_region, axis=(0, 1))
        return output

5. Softmax

　　為了完成我們的 CNN，我們需要進行具體的預測。通過 softmax 來實現，將一組數字轉換為一組概率，總和為 1。參考：Softmax function。

　　5.1 用法

　　我們將要使用一個含有 10 個節點（分別代表相應數字）的 softmax 層，作為我們 CNN 的最后一層。最后一層為一個全連接層，只是激活函數為 softmax。經過 softmax 的變換，數字就是具有最高概率的節點。

　　

　　softmax 為 13x13x8 轉換為一列節點后與 10 個節點組成一個全連接，然后 softmax 為激活函數。

　　5.2 交叉熵損失函數（Cross-Entropy Loss）

　　交叉熵用來計算概率間的距離，具體公式可參考：筆記 | 什么是Cross Entropy。

$$H(p, q)=-\sum_{x} p(x) ln(q(x))$$

　　其中：

• $p(x)$ 為真實概率
• $q(x)$ 為預測概率
• $H(p, q)$ 為預測結果與真實結果的差距

　　在我們的具體問題中，對於真實概率，只有分類正確數字對應的概率為 1，其他均為 0，因此 交叉熵損失函數 可以寫成如下形式：

$$L=-ln(p_c)$$

　　其中，$c$ 是正確分類（本例中為正確的數字），$p_c$ 是 $c$ 類的預測概率。$L$ 的值越小越好。

　　5.3 Softmax 層代碼實現

import numpy as np

class Softmax:
    # A standard fully-connected layer with softmax activation.

    def __init__(self, input_len, nodes):
        # We divide by input_len to reduce the variance of our initial values
        # input_len: 輸入層的節點個數，池化層輸出拉平之后的
        # nodes: 輸出層的節點個數，本例中為 10
        # 構建權重矩陣，初始化隨機數，不能太大
        self.weights = np.random.randn(input_len, nodes) / input_len
        self.biases = np.zeros(nodes)

    def forward(self, input):
        '''
        Performs a forward pass of the softmax layer using the given input.
        Returns a 1d numpy array containing the respective probability values.
        - input can be any array with any dimensions.
        '''
        # 3d to 1d，用來構建全連接網絡
        input = input.flatten()

        input_len, nodes = self.weights.shape

        # input: 13x13x8 = 1352
        # self.weights: (1352, 10)
        # 以上叉乘之后為 向量，1352個節點與對應的權重相乘再加上bias得到輸出的節點
        # totals: 向量, 10
        totals = np.dot(input, self.weights) + self.biases
        # exp: 向量, 10
        exp = np.exp(totals)
        return exp / np.sum(exp, axis=0)

　　至此，我們完成了我們 CNN 模型的整個 forward pass！把它們放在一起調用：

import mnist
import numpy as np

# We only use the first 1k testing examples (out of 10k total)
# in the interest of time. Feel free to change this if you want.
test_images = mnist.test_images()[:1000]
test_labels = mnist.test_labels()[:1000]

conv = Conv3x3(8)                                    # 28x28x1 -> 26x26x8
pool = MaxPool2()                                    # 26x26x8 -> 13x13x8
softmax = Softmax(13 * 13 * 8, 10) # 13x13x8 -> 10

def forward(image, label):
    '''
    Completes a forward pass of the CNN and calculates the accuracy and
    cross-entropy loss.
    - image is a 2d numpy array
    - label is a digit
    '''
    # We transform the image from [0, 255] to [-0.5, 0.5] to make it easier
    # to work with. This is standard practice.
 
   # out 為卷基層的輸出, 26x26x8
    out = conv.forward((image / 255) - 0.5)
    # out 為池化層的輸出, 13x13x8
    out = pool.forward(out)
    # out 為 softmax 的輸出, 10
    out = softmax.forward(out)

    # Calculate cross-entropy loss and accuracy. np.log() is the natural log.
    # 損失函數的計算只與 label 的數有關，相當於索引
    loss = -np.log(out[label])
    # 如果 softmax 輸出的最大值就是 label 的值，表示正確，否則錯誤
    acc = 1 if np.argmax(out) == label else 0

    return out, loss, acc

print('MNIST CNN initialized!')

loss = 0
num_correct = 0
# enumerate 函數用來增加索引值
for i, (im, label) in enumerate(zip(test_images, test_labels)):
    # Do a forward pass.
    _, l, acc = forward(im, label)
    loss += l
    num_correct += acc

    # Print stats every 100 steps.
    if i % 100 == 99:
        print(
            '[Step %d] Past 100 steps: Average Loss %.3f | Accuracy: %d%%' %
            (i + 1, loss / 100, num_correct)
        )
        loss = 0
        num_correct = 0

　　輸出結果如下所示：

MNIST CNN initialized!
[Step 100] Past 100 steps: Average Loss 2.302 | Accuracy: 11%
[Step 200] Past 100 steps: Average Loss 2.302 | Accuracy: 8%
[Step 300] Past 100 steps: Average Loss 2.302 | Accuracy: 3%
[Step 400] Past 100 steps: Average Loss 2.302 | Accuracy: 12%

　　這也比較合理，由於是通過隨機的權重初始值，目前這個 CNN 模型跟我們隨機猜測的結果類似。隨機猜測的結果是 10%。

6. 訓練概述（Training Overview）

　　訓練神經網絡一般包含兩個階段：

• forward phase: 輸入參數傳遞通過整個網絡。
• backward phase: 反向傳播更新 gradient 和 weight。

　　我們按照如上的模式來訓練 CNN。還有以下兩個方法需要使用：

• 在 forward phase 中，每一層都需要存儲一些數據（例如輸入數據，中間值等）。這些數據將會在 backward phase 中得到使用。因此每一個 backward phase 都需要在相應的 forward phase 之后運行。
• 在 backward phase 中，每一層都要獲取 gradient 並且也返回 gradient。獲取的是 loss 對於該層輸出（$\frac{\partial L}{\partial out}$）的 gradient，返回的是 loss 對於該層輸入（$\frac{\partial L}{\partial in}$）的 gradient。

　　上面兩個方法可以幫助我們更有條理且簡潔的實現訓練。訓練 CNN 的代碼大約長下面的樣紙：

# Feed forward

# image 為輸入層，28x28
# out 為卷基層輸出，26x26x8
out = conv.forward((image / 255) - 0.5)
# out 為池化層輸出，13x13x8
out = pool.forward(out)
# out 為 softmax 層輸出，10
out = softmax.forward(out)

# Calculate initial gradient
# gradient: loss 對於 softmax 輸出層的 gradient
gradient = np.zeros(10)
# ...

# Backprop

# gradient：loss 對於 softmax 輸入層的 gradient
# 輸入為 loss 對於 softmax 輸出層的 gradient
gradient = softmax.backprop(gradient)
# gradient：loss 對於池化層輸入層的 gradient
# 輸入為 loss 對於池化層輸出層的 gradient
gradient = pool.backprop(gradient)
# gradient：loss 對於卷基層輸入層的 gradient
# 輸入為 loss 對於卷基層輸出層的 gradient
gradient = conv.backprop(gradient)

7. 反向傳播：Softmax（Backprop: Softmax）

　　我們需要從最后開始朝着最前面計算，這就是 backprop 的工作原理。首先回想下交叉熵損失函數（cross-entropy loss）：

$$L=-ln(p_c)$$

　　其中，$p_c$ 是正確類 $c$ （也就是圖片中的數字）的預測概率。

　　首先我們需要計算 softmax 層的 backward phase 的輸入數據，$\frac{\partial L}{\partial out_s}$，其中 $out_s$ （下標的 $s$ 是說明 softmax 層）是指 softmax 層的輸出值：一個含有 10 個概率值的向量。由於 $p_i$ 只出現在了 loss 方程中，因此很容易計算：
$$
\begin{equation}
\frac{\partial L}{\partial out_s(i)}=
\begin{cases}
0& \text{if } i\neq c\\
-\frac{1}{p_i}& \text{if } i=c
\end{cases}
\end{equation}
$$

　　上面就是我們的初始化 gradient：

# Calculate initial gradient
# 默認都為 0
gradient = np.zeros(10)
# 只修改 label 值對應的
gradient[label] = -1 / out[label]

　　現在我們已經准備好了開始實現我們第一個 backward phase，但是我們需要首先在 forward phase 中存儲我們前面討論的相關數據。

class Softmax:
    # ...

    def forward(self, input):
        '''
        Performs a forward pass of the softmax layer using the given input.
        Returns a 1d numpy array containing the respective probability values.
        - input can be any array with any dimensions.
        '''

        # NEW ADD，13x13x8
        self.last_input_shape = input.shape

        input = input.flatten()

        # NEW ADD, 向量，1352
        self.last_input = input

        input_len, nodes = self.weights.shape

        totals = np.dot(input, self.weights) + self.biases

        # NEW ADD，softmax 前的向量，10
        self.last_totals = totals

        exp = np.exp(totals)
        return exp / np.sum(exp, axis=0)

　　接下來我們可以獲取 backprop phase 的 gradient。 我們已經獲取 softmax backward phase 的輸入 gradient：$\frac{\partial L}{\partial out_s}$。由於只有一個是有值的，其他都是 0，因此我們可以忽略除了 $out_s(c)$ 之外的其他值！

　　首先，讓我們計算 $out_s(c)$ 對於 totals （上面代碼中的，softmax 轉換前的值）的gradient。讓 $t_i$ 來表示 total 的類 $i$。然后我們可以把 $out_s(c)$ 寫作：

$$out_s(c)=\frac{e^{t_c}}{\sum_{i}e^{t_i}}=\frac{e^{t_c}}{S}$$

　　其中，$S=\sum_{i}e^{t_i}$。

　　現在，開始考慮一些類 $k$，其中 $k\neq c$。我們可以把 $out_s(c)$ 寫作：（由於只有 $out_s(c)$ 有值，因此只需考慮它就行了，其中 $e^{t_c}$ 相當於常數不用考慮）

$$out_s(c)=e^{t_c}S^{-1}$$ 

　　使用 Chain Rule 得到：

 

$$
\begin{equation}
\begin{split}
\frac{\partial out_s(c)}{\partial t_k}&=\frac{\partial out_s(c)}{\partial S}(\frac{\partial S}{\partial t_k})\\
&=-e^{t_c}S^{-2}(\frac{\partial S}{\partial t_k})\\
&=-e^{t_c}S^{-2}(e^{t_k})\\
&=\frac{-e^{t_c}e^{t_k}}{S^2}
\end{split}
\end{equation}
$$

　　上面是針對 $k\neq c$。現在讓我們算下 $k=c$ 的時候，如下所示：

$$
\begin{equation}
\begin{split}
\frac{\partial out_s(c)}{\partial t_c}&=\frac{Se^{t_c}-e^{t_c}\frac{\partial S}{\partial t_c}}{S^2}\\
&=\frac{Se^{t_c}-e^{t_c}e^{t_c}}{S^2}\\
&=\frac{e^{t_c}(S-e^{t_c})}{S^2}
\end{split}
\end{equation}
$$

　　合並如下：

$$
\begin{equation}
\frac{\partial out_s(k)}{\partial t}=
\begin{cases}
\frac{-e^{t_c}e^{t_k}}{S^2}& \text{if } k\neq c\\
\frac{e^{t_c}(S-e^{t_c})}{S^2}& \text{if } k=c
\end{cases}
\end{equation}
$$

　　如下實現：

class Softmax:
    # ...

    def backprop(self, d_L_d_out):
        '''
        Performs a backward pass of the softmax layer.
        Returns the loss gradient for this layer's inputs.
        - d_L_d_out is the loss gradient for this layer's outputs.
        '''
        # We know only 1 element of d_L_d_out will be nonzero
        for i, gradient in enumerate(d_L_d_out):
            # 找到 label 的值，就是 gradient 不為 0 的
            if gradient == 0:
                continue

            # e^totals
            t_exp = np.exp(self.last_totals)

            # Sum of all e^totals
            S = np.sum(t_exp)

            # Gradients of out[i] against totals
            # 初始化都設置為 非 c 的值，再單獨修改 c 的值
            d_out_d_t = -t_exp[i] * t_exp / (S ** 2)
            d_out_d_t[i] = t_exp[i] * (S - t_exp[i]) / (S ** 2)

            # ... to be continued

　　 我們繼續哈。我們最終是想要計算 loss 對於 weights，biases 和 input 的 gradient：

• 我們要使用 weights gradient，$\frac{\partial L}{\partial w}$，來更新層的 weights。
• 我們要使用 biases gradient，$\frac{\partial L}{\partial b}$，來更新層的 biases。
• 我們要返回 input（每一層的正向輸入） 的 gradient，$\frac{\partial L}{\partial input}$，基於 backprop 的方法，所以下一層可以使用它。

　　為了計算上面 3 個 loss gradient，我們首先需要獲取另外 3 個結果：totals（做 softmax 之前的向量，10 個元素）對於 weights，biases 和 input 的 gradient。相關公式如下：（以下為對於單獨 weight 的計算，但是代碼實現的時候是通過 matrix，相對抽象）

$$t=w*input+b$$

　　這些 gradient 很容易計算：

$$\frac{\partial t}{\partial w}=input$$

$$\frac{\partial t}{\partial b}=1$$

$$\frac{\partial t}{\partial input} = w$$

　　根據 Chain Rule 把它們放在一起：

$$\frac{\partial L}{\partial w}=\frac{\partial L}{\partial out}*\frac{\partial out}{\partial t}*\frac{\partial t}{\partial w}$$

$$\frac{\partial L}{\partial b}=\frac{\partial L}{\partial out}*\frac{\partial out}{\partial t}*\frac{\partial t}{\partial b}$$

$$\frac{\partial L}{\partial input}=\frac{\partial L}{\partial out}*\frac{\partial out}{\partial t}*\frac{\partial t}{\partial input}$$

　　其中，

• $L$：loss 函數
• $out$：做 softmax 的輸出結果，與 loss 公式直接相關的 概率
• $t$：做 softmax 的輸入參數，通過 weights，bias 以及 softmax 層的輸入來獲取

　　把它們一並放到代碼中實現如下：

class Softmax:
    # ...

    def backprop(self, d_L_d_out):
        '''
        Performs a backward pass of the softmax layer.
        Returns the loss gradient for this layer's inputs.
        - d_L_d_out is the loss gradient for this layer's outputs.
        '''
        # We know only 1 element of d_L_d_out will be nonzero
        for i, gradient in enumerate(d_L_d_out):
            if gradient == 0:
                continue

            # e^totals
            t_exp = np.exp(self.last_totals)

            # Sum of all e^totals
            S = np.sum(t_exp)

            # Gradients of out[i] against totals
            d_out_d_t = -t_exp[i] * t_exp / (S ** 2)
            d_out_d_t[i] = t_exp[i] * (S - t_exp[i]) / (S ** 2)

            # NEW ADD
            # Gradients of totals against weights/biases/input
            # d_t_d_w 的結果是 softmax 層的輸入數據，1352 個元素的向量
            # 不是最終的結果，最終結果是 2d 矩陣，1352x10
            d_t_d_w = self.last_input
            d_t_d_b = 1
            # d_t_d_input 的結果是 weights 值，2d 矩陣，1352x10
            d_t_d_inputs = self.weights

            # Gradients of loss against totals
            # 向量，10
            d_L_d_t = gradient * d_out_d_t

            # Gradients of loss against weights/biases/input
            # np.newaxis 可以幫助一維向量變成二維矩陣
            # (1352, 1) @ (1, 10) to (1352, 10)
            d_L_d_w = d_t_d_w[np.newaxis].T @ d_L_d_t[np.newaxis]
            d_L_d_b = d_L_d_t * d_t_d_b
            # (1352, 10) @ (10, 1) to (1352, 1)
            d_L_d_inputs = d_t_d_inputs @ d_L_d_t
            # ... to be continued

　　 計算出 gradient 之后，剩下的就是訓練 softmax 層。我們通過 SGD（Stochastic Gradient Decent）來更新 weights 和 bias，並返回 d_L_d_inputs：

class Softmax
    # ...

    # ADD A NEW PARAMETER - learn_rate
    def backprop(self, d_L_d_out, learn_rate):
        '''
        Performs a backward pass of the softmax layer.
        Returns the loss gradient for this layer's inputs.
        - d_L_d_out is the loss gradient for this layer's outputs.
        - learn_rate is a float
        '''
        # We know only 1 element of d_L_d_out will be nonzero
        for i, gradient in enumerate(d_L_d_out):
            if gradient == 0:
                continue

            # e^totals
            t_exp = np.exp(self.last_totals)

            # Sum of all e^totals
            S = np.sum(t_exp)

            # Gradients of out[i] against totals
            d_out_d_t = -t_exp[i] * t_exp / (S ** 2)
            d_out_d_t[i] = t_exp[i] * (S - t_exp[i]) / (S ** 2)

            # Gradients of totals against weights/biases/input
            d_t_d_w = self.last_input
            d_t_d_b = 1
            d_t_d_inputs = self.weights

            # Gradients of loss against totals
            d_L_d_t = gradient * d_out_d_t

            # Gradients of loss against weights/biases/input
            d_L_d_w = d_t_d_w[np.newaxis].T @ d_L_d_t[np.newaxis]
            d_L_d_b = d_L_d_t * d_t_d_b
            d_L_d_inputs = d_t_d_inputs @ d_L_d_t

            # NEW ADD
            # Update weights / biases
            self.weights -= learn_rate * d_L_d_w
            self.biases -= learn_rate * d_L_d_b
            # 將矩陣從 1d 轉為 3d
            # 1352 to 13x13x8
            return d_L_d_inputs.reshape(self.last_input_shape)

　　注意我們添加了 learn_rate 參數用來控制更新 weights 與 biases 的快慢。此外，我們需要將 d_L_d_inputs 進行 reshape() 操作，因為我們在 forward pass 中將 input 進行了 flatten() 操作。reshape() 操作之后，保證與原始輸入具有相同的結構。

8. 反向傳播：池化層（Backprop: Max Pooling）

　　池化層不需要訓練，因為它里面不存在任何 weights，但是為了計算 gradient 我們仍然需要實現一個 backprop() 方法。首先我們還是需要存儲一些臨時數據在 forward phase 里面。我們這次需要存儲的是 input。

class MaxPool2:
    # ...

    def forward(self, input):
        '''
        Performs a forward pass of the maxpool layer using the given input.
        Returns a 3d numpy array with dimensions (h / 2, w / 2, num_filters).
        - input is a 3d numpy array with dimensions (h, w, num_filters)
        '''
        # 存儲 池化層 的輸入參數，26x26x8
        self.last_input = input
        # More implementation
        # ...

　　在 forward pass 的過程中，Max Pooling 層選取 2x2 塊的最大值進行輸入，如下圖所示：

　　

　　backward phase 中的相同層如下圖所示：

　　 

　　每一個 gradient 的值都被賦值到原始的最大值的位置，其他的值都是 0。

　　為什么 backward phase 的 Max Pooling 層顯示如上呢？讓我們直覺思考下 $\frac{\partial L}{\partial inputs}$ （Max Pooling 的輸入數據，26x26x8）的值是多少。對於 2x2 數據塊中不是最大值的輸入像素將不會對 loss 產生任何影響，因為稍微改變這個值並不會改變輸出！換句話說，對於非最大值的像素點：$\frac{\partial L}{\partial input}=0$。另一方面，最大值的像素點會將值傳遞給輸出，所以 $\frac{\partial output}{\partial input}=1$，也就是說，$\frac{\partial L}{\partial input}=\frac{\partial output}{\partial input}$。

　　總結后就是：（output 與 input 都是相對於 Max Pooling 層來說的）

$$
\begin{equation}
\frac{\partial L}{\partial input}=
\begin{cases}
0& \text{if } input \neq max\\
\frac{\partial L}{\partial output}& \text{if } input=max
\end{cases}
\end{equation}
$$

　　代碼實現如下：

class MaxPool2:
    # ...

    def iterate_regions(self, image):
        '''
        Generates non-overlapping 2x2 image regions to pool over.
        - image is a 2d numpy array
        '''
        h, w, _ = image.shape
        new_h = h // 2
        new_w = w // 2

        for i in range(new_h):
            for j in range(new_w):
                im_region = image[(i * 2):(i * 2 + 2), (j * 2):(j * 2 + 2)]
                yield im_region, i, j

    def backprop(self, d_L_d_out):
        '''
        Performs a backward pass of the maxpool layer.
        Returns the loss gradient for this layer's inputs.
        - d_L_d_out is the loss gradient for this layer's outputs.
        '''
        # 池化層輸入數據，26x26x8，默認初始化為 0
        d_L_d_input = np.zeros(self.last_input.shape)

        # 每一個 im_region 都是一個 3x3x8 的8層小矩陣
        # 修改 max 的部分，首先查找 max
        for im_region, i, j in self.iterate_regions(self.last_input):
            h, w, f = im_region.shape
            # 獲取 im_region 里面最大值的索引向量，一疊的感覺
            amax = np.amax(im_region, axis=(0, 1))

            # 遍歷整個 im_region，對於傳遞下去的像素點，修改 gradient 為 loss 對 output 的gradient
            for i2 in range(h):
                for j2 in range(w):
                    for f2 in range(f):
                        # If this pixel was the max value, copy the gradient to it.
                        if im_region[i2, j2, f2] == amax[f2]:
                            d_L_d_input[i * 2 + i2, j * 2 + j2, f2] = d_L_d_out[i, j, f2]

        return d_L_d_input

　　對於每一個 2x2 的像素塊，我們找到 forward pass 中最大值的像素點，然后將 loss 對 output 的 gradient 復制過去 。

　　就是醬紫來弄，接下來是最后一層了。

9. 反向傳播：卷積層（Backprop: Conv）

　　終於到卷基層了：卷積層的反向傳播是 CNN 模型訓練的核心。forward phase 存儲很簡單：

class Conv3x3
    # ...

    def forward(self, input):
        '''
        Performs a forward pass of the conv layer using the given input.
        Returns a 3d numpy array with dimensions (h, w, num_filters).
        - input is a 2d numpy array
        '''
        # 輸入大數據，28x28
        self.last_input = input
        # More implementation
        # ...

　　我們主要是對卷基層的 filter 感興趣，因為我們需要跟新 filter 的 weight。我們已經得到了卷積層的 $\frac{\partial L}{\partial out}$，所以我們需要獲取 $\frac{\partial out}{\partial filters}$。為了計算這個值，我們需要問下自己：怎么樣改變 filter 的 weight 來影響 卷積層 的輸出的？

　　實際上，改變任何 filter 的 weight 都會影響到整個輸出圖片的信息，因為在卷積過程中，每一個輸出的像素都會使用每一個 filter 的 weight。為了簡單起見，我們試想下一次只有一個輸出：如何修改 filter 來改變那個具體輸出像素的值？

　　下面這個例子有助於我們思考這個問題：

　　

　　我們有一個 3x3 的圖片與一個都是 0 的 3x3 的 filter 進行卷積運算，結果只有一個 1x1 的輸出。如果我們把 filter 中間的 weight 增加到 1 呢？輸出將會隨着中心值來增加到 80：

　　

　　簡單起見，增加任何 filter 的其他權重到 1，都會最終增加相應的輸出圖片像素值！這說明一個具體的輸出像素對於具體的 filter 的 weight 的 gradient 就是對應的像素值。推導如下：

\begin{equation}
\begin{split}
out(i, j)&=convolve(image, filter)\\
&=\sum_{x=0}^{3}\sum_{y=0}^{3}image(i+x, j+y)*filter(x, y)
\end{split}
\end{equation}

$$\frac{\partial out(i, j)}{\partial filter(x, y)}=image(i+x, j+y)$$

　　如下圖所示，對於任意一個 $out(i, j)$ 都是通過 image 中的 3x3 矩陣 與 filter 的 3x3 矩陣進行點乘求和獲取的，因此對於 任意一個 $out(i, j)$ 對於 任意一個 $filter(x, y)$ 的 gradient 就是與其對應相乘的那個像素點 $image(i+x, j+y)$。

　　

　　於是，我們可以實現卷積層的 backprop 如下：

class Conv3x3
    # ...

    def backprop(self, d_L_d_out, learn_rate):
        '''
        Performs a backward pass of the conv layer.
        - d_L_d_out is the loss gradient for this layer's outputs.
        - learn_rate is a float.
        '''
        # 初始化一組為 0 的 gradient，3x3x8
        d_L_d_filters = np.zeros(self.filters.shape)

        # im_region，一個個 3x3 小矩陣
        for im_region, i, j in self.iterate_regions(self.last_input):
            for f in range(self.num_filters):
                # 按 f 分層計算，一次算一層，然后累加起來
                # d_L_d_filters[f]: 3x3 matrix
                # d_L_d_out[i, j, f]: num
                # im_region: 3x3 matrix in image
                d_L_d_filters[f] += d_L_d_out[i, j, f] * im_region

        # Update filters
        self.filters -= learn_rate * d_L_d_filters

        # We aren't returning anything here since we use Conv3x3 as
        # the first layer in our CNN. Otherwise, we'd need to return
        # the loss gradient for this layer's inputs, just like every
        # other layer in our CNN.
        return None

　　至此，我們已經實現了 CNN 的整個 backward pass。接下來我們來測試下...

　　完整代碼參考：CNN from scratch - github

10. 訓練 CNN（Training a CNN）

　　我們將要訓練我們的 CNN 模型通過幾個 epoch，跟蹤訓練中的改進，並且在另外的測試集上進行測試。下面是完整的代碼：

import mnist
import numpy as np

# We only use the first 1k examples of each set in the interest of time.
# Feel free to change this if you want.
train_images = mnist.train_images()[:1000]
train_labels = mnist.train_labels()[:1000]
test_images = mnist.test_images()[:1000]
test_labels = mnist.test_labels()[:1000]

conv = Conv3x3(8)                                    # 28x28x1 -> 26x26x8
pool = MaxPool2()                                    # 26x26x8 -> 13x13x8
softmax = Softmax(13 * 13 * 8, 10) # 13x13x8 -> 10

def forward(image, label):
    '''
    Completes a forward pass of the CNN and calculates the accuracy and
    cross-entropy loss.
    - image is a 2d numpy array
    - label is a digit
    '''
    # We transform the image from [0, 255] to [-0.5, 0.5] to make it easier
    # to work with. This is standard practice.
    out = conv.forward((image / 255) - 0.5)
    out = pool.forward(out)
    out = softmax.forward(out)

    # Calculate cross-entropy loss and accuracy. np.log() is the natural log.
    loss = -np.log(out[label])
    acc = 1 if np.argmax(out) == label else 0

    return out, loss, acc
    # out: vertor of probability
    # loss: num
    # acc: 1 or 0

def train(im, label, lr=.005):
    '''
    Completes a full training step on the given image and label.
    Returns the cross-entropy loss and accuracy.
    - image is a 2d numpy array
    - label is a digit
    - lr is the learning rate
    '''
    # Forward
    out, loss, acc = forward(im, label)

    # Calculate initial gradient
    gradient = np.zeros(10)
    gradient[label] = -1 / out[label]

    # Backprop
    gradient = softmax.backprop(gradient, lr)
    gradient = pool.backprop(gradient)
    gradient = conv.backprop(gradient, lr)

    return loss, acc

print('MNIST CNN initialized!')

# Train the CNN for 3 epochs
for epoch in range(3):
    print('--- Epoch %d ---' % (epoch + 1))

    # Shuffle the training data
    permutation = np.random.permutation(len(train_images))
    train_images = train_images[permutation]
    train_labels = train_labels[permutation]

    # Train!
    loss = 0
    num_correct = 0

    # i: index
    # im: image
    # label: label
    for i, (im, label) in enumerate(zip(train_images, train_labels)):
        if i > 0 and i % 100 == 99:
            print(
                '[Step %d] Past 100 steps: Average Loss %.3f | Accuracy: %d%%' %
                (i + 1, loss / 100, num_correct)
            )
            loss = 0
            num_correct = 0

        l, acc = train(im, label)
        loss += l
        num_correct += acc

# Test the CNN
print('\n--- Testing the CNN ---')
loss = 0
num_correct = 0
for im, label in zip(test_images, test_labels):
    _, l, acc = forward(im, label)
    loss += l
    num_correct += acc

num_tests = len(test_images)
print('Test Loss:', loss / num_tests)
print('Test Accuracy:', num_correct / num_tests)

　　例子的輸出結果如下：

MNIST CNN initialized!
--- Epoch 1 ---
[Step 100] Past 100 steps: Average Loss 2.254 | Accuracy: 18%
[Step 200] Past 100 steps: Average Loss 2.167 | Accuracy: 30%
[Step 300] Past 100 steps: Average Loss 1.676 | Accuracy: 52%
[Step 400] Past 100 steps: Average Loss 1.212 | Accuracy: 63%
[Step 500] Past 100 steps: Average Loss 0.949 | Accuracy: 72%
[Step 600] Past 100 steps: Average Loss 0.848 | Accuracy: 74%
[Step 700] Past 100 steps: Average Loss 0.954 | Accuracy: 68%
[Step 800] Past 100 steps: Average Loss 0.671 | Accuracy: 81%
[Step 900] Past 100 steps: Average Loss 0.923 | Accuracy: 67%
[Step 1000] Past 100 steps: Average Loss 0.571 | Accuracy: 83%
--- Epoch 2 ---
[Step 100] Past 100 steps: Average Loss 0.447 | Accuracy: 89%
[Step 200] Past 100 steps: Average Loss 0.401 | Accuracy: 86%
[Step 300] Past 100 steps: Average Loss 0.608 | Accuracy: 81%
[Step 400] Past 100 steps: Average Loss 0.511 | Accuracy: 83%
[Step 500] Past 100 steps: Average Loss 0.584 | Accuracy: 89%
[Step 600] Past 100 steps: Average Loss 0.782 | Accuracy: 72%
[Step 700] Past 100 steps: Average Loss 0.397 | Accuracy: 84%
[Step 800] Past 100 steps: Average Loss 0.560 | Accuracy: 80%
[Step 900] Past 100 steps: Average Loss 0.356 | Accuracy: 92%
[Step 1000] Past 100 steps: Average Loss 0.576 | Accuracy: 85%
--- Epoch 3 ---
[Step 100] Past 100 steps: Average Loss 0.367 | Accuracy: 89%
[Step 200] Past 100 steps: Average Loss 0.370 | Accuracy: 89%
[Step 300] Past 100 steps: Average Loss 0.464 | Accuracy: 84%
[Step 400] Past 100 steps: Average Loss 0.254 | Accuracy: 95%
[Step 500] Past 100 steps: Average Loss 0.366 | Accuracy: 89%
[Step 600] Past 100 steps: Average Loss 0.493 | Accuracy: 89%
[Step 700] Past 100 steps: Average Loss 0.390 | Accuracy: 91%
[Step 800] Past 100 steps: Average Loss 0.459 | Accuracy: 87%
[Step 900] Past 100 steps: Average Loss 0.316 | Accuracy: 92%
[Step 1000] Past 100 steps: Average Loss 0.460 | Accuracy: 87%

--- Testing the CNN ---
Test Loss: 0.5979384893783474
Test Accuracy: 0.78

　　我們的代碼效果不錯，實現了 78% 的准確率。

11. Keras 實現

　　通過 Keras 實現上面的功能如下：

import numpy as np
import mnist
from keras.models import Sequential
from keras.layers import Conv2D, MaxPooling2D, Dense, Flatten
from keras.utils import to_categorical
from keras.optimizers import SGD

train_images = mnist.train_images()
train_labels = mnist.train_labels()
test_images = mnist.test_images()
test_labels = mnist.test_labels()

train_images = (train_images / 255) - 0.5
test_images = (test_images / 255) - 0.5

train_images = np.expand_dims(train_images, axis=3)
test_images = np.expand_dims(test_images, axis=3)

model = Sequential([
    Conv2D(8, 3, input_shape=(28, 28, 1), use_bias=False),
    MaxPooling2D(pool_size=2),
    Flatten(),
    Dense(10, activation='softmax'),
])

model.compile(SGD(lr=.005), loss='categorical_crossentropy', metrics=['accuracy'])

model.fit(
    train_images,
    to_categorical(train_labels),
    batch_size=1,
    epochs=3,
    validation_data=(test_images, to_categorical(test_labels)),
)

　　以上代碼應用了 MNIST 的全部數據集，結果如下：

Epoch 1
loss: 0.2433 - acc: 0.9276 - val_loss: 0.1176 - val_acc: 0.9634
Epoch 2
loss: 0.1184 - acc: 0.9648 - val_loss: 0.0936 - val_acc: 0.9721
Epoch 3
loss: 0.0930 - acc: 0.9721 - val_loss: 0.0778 - val_acc: 0.9744

　　得到 97.4% 的准確率！