Digit Recognizer

手写数字识别是非常经典的入门级实践项目了，MNIST则是用来初步检验一个图像分类方法好坏的不二之选数据集。之前练习keras的时候也用CNN跑过正确率90%以上的模型，但在kaggle这个大神云集的地方，这样的score显然是不够的。

这个竞赛给我带来的收获主要有以下几点：

1、熟悉了keras/pandas/numpy等库的操作

2、学习了如何初步优化所搭建的网络结构

3、第一次实际使用数据增强来辅助图片识别任务，提升精度

下面我们一步步来看一下这个问题是如何求解的，思路主要参考了guide-to-cnns-with-data-augmentation-keras以及how-to-choose-cnn-architecture-mnist两篇kernel，令我受益匪浅。

首先导入我们需要用到的库和包：

import keras
import numpy as np
import pandas as pd

from sklearn.model_selection import train_test_split
from keras.utils import to_categorical
from keras.models import Sequential
from keras.layers.convolutional import Conv2D, MaxPooling2D
from keras.layers import BatchNormalization, Dense, Dropout, Flatten
from keras.preprocessing.image import ImageDataGenerator
from keras.callbacks import LearningRateScheduler

这里导入的时候最好按照实际使用的顺序进行导入，这样方便我们理顺思路，而且当后面有新的方法要尝试的时候也方便查找所需要的库和包是否已经导入。

首先，老规矩，先对数据进行读入和展示：

train = pd.read_csv(r'C:\Users\Sidney Nash\Desktop\data\Digits\train.csv')
test = pd.read_csv(r'C:\Users\Sidney Nash\Desktop\data\Digits\test.csv')

train.head()

可以看到数据是以像素值作为特征给出的，其中label是图片对应的数字标签。

我们将数据分解为 feature 和 label ，分别用 X 和 y 表示：

y = train["label"]
train.drop(["label"], inplace=True, axis=1) #删去label列
X = train

此时 X 和 y 表示的是训练数据集中的数据特征及标签，为了检验模型训练效果，我们需要从训练集中随机抽取一部分作为验证集。

X_train, X_val, y_train, y_val = train_test_split(X, y, train_size=0.9, random_state=42)

这里我们设置train_size = 0.9意为使用10%的训练集数据作为验证集。

接下来我们要对数据格式进行处理，首先为了使用 CNN 我们要把 784 维的向量还原成28281 的矩阵，这里的 1 代表 channel 数。

在此基础上我们对数据进行正规化，使其范围在。

X_train = X_train.reshape(-1, 28, 28, 1).astype("float32") / 255
X_val = X_val.reshape(-1, 28, 28, 1).astype("float32") / 255

然后把类别标签 y 转化成 one-hot 向量来表示。

y_train[0:3]

out:
22460    8
20828    5
32032    6
Name: label, dtype: int64

y_train = to_categorical(y_train, num_classes=10)
y_val = to_categorical(y_val, num_classes=10)

y_train[0:3]

out:
array([[0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.]], dtype=float32)

处理好数据，接下来就可以搭建模型了。我们的网络主要由 conv2D、MaxPooling2D、BatchNormalization、Dropout、Dense组成，但是凭空搭建一个网络，我们并不知道这个架构是否合理，how-to-choose-cnn-architecture-mnist一文提供了一种思路，虽然需要付出一些计算代价来做实验，且面临各个参数局部最优组合起来未必是全局最优的问题，但比起凭空搭建一个网络还是要靠谱许多的。

首先我们要确定用几个卷积层+池化层的组合是合适的，这一点其实很重要，因为这直接关系到提取特征的质量好坏。我们可以以 55 的卷积核为例，试验一下使用1、2、3个卷积层+池化层的效果。

nets = 3
model = [0] *nets

for j in range(3):
    model[j] = Sequential()
    model[j].add(Conv2D(24,kernel_size=5,padding='same',activation='relu',
            input_shape=(28,28,1)))
    model[j].add(MaxPooling2D())
    if j>0:
        model[j].add(Conv2D(48,kernel_size=5,padding='same',activation='relu'))
        model[j].add(MaxPooling2D())
    if j>1:
        model[j].add(Conv2D(64,kernel_size=5,padding='same',activation='relu'))
        model[j].add(MaxPooling2D(padding='same'))
    model[j].add(Flatten())
    model[j].add(Dense(256, activation='relu'))
    model[j].add(Dense(10, activation='softmax'))
    model[j].compile(optimizer="adam", loss="categorical_crossentropy", metrics=["accuracy"])

history = [0] * nets
names = ["(C-P)x1","(C-P)x2","(C-P)x3"]
epochs = 20
annealer = LearningRateScheduler(lambda x: 1e-3 * 0.9 ** x)

for j in range(nets):
    history[j] = model[j].fit(X_train,y_train, batch_size=80, epochs = epochs, 
                              validation_data = (X_val, y_val), callbacks=[annealer])
    print("CNN {0}: Epochs={1:d}, Train accuracy={2:.5f}, Validation accuracy={3:.5f}".format(
    names[j],epochs,max(history[j].history['acc']),max(history[j].history['val_acc'])))

out:
CNN (C-P)x1: Epochs=20, Train accuracy=0.99997, Validation accuracy=0.98833
CNN (C-P)x2: Epochs=20, Train accuracy=0.99997, Validation accuracy=0.99190
CNN (C-P)x3: Epochs=20, Train accuracy=0.99997, Validation accuracy=0.99214

可以看到，1层的效果较差，2层3层效果相差不多，3层略好。考虑到计算开销，我们使用2层是比较合适的。

当然，如果只看最终结果不够直观，我们还可以画出整个训练过程中各model的精度变化曲线。

plt.figure(figsize=(15,5))
for i in range(nets):
    plt.plot(history[i].history['val_acc'])
plt.title('model accuracy')
plt.ylabel('accuracy')
plt.xlabel('epoch')
plt.legend(names, loc='upper left')
axes = plt.gca()
axes.set_ylim([0.98,1])
plt.show()

类似上述过程，我们还可以做实验寻找合适的 feature maps 数量。由于算力有限，这里直接引用作者的结果了。

nets = 6
model = [0] *nets
for j in range(6):
    model[j] = Sequential()
    model[j].add(Conv2D(j*8+8,kernel_size=5,activation='relu',input_shape=(28,28,1)))
    model[j].add(MaxPool2D())
    model[j].add(Conv2D(j*16+16,kernel_size=5,activation='relu'))
    model[j].add(MaxPool2D())
    model[j].add(Flatten())
    model[j].add(Dense(256, activation='relu'))
    model[j].add(Dense(10, activation='softmax'))
    model[j].compile(optimizer="adam", loss="categorical_crossentropy", metrics=["accuracy"])

out:
CNN 8 maps: Epochs=20, Train accuracy=0.99979, Validation accuracy=0.98820
CNN 16 maps: Epochs=20, Train accuracy=0.99939, Validation accuracy=0.99056
CNN 24 maps: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99178
CNN 32 maps: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99199
CNN 40 maps: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99256
CNN 48 maps: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99278

综合考虑精度和计算代价，选用第一层 32 个 maps 第二层 64 个 maps 的架构比较合理。

设定好卷积层和池化层，接下来我们可以对全连接层进行实验了。首先作者实验了多 Dense 隐层的情形，发现性能相比单 Dense 隐层并没有显著提升。接下来就是寻找 Dense 层合适的节点数了。

nets = 8
model = [0] *nets

for j in range(8):
    model[j] = Sequential()
    model[j].add(Conv2D(32,kernel_size=5,activation='relu',input_shape=(28,28,1)))
    model[j].add(MaxPool2D())
    model[j].add(Conv2D(64,kernel_size=5,activation='relu'))
    model[j].add(MaxPool2D())
    model[j].add(Flatten())
    if j>0:
        model[j].add(Dense(2**(j+4), activation='relu'))
    model[j].add(Dense(10, activation='softmax'))
    model[j].compile(optimizer="adam", loss="categorical_crossentropy", metrics=["accuracy"])

out:
CNN 0N: Epochs=20, Train accuracy=0.99975, Validation accuracy=0.98885
CNN 32N: Epochs=20, Train accuracy=0.99979, Validation accuracy=0.98927
CNN 64N: Epochs=20, Train accuracy=0.99993, Validation accuracy=0.99056
CNN 128N: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99063
CNN 256N: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99092
CNN 512N: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99121
CNN 1024N: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99178
CNN 2048N: Epochs=20, Train accuracy=0.99996, Validation accuracy=0.99149

综合模型在训练集和验证集上的表现，并权衡计算代价，可以看出 128 个结点是比较合理的选择。

最后，我们还可以实验一下 Dropout 层 drop 多少比例的结点是合适的。

ets = 8
model = [0] *nets

for j in range(8):
    model[j] = Sequential()
    model[j].add(Conv2D(32,kernel_size=5,activation='relu',input_shape=(28,28,1)))
    model[j].add(MaxPool2D())
    model[j].add(Dropout(j*0.1))
    model[j].add(Conv2D(64,kernel_size=5,activation='relu'))
    model[j].add(MaxPool2D())
    model[j].add(Dropout(j*0.1))
    model[j].add(Flatten())
    model[j].add(Dense(128, activation='relu'))
    model[j].add(Dropout(j*0.1))
    model[j].add(Dense(10, activation='softmax'))
    model[j].compile(optimizer="adam", loss="categorical_crossentropy", metrics=["accuracy"])

out:
CNN D=0: Epochs=30, Train accuracy=1.00000, Validation accuracy=0.99099
CNN D=0.1: Epochs=30, Train accuracy=0.99961, Validation accuracy=0.99256
CNN D=0.2: Epochs=30, Train accuracy=0.99889, Validation accuracy=0.99271
CNN D=0.3: Epochs=30, Train accuracy=0.99732, Validation accuracy=0.99299
CNN D=0.4: Epochs=30, Train accuracy=0.99411, Validation accuracy=0.99314
CNN D=0.5: Epochs=30, Train accuracy=0.98929, Validation accuracy=0.99328
CNN D=0.6: Epochs=30, Train accuracy=0.98012, Validation accuracy=0.99206
CNN D=0.7: Epochs=30, Train accuracy=0.96441, Validation accuracy=0.98956

综合模型在训练集和验证集上的表现，可以看出 40% 是比较合理的选择。

至此，我们的模型架构基本成型了，不过，在实际搭建网络的时候，作者只是以这些信息为指导，然后根据经验来设计网络，比如作者提到：

• 可以用两个 33 的卷积层来代替一个 55 的卷积层，在感受野相似的情况下增加模型的非线性。

• 可以用 strides=2 的卷积层来代替 MaxPooling 层。至于原因，我认为知乎上的一个回答提供了一种合理的解释：1、conv 可以看做是 pooing 的超集，mean pooling 是权重为 1/num 的 conv ，max pooling 是只在最大值的位置权值是 1 的conv。所以理论上任何 pooling 层都可以换成对应的 conv 层；2、用 pooling 层其实这是人为添加的比较强的先验，比如第一层卷积之后加 max-pooling 可以认为是为了降低一些噪声干扰，在最后一层添加 mean-pooling 则可以认为是对 feature 的加权。因此我们使用 conv 代替 pooling 的一个好处就是，可以让模型自己学习什么时候应该用何种形式的 pooling 。

• 为了提升训练速度及稳定性，在层与层之间加入 BatchNormalization 层。

基于上述理由，作者最终给出的网络架构如下：

model = Sequential()

#用两个 3*3 的卷积层来代替一个 5*5 的卷积层
model.add(Conv2D(32,kernel_size=3,activation='relu',input_shape=(28,28,1)))
model.add(BatchNormalization())
model.add(Conv2D(32,kernel_size=3,activation='relu'))
model.add(BatchNormalization())
#用一个conv(strides=2)来代替MaxPooling
model.add(Conv2D(32,kernel_size=5,strides=2,padding='same',activation='relu'))
model.add(BatchNormalization())
model.add(Dropout(0.4))

model.add(Conv2D(64,kernel_size=3,activation='relu'))
model.add(BatchNormalization())
model.add(Conv2D(64,kernel_size=3,activation='relu'))
model.add(BatchNormalization())
model.add(Conv2D(64,kernel_size=5,strides=2,padding='same',activation='relu'))
model.add(BatchNormalization())
model.add(Dropout(0.4))

model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(BatchNormalization())
model.add(Dropout(0.4))
model.add(Dense(10, activation='softmax'))

model.compile(optimizer="adam", loss="categorical_crossentropy", metrics=["accuracy"])

用上述模型直接在训练集上进行训练，已经可以得到大概 99.5% 的 score 。这是很不错的表现了，在此基础上，我们可以使用 Data Augmentation 的方法使得 score 提升至 99.6%~99.7% 。使用 Ensemble 的方法同样可以达到差不多的 score ，但对神经网络的 Ensemble 需要的计算代价太大，这里不予考虑。

数据增强是如何进行的呢？我们通过下面的可视化过程来看一下。首先我们选取一个训练样本作为例子。

data = X_train[0].reshape(28,28)
plt.imshow(data)

接下来我们对这张手写数字 8 的图片进行数据增强。

datagen = ImageDataGenerator(
                            rotation_range=40,
                            width_shift_range=0.2,
                            height_shift_range=0.2,
                            shear_range=0.2,
                            zoom_range=0.2,
                            horizontal_flip=True,
                            fill_mode='nearest')

plt.figure(figsize=(8,8))
i=0
for batch in datagen.flow(data, batch_size=1):
    plt.subplot(8,8,i+1)
    plt.axis("off")
    imgplot = plt.imshow(batch[0].reshape(28,28))
    i += 1
    if i == 64:
        break
    
plt.show()

可以看到，我们可以对一个样本进行旋转、平移、缩放等操作生成新的增强样本供模型学习，这样训练出的模型对数字的位置以及方向等与分类无关的特征更加稳健，可以学到真正对分类有用的特征。

当然在进行数据增强时要注意参数的设定，比如最大旋转角度要合理，如果设置为90那么数字8就可能变成无穷符号，反而会对分类结果产生不好的影响。

最终模型的训练过程如下：

datagen = ImageDataGenerator(rotation_range=10,
                             width_shift_range=0.1,
                             height_shift_range=0.1)

datagen.fit(X_train)

annealer = LearningRateScheduler(lambda x: 1e-3 * 0.9 ** x)

model.fit_generator(datagen.flow(X_train, y_train, batch_size=64),
                    epochs = 50, validation_data = (X_val, y_val),
                    verbose = 1, steps_per_epoch=X_train.shape[0] / 64,
                    callbacks=[annealer])

最终模型在测试集上的表现达到了0.99657 ，我个人认为是比较满意的一个 score 了，或许多跑几次精度可以达到0.997，不过算力有限我就不尝试了。

最终成绩停留在 Top 11%，让我有点惊讶，看了眼 leaderboard 上前 93 名都达到了 100% 的精度，简直开挂……后来看了某大神 kernel 的揭秘才知道这些人是用了整个 MNIST 数据集而不是比赛中提供的 train_data 进行训练的，而整个 MNIST 数据集是包含 kaggle 上的test_data的，而且这些人直接用了 KNN ，简直有毒……不禁回忆起之前 Titanic 中直接利用遇难者名单进行匹配的操作……果然用来练手的比赛就是会有这样的问题啊……

总之这场比赛至此就告一段落了，虽然没拿到 top10% 但基本方法思路应该已经是较优的了。