ベイズ最適化シリーズ（3） －XGBoostのハイパーパラメータ探索－

今回から「ベイズ最適化」関連は、ベイズ最適化シリーズと題して
お送りしたいと思います。今回は第３回目です。

ベイズ最適化シリーズ（1）　－ベイズ最適化の可視化－
ベイズ最適化シリーズ（2）　－アンサンブル学習（Voting）の最適化－

今回は、やる♪やる♪といいながら、一か月ほど保留になっていた
XGBoostのハイパーパラメータ探索を行います。

実は、ずっと放置していたわけではなく、しばらくハマっていました。
詳細は、後ほど述べます。

#XGBoostとは
私が説明するよりも良い記事がたくさんありますので、そちらをご覧ください。
例えば、以下の記事が挙げられます。
https://qiita.com/yh0sh/items/1df89b12a8dcd15bd5aa

#ハイパーパラメータ
XGBoostには10種類以上のハイパーパラメータがあります。
ハイパーパラメータの設定によって、精度が数パーセント変わってくる
と言われています。

このハイパーパラメータの探索は、総当たりによる探索や、徐々に範囲を絞っていく手法も
ありますが、時間がかかり過ぎたり、面倒くさい作業が増えたりするデメリットもあります。

#ベイズ最適化によるハイパーパラメータ探索
ベイズ最適化を使えば、短時間で良いハイパーパラメータを探してくれます。

ただし、ベイズ最適化はガウス過程に基づいて探索するため、その過程から
外れる場合は、うまくいかない可能性があります。例えば、多峰性のある
ノイズが入る場合は苦手です。

いずれにしても、デフォルトのモデルで学習させるよりは、ベイズ最適化で
見つけたモデルで学習させた方が、はるかに良い精度が出ると思います。

さっそく実装します。使うデータはお馴染みのMNISTです。

import numpy as np
import pandas as pd

from keras.datasets import mnist

#MNIST データロード
(X_train, y_train), (X_test, y_test) = mnist.load_data()
X_train = X_train.reshape(60000, 784)
X_test = X_test.reshape(10000, 784)

x_train = X_train.astype('float32')
x_test = X_test.astype('float32')
x_train /= 255
x_test /= 255

x_train = x_train[:1000]
x_test = x_test[:500]
y_train = y_train[:1000]
y_test = y_test[:500]

XGBoostを定義し、デフォルトのハイパーパラメータで精度を出してみます。

import xgboost as xgb

from sklearn.metrics import accuracy_score
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score

#XGBoost
Xgb = xgb.XGBClassifier()
Xgb.fit(x_train,y_train)

xgb_predict = Xgb.predict(x_test)
print("Default XGBoost")
print(accuracy_score(y_test, xgb_predict))

Default XGBoost
0.842

精度は84.2％と出ました。

次に、GPyOptのベイズ最適化を使って、ハイパーパラメータを探索します。

#ベイズ最適化
import GPy
import GPyOpt
import matplotlib.pyplot as plt

def f(x):
    model = xgb.XGBClassifier(num_class = 10, 
                              eta = float(x[:,0]),
                              max_depth = int(x[:,1]),
                              min_child_weight = float(x[:,2]),
                              subsample = float(x[:,3]),
                              colsample_bytree = float(x[:,4]),
                              gamma = float(x[:,5]),
                              n_estimators = int(x[:,6]),
                              learning_rate = float(x[:,7]),
                              reg_lambda = float(x[:,8]),
                              reg_alpha = float(x[:,9]),
                              num_boost_round = int(x[:,10]))

    # CV
    kfold = KFold(n_splits=10, random_state=7)
    results = cross_val_score(model, x_train, y_train, cv=kfold)
    score = results.mean()*100
    print("Accuracy: %.2f%%" % (results.mean()*100))

    return -score

bounds = [{'name': 'eta', 'type': 'continuous', 'domain': (0,5)},
          {'name': 'max_depth', 'type': 'continuous', 'domain': (150,300)},
          {'name': 'min_child_weight', 'type': 'continuous', 'domain': (1,5)},
          {'name': 'subsample', 'type': 'continuous', 'domain': (0,1)},
          {'name': 'colsample_bytree', 'type': 'continuous', 'domain': (0.1,1)},
          {'name': 'gamma', 'type': 'continuous', 'domain': (0,2)},
          {'name': 'n_estimators', 'type': 'continuous', 'domain': (10,50)},
          {'name': 'learning_rate', 'type': 'continuous', 'domain': (0,1)},
          {'name': 'reg_lambda', 'type': 'continuous', 'domain': (0,1.1)},
          {'name': 'reg_alpha', 'type': 'continuous', 'domain': (0,1.1)},
          {'name': 'num_boost_round', 'type': 'continuous', 'domain': (200,400)}]

print("Bayesian Optimization")
myBopt = GPyOpt.methods.BayesianOptimization(f=f, domain=bounds)
myBopt.run_optimization(max_iter=50)

#探索履歴描画
result_z = myBopt.Y

plt.figure()
plt.plot(-result_z)
plt.xlabel("epochs")
plt.ylabel("CV Accuracy")
plt.show()

#最適なパラメータの表示
print("best parameters =")
print(myBopt.x_opt)
x = myBopt.x_opt

Xgb = xgb.XGBClassifier(num_class=10, eta=x[0], max_depth=int(x[1]),min_child_weight=x[2],
                        subsample=x[3], colsample_bytree=x[4], gamma=x[5], n_estimators=int(x[6]),
                        learning_rate=x[7], reg_lambda=x[8], reg_alpha=x[9],num_boost_round=int(x[10]))
Xgb.fit(x_train,y_train)

#最適なパラメータで再度学習、結果表示
xgb_predict = Xgb.predict(x_test)
print("Optimized XGBoost")
print(accuracy_score(y_test, xgb_predict))

Optimized XGBoost
0.848

精度は84.8％となり、デフォルト値に対し0.6％改善しました。
改善幅は思った以上に小さかったです。

#学習のコツ
MNISTを使ってハイパーパラメータ探索をしましたが、いくつか注意すべきポイントが
あったので記しておきます。

・探索するパラメータの範囲を大きく、そしてイタレーションの数を大きくとれば、
　一発でベストなパラメータを見つけてくれると思ったのですが、現実は甘くなかったです。
　結論から言うと、（パラメータの範囲を大きく、イタレーションは200くらい）→
　（範囲を狭めていく、イタレーションは50くらい）と段階的に絞っていかないと、
　デフォルトのパラメータに負けてしまいます。（かれこれ50回以上は負けました。）

　理想像からすれば、手間は増え、かつ探索時間も増えますが、既存の手法に比べれば
　手間は減ります。感覚的には、（既存：20回作業）→（ベイズ最適化：5回作業）になる
　レベルです。

・ベイズ最適化で観察する関数は、クロスバリテーションの分類精度を使いました。
　具体的には、下記のscikit-learnのcross_val_scoreを使いました。学習データを
　分けても良いのですが、クロスバリテーションのスコアを使えば間違いないです。

from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score

kfold = KFold(n_splits=10, random_state=7)
results = cross_val_score(model, x_train, y_train, cv=kfold)
score = results.mean()*100