Pythonでベイズ最適化を使ってハイパーパラメータを探索するライブラリ実装のメモ

はじめに

ベイズ最適化(参考：ベイズ最適化入門, 機械学習のためのベイズ最適化入門)を使うと、機械学習の時の各種Try&Errorで決めざるを得ないようなハイパーパラメータの探索を効率よく実施できる可能性があります。

考え方などは最近色々解説記事が増えてきたおかげで理解はできるのですが、GridSearchのライブラリみたいな形ではWeb上で見つけられなかったので、今回作りました。きっと車輪の再発明なのだと思うのですが、まあ再発明は勉強にはなるので良しとします。

今回使っている各種Version

• Python3.5
• numpy==1.11.1
• scikit-learn==0.17.1

コード

bayesian_optimizer.py

from itertools import product
from sklearn.gaussian_process import GaussianProcess

# The MIT License (C) 2016 mokemokechicken


class BayesianOptimizer:
    x_list = None
    y_list = None
    yielding_index = None
    k_band = 5
    verbose = False

    def __init__(self, params):
        self.params = params
        self.keys = []
        self.values = []
        for k, v in sorted(self.params.items()):
            self.keys.append(k)
            self.values.append(v)

    @property
    def n_pattern(self):
        return len(list(product(*self.values)))

    def output(self, *args, **kwargs):
        if self.verbose:
            print(*args, **kwargs)

    def supply_next_param(self, max_iter=None):
        self.x_list = []
        self.y_list = []
        all_parameters = list(product(*self.values))  # [(0.01, [0, 0], 0, [10, 10]), (0.01, [0, 0], 0, [15, 15]), ...
        index_space = [list(range(len(v))) for v in self.values]  # [[0], [0, 1, 2], [0], [0, 1, 2]]
        all_index_list = list(product(*index_space))  # [(0, 0, 0, 0), (0, 0, 0, 1), ...

        # examine 2 random points initially
        idx = list(range(len(all_index_list)))
        np.random.shuffle(idx)
        searched_index_list = []
        for index in idx[:2]:
            param = self.to_param(all_parameters[index])
            self.yielding_index = all_index_list[index]
            searched_index_list.append(index)
            yield param

        # Bayesian Optimization
        max_iter = int(min(max_iter or max(np.sqrt(self.n_pattern)*4, 20), self.n_pattern))  # 最大探索回数を適当に計算。
        for iteration in range(max_iter):
            k = 1 + np.exp(-iteration / max_iter * 3) * self.k_band  # kの値を徐々に減らして 探索重視 → 活用重視 にしてみる
            gp = self.create_gp_and_fit(np.array(self.x_list), np.array(self.y_list))

            mean_array, mse_array = gp.predict(all_index_list, eval_MSE=True)
            next_index, acq_array = self.acquisition(mean_array, mse_array, k, excludes=searched_index_list)

            self.output("--- Most Expected Predictions")
            for acq, ps in sorted(zip(acq_array, all_parameters), reverse=True)[:3]:
                self.output("%.2f: %s" % (acq, list(zip(self.keys, ps))))
            self.output("--- Past Best Results")
            for acq, vs, ps in self.best_results(3):
                self.output("%.2f: %s" % (acq, list(zip(self.keys, vs))))

            if next_index in searched_index_list:
                break
            searched_index_list.append(next_index)
            self.yielding_index = all_index_list[next_index]
            yield self.to_param(all_parameters[next_index])

    @staticmethod
    def create_gp_and_fit(x, y, max_try=100):
        # この辺怪しい
        theta0 = 0.1
        for i in range(max_try+1):
            try:
                gp = GaussianProcess(theta0=theta0)
                gp.fit(x, y)
                return gp
            except Exception as e:
                theta0 *= 10
                if i == max_try:
                    print(theta0)
                    raise e

    def to_param(self, row):
        return dict(zip(self.keys, row))

    def report(self, score):
        self.x_list.append(self.yielding_index)
        self.y_list.append(score)

    def best_results(self, n=5):
        index_list = []
        param_list = []
        for xs in self.x_list:
            values = [self.values[i][x] for i, x in enumerate(xs)]
            index_list.append(values)
            param_list.append(self.to_param(values))
        return sorted(zip(self.y_list, index_list, param_list), reverse=True)[:n]

    @staticmethod
    def acquisition(mean_array, mse_array, k, excludes=None):
        excludes = excludes or []
        values = mean_array + np.sqrt(mse_array) * k
        for_argmax = np.copy(values)
        for ex in excludes:
            for_argmax[ex] = -np.Inf
        return np.argmax(for_argmax), values

使い方

基本的な使い方は以下のようになります。

params = {
    "parameter1": [10, 20, 25, 50],
    "parameter2": ['p1', 'p2'],
    ...
}

bo = BayesianOptimizer(params)  # パラメータの組み合わせを辞書で渡す

for param in bo.supply_next_param():  # 次調べたほうが良いparamをdictで渡される
    y = unknown_score_function(param)              # なんらかの方法でその param を評価する
    bo.report(y)                                    # 評価値を教えてあげる

print(bo.best_results())                            # 一応一番良いScoreの値やparamは覚えてくれている

２次元空間でのデモ

ちゃんと動くか Jupyter notebook上で試したものを以下に置いてます。
https://gist.github.com/mokemokechicken/7f3cf33c71d33bf0ec242c37cb3f2a75

探索空間

以下の様な 50 x 50 平面に sin, cos で起伏を複数作った後、ノイズを乗っけてちょっと凸凹した２次元空間を用意します。

赤い点がある (row=3, col=29) の部分が最も値の大きい地点です。

探索

こんな感じで進んでいきます。
探索回数は200回です。

序盤

まずは全体的に空間を調べています。少し囲碁っぽいｗ。

序盤２

運良く、真ん中上にある本命の近くを何度か調べたようです。

中盤

本命の地点は結構調査済みになりました。

中盤２

もう一つの良さそうな左下の山も丹念にチェック。真ん中上の本命周辺も結構網羅。

終盤

他に良さそうなところがないか空白地点をチェック。

最終状態

他に特に良さそうなところも無く終了。まあ、人間がやってもこんな感じになりそうです。

所感

• 200回の探索で50x50のパラメータ空間から良いポイントを見つけられました。必ずしも最善の地点を見つけられはしませんが、センスの良い探索をしてくれています。
• 今回のように探索空間が凸凹していると(?) sklearnのGaussianProcessの fit() で Exception が発生することがしばしばあります。 「theta0 や thetaU を 大きくしたらどうだ」と言われるので、今のところそれで対応しています（それで良いのかよくわかってないです）。実際に機械学習のハイパーパラメータの評価のケースでは、結構ノイズも大きいのでこの問題はしばしば起こります。
• 今回の実装では、パラメータの「値」ではなく「配列上の位置」を使って共分散を計算しているので、文字列パラメータにも対応していますが、順序が無い値(例えば、 ['dog', 'cat', 'bird'] とか)だと少し変な感じです。

さいごに

作ってみると、結構色々なシーンで使えそうな感じです。
人生もこんな感じだよなぁ。