ニューラルネットワークのハイパーパラメーターをベイズ最適化実装　（Chainer+GPyOpt）

#ニューラルネットワークのハイパーパラメーターを自動的に最適化しましょう

ニューラルネットワークのハイパーパラメーター

皆さんはどのように決めていますか？ニューラルネットワークは計算に時間がかかるので，勘や経験に頼ることが多いでしょう。しかし，ある程度の時間内に終わるのであれば（その計算が多くとも30分〜1時間程度で済むのであれば），ベイズ最適化という手があります。ベイズ最適化は最適化したいブラックボックスが，計算に時間がそれなりにかかるものにも威力を発揮します。

おととい，[Pythonでベイズ最適化を行うGPyOptについての記事](http://qiita.com/YoshikawaMasashi/items/51b82a7b990d51bd98cd)を書きました。今回は，この方法を使って，ニューラルネットワークのハイパーパラメーターをMNISTの正答率が高くなるように，最適化しました。実際にPythonでコーディングしてみましたので，みてみてください。

#さっそく結果
まあ，まず最初に結果を示します。

今回行うタスクはMNIST(手書き文字認識）の識別率を向上させることです。

MNISTを訓練データとテストデータに分けて，訓練データでニューラルネットワークを学習させます。そのあとにテストデータの識別を行いますが，この誤答率を下げることを行います。

最適化するハイパーパラメーターは以下の通りです。

[追記 2017/2/8]ステップ数をエポック数に修正したコードをまわしました．

| ハイパーパラメーター | 範囲 |
|:-----------|:------------|
| エポック数       |        1~20 |
| ドロップアウト率       |        0~0.9 |
| 隠れ層のチャンネル数（隠れ層は入力層の何倍のユニットにするか） |          0.51~5 |
| バッチサイズ     |      50,100,200 |
| 畳み込みにするか       |       True,False |
| 層の数    |     1,2,3 |
| residual にするか       |       True,False |
| 畳み込みのフィルターサイズ    |     5,7,9,11,13 |

まず比較実験としてオーソドックスなニューラルネットワークで行なったところ，誤答率は2.58%でした。

そしてベイズ最適化(20回のランダムサンプリング+20回のベイズ最適化のイテレーション)したところ誤答率は1.18%になりました。
最適化された後のハイパーパラメーターは以下のとおりです。

| ハイパーパラメーター | 範囲 | 比較実験 |最適化結果|
|:-----------|:------------|:------------|:------------|
| エポック数       |        1~20 |10|20(上限)|
| ドロップアウト率       |        0~0.9 |0|0|
| 隠れ層のチャンネル数（隠れ層は入力層の何倍のユニットにするか） | 0.51~5 | 1|5(上限)|
| バッチサイズ     |      50,100,200 |100|200(上限)|
| 畳み込みにするか       |       True,False |False|True|
| 層の数    |     1,2,3 |2|3(上限)|
| residual にするか       |       True,False |False|True|
| 畳み込みのフィルターサイズ    |     5,7,9,11,13 |-|9|


#結果を踏まえて
エポック数や隠れ層のチャンネル数は上限になりました．

まあ，今回はパソコンのスペックの問題があって，エポック数を少なめにしましたね（それでも20イテレーションしか回らなかった，特に後の方はエポック数とバッチサイズが上限のものばかり試していたのでめっちゃ時間がかかった）。

エポック数などの上限をあげれば，もっと精度は高くなると思われます。

まあ，それでも最適化されているようでよかったです。

#ニューラルネットワークの方の実装

chainerで実装しました。

test_mnist_nn関数というものを作りました。この関数はニューラルネットワークのハイパーパラメーターを入力にとって，誤答率を出力する関数です。この誤答率を最小にするように，ベイズ最適化を行います。

```py3:mnist_nn.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
パラメーターを最適化する，MNIST用のニューラルネットワークです．
最適化するパラメーターは
epoch_num: エポック数(1~20)
dropout_r:ドロップアウト率(0~0.9)
h_cnl: 隠れ層のチャンネル数（隠れ層は入力層の何倍のユニットにするか）(0.51~5)
batch_size:バッチサイズ(50,100,200)
conv: 畳み込みにするか(True or False)
layers: 層の数 (1~3)
residual: residual にするか(True or False)
conv_ksize: 畳み込みのフィルターサイズ(5,7,9,11,13)
"""

import numpy as np
import chainer
from chainer import optimizers,datasets
import chainer.functions as F
import chainer.links as L

class mnist_nn(chainer.Chain):
    '''
    mnist分類タスクを行うニューラルネットワーク
    '''
    def __init__(self, dropout_r=0.5, h_cnl=0.51, conv=True, layers=3, residual=True, conv_ksize=9):
        '''
        パラメーターの初期化と
        層の作成（Linkの追加）
        を行います
        '''
        super(mnist_nn, self).__init__()
        self.dropout_r = dropout_r
        self.h_cnl = h_cnl
        self.conv = conv
        self.layers = layers
        self.residual = residual
        self.conv_ksize = conv_ksize

        if conv:
            p_size = int((conv_ksize-1)/2)
            self.add_link('layer{}'.format(0), L.Convolution2D(1,round(h_cnl),(conv_ksize,conv_ksize),pad=(p_size,p_size)))
            for i in range(1,layers):
                self.add_link('layer{}'.format(i), L.Convolution2D(round(h_cnl),round(h_cnl),(conv_ksize,conv_ksize),pad=(p_size,p_size)))
            self.add_link('fin_layer', L.Convolution2D(round(h_cnl),10,(28,28)))
        else:
            self.add_link('layer{}'.format(0), L.Linear(784,round(784*h_cnl)))
            for i in range(1,layers):
                self.add_link('layer{}'.format(i), L.Linear(round(784*h_cnl),round(784*h_cnl)))
            self.add_link('fin_layer', L.Linear(round(784*h_cnl),10))

    def __call__(self, x, train=False):
        '''
        ニューラルネットワーク本体です
        '''
        if self.conv:
            batch_size = x.shape[0]
            x = x.reshape(batch_size, 1, 28, 28)
            h = chainer.Variable(x)
            h = F.dropout(F.relu(self['layer{}'.format(0)](h)),train=train,ratio=self.dropout_r)
            for i in range(1,self.layers):
                if self.residual:
                    h = F.dropout(F.relu(self['layer{}'.format(i)](h)),train=train,ratio=self.dropout_r) + h
                else:
                    h = F.dropout(F.relu(self['layer{}'.format(i)](h)),train=train,ratio=self.dropout_r)
            h = self['fin_layer'](h)[:,:,0,0]
        else:
            h = chainer.Variable(x)
            h = F.dropout(F.relu(self['layer{}'.format(0)](h)),train=train,ratio=self.dropout_r)
            for i in range(1,self.layers):
                if self.residual:
                    h = F.dropout(F.relu(self['layer{}'.format(i)](h)),train=train,ratio=self.dropout_r) + h
                else:
                    h = F.dropout(F.relu(self['layer{}'.format(i)](h)),train=train,ratio=self.dropout_r)
            h = self['fin_layer'](h)
        return h
    def loss(self,x,t):
        '''
        誤差関数は交差エントロピーです
        '''
        x = self.__call__(x, True)
        t = chainer.Variable(t)
        loss = F.softmax_cross_entropy(x,t)
        return loss


def test_mnist_nn(epoch_num=10, dropout_r=0.5, h_cnl=0.51, batch_size=100, conv=True, layers=3, residual=True, conv_ksize=9):
    '''
    パラメーター入れたら，MNISTの誤答率を出す関数
    '''
    #一応シードは固定する（これによって，パラメーターに対して出力はexactに定まる
    np.random.seed(1234)

    #モデル最適化の準備
    model = mnist_nn(dropout_r, h_cnl, conv, layers, residual, conv_ksize)
    optimizer = optimizers.Adam()
    optimizer.setup(model)

    #データの用意
    train, test = datasets.get_mnist()
    trainn, testn = len(train), len(test)

    #ログ出力用のlinspace
    logging_num = np.linspace(0,epoch_num*int(trainn/batch_size),11).astype(int)

    iter = 0
    for epoch in range(epoch_num):
        #バッチ作成
        batch_idx = np.array(range(trainn))
        np.random.shuffle(batch_idx)
        for batch_num in np.hsplit(batch_idx,int(trainn/batch_size)):
            batch = train[batch_num]
            x,t = batch[0],batch[1]

            #学習
            model.zerograds()
            loss = model.loss(x,t)
            #念のため
            if np.isnan(loss.data):
                print("pass")
                continue
            loss.backward()
            optimizer.update()

            #ログ出力
            if iter in logging_num:
                print(str(np.where(logging_num==iter)[0][0]*10)+'%','\tcross entropy =',loss.data)
            iter+=1


    #性能評価(メモリ節約のため，ミニバッチで行う)
    batch_idx = np.array(range(testn))
    false_p = []
    for batch_num in np.hsplit(batch_idx,100):
        batch = test[batch_num]
        x,t = batch[0].reshape(len(batch_num), 1, 28, 28),batch[1]
        res = model(x).data.argmax(axis=1)
        false_p.append(np.mean(t!=res))
    false_p = np.mean(false_p)

    print("False {:.2f}%".format(false_p*100))
    return false_p

if __name__ == '__main__':
    test_mnist_nn(epoch_num=10, dropout_r=0, h_cnl=1, batch_size=100, conv=False, layers=2, residual=False, conv_ksize=9)

```

#ベイズ最適化の方の実装

GPyOptを使いました。

最初に20個サンプリングして，その後100イテレーションを行なって最適化しようとしました（実際100イテレーションも回らなかった。）

```py3:bayesian_opt.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jan 27 20:09:10 2017

@author: marshi
"""

import GPyOpt
import numpy as np
import mnist_nn
import os

def f(x):
    '''
    mnist_nn.test_mnist_nnのラッパー関数
    _x[0](1~20) -> epoch_num = int(_x[0])
    _x[1](0~0.9) -> dropout_r = np.float32(_x[1])
    _x[2](0.51~5) -> h_cnl = np.float32(_x[2])
    _x[3](50,100,200) -> batch_size = int(_x[3])
    _x[4](0,1) -> conv = bool(_x[4])
    _x[5](1,2,3) -> layers = int(_x[5])
    _x[6](0,1) -> residual = bool(_x[6])
    _x[7](5,7,9,11,13) -> conv_ksize = int(_x[7])
    '''
    ret = []
    for _x in x:
        print(_x)
        _ret = mnist_nn.test_mnist_nn(epoch_num = int(_x[0]),
                                        dropout_r = np.float32(_x[1]),
                                        h_cnl = np.float32(_x[2]),
                                        batch_size = int(_x[3]),
                                        conv = bool(_x[4]),
                                        layers = int(_x[5]),
                                        residual = bool(_x[6]),
                                        conv_ksize = int(_x[7]))
        ret.append(_ret)
    ret = np.array(ret)
    return ret

#それぞれの変数の領域を指定
bounds = [{'name': 'epochs', 'type': 'continuous', 'domain': (1,20)},
          {'name': 'dropout_r', 'type': 'continuous', 'domain': (0.0,0.9)},
          {'name': 'h_cnl', 'type': 'continuous', 'domain': (0.51,5)},
          {'name': 'batch_size', 'type': 'discrete', 'domain': (50,100,200)},
          {'name': 'conv', 'type': 'discrete', 'domain': (0,1)},
          {'name': 'layers', 'type': 'discrete', 'domain': (1,2,3)},
          {'name': 'residual', 'type': 'discrete', 'domain': (0,1)},
          {'name': 'conv_ksize', 'type': 'discrete', 'domain': (5,7,9,11,13)}]

#ベイズ最適化オブジェクト作成
#前にセーブしたX,Yがあるなら途中から行う
#ないなら最初に数サンプルをランダムにサンプリングする
filename = "XY.npz"
if os.path.exists(filename):
    XY = np.load(filename)
    X,Y = XY['x'],XY['y']
    myBopt = GPyOpt.methods.BayesianOptimization(f=f, domain=bounds, X=X, Y=Y)
else:
    myBopt = GPyOpt.methods.BayesianOptimization(f=f, domain=bounds)

#ベイズ最適化100ステップとステップごとの結果の表示
for i in range(1000):
    print(len(myBopt.X),"Samples")
    myBopt.run_optimization(max_iter=1)
    print(myBopt.fx_opt)
    print(myBopt.x_opt)
    #逐次セーブ
    np.savez(filename, x=myBopt.X, y=myBopt.Y)

```