Kaggle / MNIST をサポートベクターマシンで頑張る

はじめに

Kaggle / MNIST は, 28x28 の手書き数字の画像を学習するミッションである. 画像データの学習というと畳み込みニューラルネットワーク (Convolutional Neural Network, CNN) に手を出したくなるが, ここではあえてサポートベクターマシンでどこまで精度を上げられるかトライした.

精度的には (私の機械学習の実力では) CNN に劣るが, ネットで調べても Kaggle / MNIST をサポートベクターマシンでトライした事例があまり無かったので, 「これくらいの精度が出るよ」という記録として残しておく. 結果的には, 正答率は 0.98375 となった.

精度を上げるためにやってみたことは, 以下の 3 つである. これらについて説明する.

特徴量を増やす
役に立たない (と思われる) 特徴量を減らす
ハイパーパラメタを調整する

対象読者

kaggle を触ったことがある人 (データを読み込んで, 何らかの処理をして, submit したことがある人)
サポートベクターマシンについて, 何となく知っている人

データを読み込む

まず, 型通りにデータを train_data, test_data に読み込んで, データの長さを train_data_len. test_data_len に保存する.

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)

# Input data files are available in the read-only "../input/" directory
# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory

import os
for dirname, _, filenames in os.walk('/kaggle/input'):
    for filename in filenames:
        print(os.path.join(dirname, filename))

# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using "Save & Run All" 
# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session


# load data
train_data = pd.read_csv("/kaggle/input/digit-recognizer/train.csv")
test_data = pd.read_csv("/kaggle/input/digit-recognizer/test.csv")

train_data_len = len(train_data)
test_data_len = len(test_data)

print("Length of train_data ; {}".format(train_data_len))
print("Length of test_data ; {}".format(test_data_len))

# Length of train_data ; 42000
# Length of test_data ; 28000

続いて,

train_data_y ; ラベルを保存
train_data_x; 元のデータからラベルを drop したものを保存

train_data_y = train_data["label"]
train_data_x = train_data.drop(columns="label")

データの特徴量を増やす

数字の画像を眺める (生データを眺めるのは重要!) と, 例えば「1」と「8」では文字が書かれていピクセルの数 (面積) が大きく異なる. なので, 全ピクセル値の 平均値 は新たな特徴量になると考えられる. 同様に, 全ピクセル値の分散も新たな特徴量になると考えられる.

また, 28x28 のエリアの, 上半分・下半分のそれぞれの平均値と分散も特徴量になりうる. という感じで, 以下の領域で平均値と分散を求め, 新たな特徴量とする.

28x28 全体
上半分, 下半分
左半分, 右半分
1/4 のエリア
1/9 のエリア
1/16 のエリア

各エリアのピクセルを抽出するには pandas.DataFrame.query() を使うのが簡単だが, そのために元のデータの 行と列を入れ替え てからデータ処理する. またデータ処理のために, 通し番号を付けた no という列を追加する.

df = pd.concat([train_data_x, test_data]) # 訓練データとテストデータを一括で処理する
df_T = df.T                               # 行と列を入れ替え
df_T["no"] = range(len(df_T))             # no の列を追加

続いて, 28x28 全体の平均と分散を計算して, 新しい特徴量 (a_mean, a_std) を作る.

df_T.loc["a_mean"] = df_T.mean()
df_T.loc["a_std"] = df_T.std()

次は, 上半分・下半分の平均と分散を計算する. 文字列で抽出条件を作っている (便利!).

# horizontal

for i in range(2):
    q = 'no < 28*28/2' if i == 0 else 'no >= 28*28/2'
    df_T.loc["b{}_mean".format(i)] = df_T[:784].query(q).mean()
    df_T.loc["b{}_std".format(i)] = df_T[:784].query(q).std()

そして, 左半分・右半分の平均と分散を計算する. 抽出条件が少しずつややこしくなるが, no を 28 で割った余りが 14 未満か以上かで, 左半分と右半分を分けている.

# vertical

for i in range(2):
    q = 'no % 28 < 14' if i == 0 else 'no % 28 >= 14'
    df_T.loc["c{}_mean".format(i)] = df_T[:784].query(q).mean()
    df_T.loc["c{}_std".format(i)] = df_T[:784].query(q).std()

1/4, 1/9, 1/16 のエリアについても, 平均と分散を計算する.

# mean and std of 1/4 area
for i in range(2):
    qi = 'no < 28*28/2' if i == 0 else 'no >= 28*28/2'
    for j in range(2):
        qj = 'no % 28 < 14' if j == 0 else 'no % 28 >= 14'
        q = qi + " & " + qj
        num = i * 2 + j
        df_T.loc["d{}_mean".format(num)] = df_T[:784].query(q).mean()
        df_T.loc["d{}_std".format(num)] = df_T[:784].query(q).std()

# mean and std of 1/9 area
for i in range(3):
    if i == 0:
        qi = 'no < 262'
    elif i == 1:
        qi = "262 <= no < 522"
    else:
        qi = "522 <= no < 784"
        
    for j in range(3):
        if j == 0:
            qj = 'no % 28 < 9'
        elif j == 1:
            qj = '9 <= no % 28 < 18'
        else:
            qj = '18 <= no % 28'
        
        q = qi + " & " + qj
        num = i * 3 + j
        df_T.loc["e{}_mean".format(num)] = df_T[:784].query(q).mean()
        df_T.loc["e{}_std".format(num)] = df_T[:784].query(q).std()

# mean and std of 1/16 area
for i in range(4):
    qi = '{0} <= no < {1}'.format(28*28/4*i, 28*28/4*(i+1))
        
    for j in range(4):
        qj = '{0} <= no % 28 < {1}'.format(28/4*j, 28/4*(j+1))
        
        q = qi + " & " + qj
        num = i * 4 + j
        df_T.loc["f{}_mean".format(num)] = df_T[:784].query(q).mean()
        df_T.loc["f{}_std".format(num)] = df_T[:784].query(q).std()

最後に, no の列を drop して, 行と列をもとに戻す

df_T.drop(columns="no", inplace=True)
df = df_T.T

特徴量を減らす

数字の画像を見ていると, 例えば左上のピクセルは どの数字でも 0 なんじゃないか と思えてくる. すべての画像において値が同じピクセルは情報を持ってないと判断して, drop することにした (CNN の場合は, 周囲のピクセルとの関係が大事なので, こういうことはできないが).

各ピクセルの最大と最小を計算して, それが等しい列は drop することにした (それ以外にも, 例えば, 分散が小さい列を drop するなどの方法もあるかもしれない). 結果的に, 65 の列を drop した.

# drop columns if all values are the same
drop_col = []               # 空のリストを準備

for c in df.columns:        # 列ごとに
    col_max = df[c].max()
    col_min = df[c].min()
    if col_max == col_min:  # 最大と最小が等しければ drop する
        drop_col.append(c)
        
print("# of dropping columns ; {}".format(len(drop_col)))
        
df.drop(drop_col, axis=1, inplace=True)

# number of dropping columns ; 65

データの値を調整する

ハイパーパラメタを調整する前に, データの値を 0 ～ 1 に変換しておく. そうしないと, 後々の計算で長い時間がかかってしまう.

ここでは sklearn.preprocessing.MinMaxScaler を使った. これの transform() メソッドは ndarray を返すので, それを pandas.DataFrame に変換している.

# scaling
from sklearn import preprocessing
mmscaler = preprocessing.MinMaxScaler()
mmscaler.fit(df)

df_scaled = pd.DataFrame(mmscaler.transform(df), columns=df.columns, index=df.index)

# separate df_scaled into train_data and test data
train_data_x = df_scaled[:train_data_len]
test_data = df_scaled[train_data_len:]

ハイパーパラメタを調整する

サポートベクターマシンは sklearn.svm.SVC() を使う. ハイパーパラメタの調整には sklearn.model_selection.GridSearchCV() を使用し, 少量のデータを使いながら, 少しずつ C と gamma の範囲を狭めていった.

これら以外にも kernel や decision_function_shape などのハイパーパラメタもある. これらを振った計算もしたのだが, ほぼ下記の値が最適と出ていたので, 以降は下記の値を使用している.

kernel="rbf"
decision_function_shape="ovo"

スクリプト例は, 以下の通り.

# obtain small size of data
from sklearn.model_selection import train_test_split
train_data_x_sub, x_test, train_data_y_sub, y_test = train_test_split(train_data_x, train_data_y,
                                                                      train_size=3000, test_size=100,
                                                                      random_state=1)

from sklearn.model_selection import GridSearchCV
from sklearn import svm
from sklearn import metrics

param_grid = {"C": [10 ** i for i in range(-5, 6, 2)],
              "gamma" : [10 ** i for i in range(-5, 6, 2)]
             }

model_grid = GridSearchCV(estimator=svm.SVC(kernel="rbf", decision_function_shape="ovo", random_state=1),
                 param_grid = param_grid,   
                 scoring = "accuracy",  # metrics
                 verbose = 2,
                 cv = 4)              # cross-validation

model_grid.fit(train_data_x_sub, train_data_y_sub)
model_grid_best = model_grid.best_estimator_ # best estimator
print("Best Model Parameter: ", model_grid.best_params_)

# Best model parameter : {'C': 10, 'gamma': 0.001}

print('Train score: {}'.format(model_grid_best.score(train_data_x_sub, train_data_y_sub)))
print('Cross Varidation score: {}'.format(model_grid.best_score_))

# Train score: 0.9593333333333334
# Cross Varidation score: 0.9063333333333333

# check calculation results
means = model_grid.cv_results_['mean_test_score']
stds = model_grid.cv_results_['std_test_score']
for mean, std, params in zip(means, stds, model_grid.cv_results_['params']):
    print("%0.4f (+/-%0.04f) for %r" % (mean, std * 2, params))

prediction = model_grid_best.predict(train_data_x_sub)
co_mat = metrics.confusion_matrix(train_data_y_sub, prediction)
print(co_mat)

print('Total Train score: {}'.format(model_grid_best.score(train_data_x, train_data_y)))

# Total Train score: 0.9251904761904762

上記のスクリプトにおいて, 以下の部分は各パラメータの組み合わせにおける計算結果を表示させている. これを見ながら, 次の計算の C と gamma の範囲を絞り込んでいった.

# check calculation results
means = model_grid.cv_results_['mean_test_score']
stds = model_grid.cv_results_['std_test_score']
for mean, std, params in zip(means, stds, model_grid.cv_results_['params']):
    print("%0.4f (+/-%0.04f) for %r" % (mean, std * 2, params))

ちなみに, 私の場合のパラメータの絞り込み方は, 以下の通りで行った. ランダムサーチを使って, もっと時間短縮する方法があるかもしれない.

No	データ数	cv	`C`	`gamma`	最適の `C`	最適の `gamma`	CV score
1	3000	4	`[10 ** i for i in range(-5, 6, 2)]`	`[10 ** i for i in range(-5, 6, 2)`	10	0.001	0.9063
2	3000	5	`[0.3, 1, 3, 10, 30, 100, 300]`	`[0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03]`	3	0.03	0.9443
3	8000	5	`[3, 10, 30, 100, 300]`	`[0.025, 0.03, 0.04, 0.05]`	10	0.025	0.9695
4	8000	5	`[5, 8, 10, 15, 20]`	`[0.015, 0.02, 0.024, 0.028]`	5	0.028	0.9694

今見ると, CV score は一番最後よりも, その 1 つ手前のほうがちょっと良かったのか… :-p

結果

得られた最適の C と gamma を使って, 全データで学習する.

from sklearn import svm
from sklearn import metrics

clf = svm.SVC(C=5, gamma=0.028, decision_function_shape="ovo", kernel="rbf", verbose=2)
clf.fit(train_data_x, train_data_y)
prediction = clf.predict(train_data_x)
accuracy_score = metrics.accuracy_score(train_data_y, prediction)
print(accuracy_score)
# Accuracy : 0.9999761904761905

co_mat = metrics.confusion_matrix(train_data_y, prediction)
print(co_mat)

prediction = clf.predict(test_data)
output = pd.DataFrame({"ImageId" : np.arange(1, 28000+1), "Label":prediction})
output.head()

output.to_csv('digit_recognizer_SVM7a.csv', index=False)
print("Your submission was successfully saved!")

得られた結果の正答率は 0.98375. 上位 50% を, ちょっと下回る.

ちなみに, sklearn.svm.SVC() でデフォルトのハイパーパラメータを使うと 0.97571. 特徴量を増減させずに, ハイパーパラメタだけ最適化すると 0.98207. なので, 特徴量の増減はちょっとは効果がある (苦しいか).

これ以上を狙うなら, ささっと CNN をトライするほうが良いだろう.

参考

Github に上げたスクリプト (digit-recognition_SVM7a.py)
pandas.DataFrameの行を条件で抽出するquery
sklearn.preprocessing.MinMaxScaler の本家マニュアル
MNIST database - wikipedia ; ここによるとサポートベクターマシンでの例として error rate = 0.56 と出ているが, 詳細はよく分からない.