GPyOpt: Digital Recognizer(MNIST) CNN Keras ハイパーパラメータ最適化

引き続きハイパーパラメータ最適化として、今回はGPyOptを使ってみました。これまでHyperas、Hyperoptを試して見ましたが、ベイズ最適化でも採用しているアルゴリズムが微妙に違うようなので試してみたという感じです。
個人的にはHyperoptが一番使いやすく感じましたが、GPyOptは以前scikit-learnで試したベイス最適化に近いアルゴリズムだったのでもう少し理解を深めてみようかと。
まだ手探り段階なので、使い方に関しては後から追記するかもしれません。


使い方：

ハイパーパラメータの設定：
まずはMNISTモデルにおけるハイパーパラメータの設定からです。今回はやや少なめで。

最適化するハイパーパラメータ：
・各層のドロップアウト率：連続値
・Dense層出力ユニット数：離散値
・validation_splitの比率：連続値

GPyOptではハイパーパラメータを以下のようなフォーマットで書きます。

params = [
    {'name': 'Dropout_0',        'type': 'continuous',  'domain': (0.0, 0.5)},
    {'name': 'Dropout_1',        'type': 'continuous',  'domain': (0.0, 0.5)},
    {'name': 'Dropout_2',        'type': 'continuous',  'domain': (0.0, 0.5)},
    {'name': 'Dropout_3',        'type': 'continuous',  'domain': (0.0, 0.5)},
    {'name': 'Dense_0',          'type': 'discrete',    'domain': (128, 256, 512)},
    {'name': 'Dense_1',          'type': 'discrete',    'domain': (64,128, 256)},
    {'name': 'validation_split', 'type': 'continuous',  'domain': (0.1, 0.3)}
]

リスト化されたディクショナリーで（ここを参考に）、

・'name' : パラメータ名
・'type' : 'continuout'（連続値）、'discrete'（離散値）、'categorical'（分類値）
・'domain' : 適用範囲または選択肢を()で括る

となるようです。


CNNモデルの構築：
次にモデルを構築します。前回同様MNIST分類用のCNNを使います。このモデルからはベイズ最適化するための評価値となるlossかaccが求められればいいのですが、

・loss
・acc
・model
・history

の４種類を戻り値にしておきました。
model.fit()させてEarlyStoppingで打ち切りになった最後のval_lossとval_accの値を参照しています。

    loss = hist.history['val_loss'][-1]
    acc = hist.history['val_acc'][-1]

modelやhistoryは不要ですが、後から参照するかもしれないので一応入れておきました（使うかどうかは分からない）。

return loss, acc, model, hist

上記パイパーパラメータに対応する変数部分には、x[:, 0]などと引数にインデックス番号をつけるようですが、どれが何番目かはわかりにくいのでハイパーパラメータの'name'から参照できる関数をつくってみました。

model.add(Dropout(Param('Dropout_0'), seed=seed))

このように書き込めばx[:, 0]へ自動変換してくれます。Hyperoptなどでもディクショナリーのキーを使っていたので、このほうが個人的には使いやすいかと（リスト内容を変えた場合にインデックス番号だと、他の番号も変わってしまうのが面倒なので）。
注意点として、最初に書いたハイパーパラメータはリストであるのに対して、この変数は２次元のndarrayに変換されてから代入されるようです。この変換関数は以下（cnn_model関数内）。

    def Param(p_name):
        p_index = [p['name'] for p in params].index(p_name)
        p_type = params[p_index]['type']
        
        if type(x) is np.ndarray:
            if p_type == 'continuous':
                return float(x[:, p_index])
            else:
                return int(x[:, p_index])
        else: # list
            if p_type == 'continuous':
                return float(params[p_index]['domain'])
            else:
                return int(params[p_index]['domain'])

後で最適化されたハイパーパラメータリストを直接渡せるようにしてあります。引数がndarrayならx[:,0]のような2次元ndarray、listならlist内のスカラー値へ変換後代入。また今回の場合、離散値はすべて整数だったのでintかfloatかも振り分けています。


フィッティング関数：
上記CNNモデルを以後のベイズ最適化関数GPyOpt.methods.BayesianOptimization()に直接渡してもいいのですが、CNNモデルからは４種類の値を出力することにしたので、このf(x)関数を間にはさんで必要な評価値だけを渡せるようにしました。今回はaccを評価値として渡すことにし、最小化するためにマイナス反転して-accにしています。
前述のように引数のxは二次元のnumpy.ndarrayになるようです。今回は7種類のハイパーパラメータがあるので、x.shapeは(1,7)になります。最初に設定したハイパーパラメータはリストでありndarrayではないので、この辺をいじる場合は変換するなどの工夫が必要です（このサンプルを参照）。
実際は、

def f(x):
    x = np.atleast_2d(x)
    fs = np.zeros((x.shape[0],1))
    for i in range(x.shape[0]):
        loss, acc, model, hist = cnn_model(x)
        fs[i] += np.log(acc)*(-1)
    return fs

このように書いたほうがいいのかもしれませんが、戻り値は1次元のndarrayだったので、今回は省略して以下のようにしました。対数変換したほうがいいのかもしれませんが効果の違いは検証していません。

def f(x):
    loss, acc, model, hist = cnn_model(x)
    return -acc

ベイズ最適化関数：
GPyOpt.methods.BayesianOptimization()に先程のf(x)関数とハイパーパラメータリストparamsを渡し、その他初期探索値や獲得関数などを決めます。獲得関数はデフォルトではEIになっていますが'EI_MCMC'を選んでみました。'EI_MCMC'を選択する場合は、model_typeで'GP_MCMC'を選んでおかなければいけないようです。
initial_design_numdataは20に設定しましたが、これはどのくらいがいいのかは不明（デフォルト：5）。探索する前のランダムな開始点の数なのかもしれませんが、今回の7次元に対してどのくらいが適当なのか？探索点は徐々に追加されながらフィッティングしていくと思うのでデフォルトの5でもいいのかもしれません。入れた回数だけループするようです（20回で約1時間）。
こまかな設定がいくつかありますが、まだ使いながら試している段階です。

次に、run_optimization(max_iter=50)で最適化が始まります。イテレーションを50回に設定しました。7種類のハイパーパラメータに対してどのくらいが適当なのかはまだ不明（ハイパーハイパーパラメータ）。50回で約4時間かかりました。
ループが終了すれば最適なハイパーパラメータが見つかったことになります。設定した回数より早く終わることもあります。


最適化されたハイパーパラメータの取得：
以下で結果を取得することができます。

x_best = opt.x_opt
print([i for i in x_best])

y_best = opt.fx_opt
print(y_best)

そうすると、

[0.1732254530746627, 0.39555160207057505, 0.14877909656106353, 0.07323704794308367, 128.0, 128.0, 0.1471216716379693]
-0.9945388349514563

と値が出てきて、最初のリストが最適化された各ハイパーパラメータ。
下の値はそのときのロス値。今回はaccをマイナス反転してあるのでaccの値と同じ。精度0.994以上でているのでまあまあの結果です。


最適化されたハイパーパラメータをモデルに適用：
上記結果と同時にベストモデルやベストウェイトを直接取り出したいのですが、そのような方法がGPyOptにはないようなので、最適化されたハイパーパラメータをCNNモデルに入れ直して再度訓練させてみました。
一応、上記ハイパーパラメータリストを元々のディクショナリー型のリストへ移し替えてからCNNモデルに渡しています。CNNモデルの引数がlistの場合はスカラー値を各変数に代入するような関数にしています。
CNNモデルはEarlyStopping機能をつけているので15ループで収束してくれました（4分25秒）。
このモデルを利用して提出用データを予測します。


まとめと結果（スコア）：
最終的にスコアは0.99457でした。まあまあいい結果です（それでも手動調整のベストスコアである0.99528には達していない）。約6時間でこの結果ですが、もっと回せば向上するかはわからないです。これ以上のスコアを出すには、data augmentationでデータを水増しするなど必要かもしれません。
GpyOptはHyperoptに比べるとやや使いにくいという印象でした（サンプルも少ない）。しかしやりたいことに応じて使いやすく改造すればいいのかもしれません。もともとのアルゴリズム自体は優れていると思うので、いくつかを同時に試して結果的にいい方を選ぶ感じでしょうか。時間的にもHyperoptのほうが速いかもしれませんが、どのライブラリであっても数時間はかかるので時間よりも精度がでるほうがいいと思います（仕事で使っているわけではないので）。
このほか気になるライブラリとして、Skopt、Kopt、PyBO、SpearMintなどありますが、とりあえずはもう十分かと。

これまでは機械学習理論やアルゴリズムの種類を覚えていくことが面白かったのですが、Kaggleをきっかけにスコア（精度）を少しでもあげようとすることにも興味を持てたのはよかったです。実際使ってみて、その結果から次にどうすればいいのかという具体的な疑問が次のモチベーションになるので、より理解も深まりつつ面白くなっていく感じです。

追記：
その後、4つのDropout率だけをハイパーパラメータとして最適化した結果スコア:0.99557まで向上（これまでのベストスコアは0.99524）。
その他のハイパーパラメータは以下のように固定。
validation data：test_size=0.15
Dense_0 output units: 256
Dense_1 output units: 128
batch_size=32

そして最適化においては以下の探索回数に設定。
initial_design_numdata=30（2h 29mins）
max_iter=100(stop at 52: 7h 47mins)
max_iterは最大100回に設定しましたが途中52回で収束し停止しました。
合計で10時間30分（GTX1060で）。

関連：
Kaggle Digital Recognizer(MNIST): Hyperopt + Data Augmentation
Kaggle Digital Decognizer(MNIST): Keras, fit_generator() + hyperopt

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

Posted with Amakuri at 2018.12.21

須山 敦志

講談社

販売価格 ￥3,024(2018年12月21日20時40分時点の価格)

Amazonで詳細を見る

Bayesian Optimization / Gaussian Process

前回、matplotlibのアニメーションについては大体使えるようになったので、「ベイズ最適化・ガウス過程」を実験してみました。いろいろと専用のライブラリがあるのですが（GPy、GPFlow、Edward、Pyroなど）、今回はscikit-learn Gaussian Processesを使っています。

ベイズ最適化はハイパーパラメータチューニングに使われるようで、ブラックボックス関数における最大値（あるいは最小値）を全領域をスキャンしなくても探し出すことができる方法のようです。

今回の場合は、グラフ内の赤破線の最大値がどこにあるか探索しています。赤破線を見ると、左と右にピークが二箇所ありますが、微妙に右側のほうが高いという設定にしてあります。それを間違えないように探し出せるかという実験です。

右側の赤丸が正解（最大値）の地点です。この地点は、x座標-2〜2までの範囲を0.1刻みで200回スキャンし、その結果から最大値を求めています。

しかし、ベイズ最適化の場合は、200回スキャンしなくても約10回程度で、どこが最大値なのかを効率よく探してくれるようです。

今回の手順として：
・カーネル：Matern、獲得関数：UCBを選択

・xの範囲（-2〜2）からランダムにまず2点（青丸：start印）選んで初期観測値とする

・その観測値から、全域のμ（青破線）とσ（青塗）を求める

・μとσを獲得関数（UCB）に代入

・獲得関数の結果（緑実線）から最大値の箇所を見つける

・獲得関数の最大値から、それに相当するxの位置を見つける

・xの位置から、それに相当するyの位置を求める

・ここで求めたxとyの座標を次回の探索箇所にする

・次回探索箇所を観測値リストに加える

・あとは繰り返し

グラフ下部（緑実線）にあるのが獲得関数UCBの値です。緑実線の一番高いところが、次の探索地のx座標となります。

今回使用したUCBの中身は単純で、

UCB最大値 = mu + kappa * sigma

となっています。係数kappa=1.0の場合は、μ値（青破線）にσ値（青塗）を足した値がそのままUCB最大値となるので、グラフ内の青塗を含めた一番高い箇所が次の探索箇所となります。μ値が高いところほど選ばれがちですが、kappaを大きくするとσ値も大きくなるので、μ値が低い箇所でも探索される可能性があがります。

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
import matplotlib.animation as animation
from IPython.display import HTML

np.random.seed(111)

def f(x):
    return -np.sin(x*4.5)*x*0.3+0.02*x                  # 真の関数（赤破線表示）

x = np.linspace(-2,2,200)                               # x値範囲
y = f(x)                                                # 真のy値

fig = plt.figure(figsize=(12,6))                        # 画面サイズ設定
ax = plt.subplot()                                      # サブプロットを用意
ax.axis([x.min(),x.max(),y.min()*2,y.max()*2])          # グラフ範囲設定（x,yの値に依存）
ax.set_aspect(1.0)                                      # アスペクト比設定

ax.plot(x,y,'r--',alpha=0.8)                            # 真の関数曲線をプロット（赤破線）
ax.plot(x[np.argmax(y)],y.max(), 'ro',alpha=0.8)        # y正解値；x全域スキャンによるy最大値（赤丸）

init_pt = 2                                             # 初期観測値個数
xi = np.random.uniform(x.min(),x.max(),init_pt)         # 初期観測値：xの範囲から2箇所ランダム選択
yi = f(xi)                                              # 初期観測値：yの値
ax.plot(xi,yi,'bo')                                     # 初期観測値プロット（青丸）
for i in range(len(xi)):
    ax.text(xi[i],yi[i],'  start'+str(i),color='b')     # 初期観測値：'start'文字表示

def UCB(mu_x, sigma_x, kappa=1.0):                      # 今回使用する獲得関数：upper_confidence_bound
    return mu_x + kappa * sigma_x                       # kappa増大-->σ増大-->探索領域拡大

kernel = Matern(length_scale_bounds="fixed")            # 今回使用するカーネル：Matern
gp = GaussianProcessRegressor(kernel=kernel,alpha=0.01) # scikit-learnガウス過程回帰オブジェクト

ims = []                                                # アニメ用リスト
itr = 15                                                # 探索ループ回数

for i in range(itr):
    gp.fit(xi[:,None],f(xi))                                 # 観測値フィッティング（xiをリシェイプして代入）
    mu_x, sigma_x = gp.predict(x[:, None], return_std=True)  # μ(x)とσ(x)を獲得（事後分布）

    a_x = UCB(mu_x, sigma_x, kappa=2.0)                      # UCBでa(x)最大値を得る
    index = np.argmax(a_x)                                   # a(x)最大値のindexを得る
    x_next = x[index]                                        # a(x)最大値のindexから次のx探索値を決定
    y_next = f(x_next)                                       # y探索値決定
    
    P1 = ax.plot(x, mu_x, 'b--')                                               # μ値（青破線）
    P2 = ax.fill_between(x, mu_x-sigma_x, mu_x+sigma_x, color='b', alpha=0.2)  # σの範囲（青塗）
    P3 = ax.plot(xi,yi,'bx',markersize=8)                                      # 探索済み箇所（青クロス）
    P4 = ax.plot(x_next, y_next, 'mo')                                         # 現在探索箇所（マジェンタ丸）
    P5 = ax.plot([x_next,x_next],[y.min()*2, y_next] ,'m-', alpha=0.6)         # 現在探索箇所（マジェンタ縦線）
    P6 = ax.plot(x_next,(a_x[index]-a_x.min())/2+y.min()*2, 'go',alpha=0.8)    # a(x)最大値：σ（青丸）
    P7 = ax.plot(x, (a_x-a_x.min())/2+y.min()*2, 'g-')                         # 獲得関数値：a(x)
    P8 = ax.plot([xi,xi],[[y.min()*2]*len(yi),yi], 'b-', alpha=0.2)            # 探索済み箇所（青縦線）
       
    xi = np.r_[xi,x_next]                                  # 現在探索箇所をリストに加える（r_: concatenateの代用）
    yi = np.r_[yi,y_next]
    
    ims.append(P1 + [P2] + P3 + P4 + P5 + P6 + P7 + P8)    # アニメ用リストへ各プロットを入れる
    
print('x :',x[y.argmax()], '   y :',y.max())               # 正解値出力（y値）
print('xi:',xi[yi.argmax()], '   yi:',yi.max())            # 最終最大値出力（yi値）

ani = animation.ArtistAnimation(fig, ims, interval=1000)   # アニメ関数
#ani.save("gp_sim.gif", writer = "imagemagick")            # GIFアニメ保存
plt.close()
HTML(ani.to_html5_video())                                 # HTML5 Video 出力（mp4ダウンロード可）
#HTML(ani.to_jshtml())                                     # JavascriptHTML 出力

GaussianProcessRegressor()にはカーネルを指定する箇所があり、今回は、scikit-learnのMaternを使っています。ちなみにカーネルを指定せずデフォルトのままだと、以下のような感じになります。

σ（青塗）の領域が結構変わります。
scikit-learnのこのページにカーネルの違いについて書いてあります。

ハイパーパラメータチューニングと言っても、カーネルの選び方でもかなり違いがでてきます。
また、獲得関数もUCB以外にいくつかあるので、それも選ばなければいけないという感じです。


参考にしたサイト：
Scikit-learn Gaussian Processes
llSourcell/hyperparameter_optimization_strategies
ベイズ最適化入門
機械学習のハイパーパラメータ探索 : ベイズ最適化の活用

関連：
GPyOpt: Digital Recognizer(MNIST) CNN Keras ハイパーパラメータ最適化