実用的に使える多次元入力ベイズ最適化をPython『Bayesian Optimization』でやってみた【GPyOptなし】

テクノロジー

Pythonで簡単にベイズ最適化

Pythonでベイズ最適化をする方法には,いくつかあります.ひとつは,自分で一からアルゴリズムを書く方法で,もう一つは,ライブラリを使って勝手に計算してもらう方法です.

Pythonで使えるベイズ最適化のライブラリには,Bayesian Optimization,GPy,GPyOpt,scikit-optimizeなどがあります.

今回は,個人的に最も簡単に使えるBayesian Optimizationを使いました.

GPyOptが良さそうだったんですが,インストールに引っかかり,心が折れたので簡単にインストールできたBayesian Optimizationを使いました.

環境

Windows10

(anacondaは無し)

他の環境やインストール方法なら,GPyOptのインストールもできたかもしれません(インストール方法は後述).

実行結果

結果から先に.こんな感じでベイズ最適化できます.

一変数関数

緑の測定データからベイズ最適化.赤色の点がベイズ最適化から得られた最大点で,本来の関数f(x)=-x^2の最大値からは少しずれている.オレンジ色の線は予測された関数.

二変数関数

黒の点が測定データ,緑の点がベイズ最適化から得られた最大点.真の関数はf(x_1, x_2) = -x_1^2 -3x_2^2である.

今回やったこと

大学院の研究でデータ点から関数を予測して最適化しようと思い,ネットでベイズ最適化のプログラムを拾おうとしたんですが,ネットで調べた記事は,ほとんどが「与えた関数をベイズ最適化」するものでした.

実際には与えられた関数は分からないわけですから,実用的に使えるように,『データ点を与えてベイズ最適化したい!』と思い,この記事を書きました.

ここで言う『データ点』とは,変数とそれに対する正解(ターゲット点)のことです.

データ点からベイズ最適化

実験を行うとき,f(x)の形は分からないけど,何回かの実験を行ってデータ点を取得し,f(x)が最大(または最小)になるような点x_maxを見つけます.

例えば実験でx=5をやってみたらf(5)=15でした,x=10をやってみたらf(10)=8でした,みたいなのを何回かやって,そのデータ点のみから最大のx_maxおよびf(x_max)を見つける感じです.

ちなみに,「次にどのxを実験すればよいのか」はベイズ最適化のプログラムが決めます(獲得関数を設定する).

・パラメータを多次元入力でもやりました.

・最適化の様子をアニメーションにしました.

Bayesian Optimizationをインストール

pipが使える環境なら,これでインストールできます.

 pip install bayesian-optimization

また,今回使用するライブラリは以下の通りなので,プログラムを実行する前にインストールしておいてください.

・numpy

・matplotlib

ベイズ最適化プログラム

主に参考にした記事は,こちらです.

ベイズ最適化シリーズ(5) -リアルな実験で使ってみたい- - Qiita
ベイズ最適化シリーズ5回目。 今回は、リアルな実験でベイズ最適化を使う方法をご紹介します。 ベイズ最適化シリーズ(1) -ベイズ最適化の可視化- ベイズ最適化シリーズ(2) -アンサンブル学習(Voting)の最適化- ベイズ最適化...

こちらの記事もデータ点からベイズ最適化を行っています.こちらの記事では測定データ(結果)を手入力してますが,今回は正解の関数に乱数を加えた関数で測定データを疑似的に与えています.

つまり,正解の関数をピッタリ最適化できるとは限らず,ベイズ最適化の性能次第ということです.

この部分を,自分が与えたい正解データを生成するプログラムを書くか,正解データを手入力するように書きかえれば自由にベイズ最適化を行うことができます.

プログラムは以下です.コピペしてそのまま実行してください.

一変数のベイズ最適化

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import math
import numpy as np
from bayes_opt import BayesianOptimization
from matplotlib import pyplot as plt
import random

#予備実験のデータ
initial_x = np.array([-1, 0.5])
initial_y = np.array([1, 0.25])

initial_x = initial_x.reshape((-1,1))
initial_y = initial_y.reshape((-1,1))

#実験回数
exp_num = 5

try_ = 1

def f(x):
    global try_

    print("Try : " , try_ , " / {0} , next x is ".format(exp_num) , x)
    #score = float(input("Input y : "))
    score = -x**2 + 10*random.uniform(-1,1)
    print('Input y :', score)

    try_ += 1

    return score

def plot_bo(bo):
    # プロット範囲 (決め打ち)
    X = [x for x in np.arange(-5, 15, 0.1)]

    # 真の関数
    y = [-x**2 for x in X]
    plt.plot(X, y, label='true')

    # サンプル点
    #print(bo.res[0])
    xs = [bo.res[p]['params']['x'] for p in range(len(bo.res))]
    ys = [bo.res[p]['target'] for p in range(len(bo.res))]
    plt.scatter(xs, ys, c='green', s=20, zorder=10, label='sample')

    # 予測結果
    mean, sigma = bo._gp.predict(np.array(X).reshape(-1, 1), return_std=True)
    plt.plot(X, mean, label='pred')  # 推定した関数
    plt.fill_between(X, mean + sigma, mean - sigma, alpha=0.1)  # 標準偏差

    # 最大値
    max_x = bo.max['params']['x']
    max_y = bo.max['target']
    print('max = ', max_x, max_y)
    plt.scatter([max_x], [max_y], c='red', s=50, zorder=10, label='pred_max')


def main():
    # 探索するパラメータと範囲を決める
    pbounds = {
        'x': (-5, 15),
    }

    # 探索対象の関数と、探索するパラメータと範囲を渡す
    bo = BayesianOptimization(f=f, pbounds=pbounds)
    # 最大化する
    bo.maximize(init_points=3, n_iter=10)

    # 結果をグラフに描画する
    plot_bo(bo)

    # グラフを表示する
    plt.legend()
    plt.grid()
    plt.show()


if __name__ == '__main__':
    main()

これを実行すると,以下の実行結果が得られます(ただし,乱数が入っているので全く同じはならないです).

実行結果

|   iter    |  target   |     x     |
-------------------------------------
Try :  1  / 5 , next x is  11.259746731423292
Input y : -123.43686037691553
|  1        | -123.4    |  11.26    |
Try :  2  / 5 , next x is  10.17666985653311
Input y : -110.96666257730458
|  2        | -111.0    |  10.18    |
Try :  3  / 5 , next x is  1.893325699549166
Input y : 6.168510513881696
|  3        |  6.169    |  1.893    |
Try :  4  / 5 , next x is  -5.0
Input y : -28.457911018076135
|  4        | -28.46    | -5.0      |
Try :  5  / 5 , next x is  -0.9393494872633938
Input y : -7.494087232595987
|  5        | -7.494    | -0.9393   |
Try :  6  / 5 , next x is  3.8837562094488702
Input y : -14.850130945917755
|  6        | -14.85    |  3.884    |
Try :  7  / 5 , next x is  0.9116311089700233
Input y : -3.57363565519136
|  7        | -3.574    |  0.9116   |
Try :  8  / 5 , next x is  -2.8434605290266846
Input y : -14.069363342763157
|  8        | -14.07    | -2.843    |
Try :  9  / 5 , next x is  2.579734536274863
Input y : -0.6395105743239551
|  9        | -0.6395   |  2.58     |
Try :  10  / 5 , next x is  6.339581630124683
Input y : -31.424942315414583
|  10       | -31.42    |  6.34     |
Try :  11  / 5 , next x is  15.0
Input y : -224.12057472706056
|  11       | -224.1    |  15.0     |
Try :  12  / 5 , next x is  1.6599607452794116
Input y : -11.12890254409341
|  12       | -11.13    |  1.66     |
Try :  13  / 5 , next x is  7.965308786494606
Input y : -58.939492177353785
|  13       | -58.94    |  7.965    |
=====================================
x_max =  1.893325699549166 max_y =  6.168510513881696

緑のデータ点から予想される最大の点はx = 1.89のときy = 6.17であるという解が得られました.

以降で説明するように,データ点に加えた乱数(真の関数からの差)を結構大きくしているので,真の最適解x = 0, y = 0から結構外れてます.

ちなみに,探索点の順番をアニメーションにするとこんな感じです.

アニメーションを実行するには,コード中に

import matplotlib.animation as animation

および

def plot_ani(bo):
    fig = plt.figure()
    ax = fig.add_subplot(111)

    X = [x for x in np.arange(-5, 15, 0.1)]
    # 真の関数
    y = [-x**2 for x in X]
    plt.plot(X, y, label='true')

    #最大値
    max_x = bo.max['params']['x']
    max_y = bo.max['target']
    print('max = ', max_x, max_y)

    # サンプル点
    xs = [bo.res[p]['params']['x'] for p in range(len(bo.res))]
    ys = [bo.res[p]['target'] for p in range(len(bo.res))]

    line, = ax.plot([], [], '.g', markersize=15)

    def init():
        line.set_data([], [])
        return line

    def animate(i):
        line.set_data(xs[i],ys[i])
        ax.scatter(xs[i], ys[i] , c='black', s=15)
        return line

    ani = animation.FuncAnimation(fig, animate, len(bo.res), interval=500)

    plt.show()
    ani.save("bayes_ani_2D.gif", writer='pillow')

を加えて,def main():内の最後に

#アニメーション--------------
plot_ani(bo)

を加え,以下のように

# グラフを表示する
#plt.legend()
#plt.grid()
#plt.show()

とコメントアウトしてください.

プログラム解説

正解データ生成関数

def f(x):
    global try_

    print("Try : " , try_ , " / {0} , next x is ".format(exp_num) , x)
    #score = float(input("Input y : "))
    score = -x**2 + 10*random.uniform(-1,1)
    print('Input y :', score)

    try_ += 1

    return score

測定したデータ点を生成する部分です.xを入力すると,正解データ点が返されます.

ここでは,正解データは真の関数に乱数を加えたscore = -x*2 + 10*random.uniform(-1,1)で生成しました.乱数は±10の値を取るので,そこそこ-x*2から離れた値が生成されます.

scoreを好きな関数にすることで,アレンジできます.もちろん,

score = float(input("Input y : "))

として実験結果データなどを手入力しても構いません.

プロット関数

def plot_bo(bo):
    # プロット範囲 (決め打ち)
    X = [x for x in np.arange(-5, 15, 0.1)]

    # 真の関数
    y = [-x**2 for x in X]
    plt.plot(X, y, label='true')

    # サンプル点
    #print(bo.res[0])
    xs = [bo.res[p]['params']['x'] for p in range(len(bo.res))]
    ys = [bo.res[p]['target'] for p in range(len(bo.res))]
    plt.scatter(xs, ys, c='green', s=20, zorder=10, label='sample')

    # 予測結果
    mean, sigma = bo._gp.predict(np.array(X).reshape(-1, 1), return_std=True)
    plt.plot(X, mean, label='pred')  # 推定した関数
    plt.fill_between(X, mean + sigma, mean - sigma, alpha=0.1)  # 標準偏差

    # 最大値
    max_x = bo.max['params']['x']
    max_y = bo.max['target']
    print('max = ', max_x, max_y)
    plt.scatter([max_x], [max_y], c='red', s=50, zorder=10, label='pred_max')

単に,グラフを表示させるための関数です.プログラムの最後でplot_bo(bo)で呼び出します.

xs = [bo.res[p]['params']['x'] for p in range(len(bo.res))]

max_x = bo.max['params']['x']

の意味が分からない場合は,

#print(bo.res[0])

のコメントアウトを外して(#を消去して),bo.resの中身を見てみれば良いでしょう.

mean, sigma = bo._gp.predict(np.array(X).reshape(-1, 1), return_std=True)

とするだけで平均と分散を計算できるというのはありがたいですね.

ベイズ最適化

# 探索するパラメータと範囲を決める
pbounds = {
    'x': (-5, 15),
}
# 探索対象の関数と、探索するパラメータと範囲を渡す
bo = BayesianOptimization(f=f, pbounds=pbounds)
# 最大化する
bo.maximize(init_points=3, n_iter=10)

BayesianOptimizationの引数は(f=正解データ点を生成する関数,pbounds=探索範囲)とします.

pboundsは,辞書型にする必要があるようです.

bo.maximize(init_points=3, n_iter=10)のinit_pointsは初期観測点の数,n_iterは何点評価するかを意味します.n_iterを大きくしすぎると過学習になってしまうと思います.

ちなみに,maximizeしかなく,minimizeは用意されていないようです.minimizeにする方法はあるようですが,単純に関数にマイナスをつけておけばOKでしょう.

多変数関数のベイズ最適化

多変数関数と言っても,コードを少し変えるだけで実行できます.

プログラムは以下です.

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import math
import numpy as np
from bayes_opt import BayesianOptimization
from matplotlib import pyplot as plt
import random
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

#実験回数
exp_num = 5

try_ = 1

x1_max = 10
x1_min = -10
x2_max = 10
x2_min = -10

def f(x1, x2):
    global try_

    #print("Try : " , try_ , " / {0} , next x is ".format(exp_num) , x1)
    #score = float(input("Input y : "))
    score = - x1**2 - 3*x2**2 + random.uniform(-0.1,0.1)
    print('Input y :', score)

    try_ += 1

    return score

def plot_bo(bo):
    # プロット範囲 (決め打ち)
    X1 = [x1 for x1 in np.arange(x1_min, x1_max, 0.1)]
    X2 = [x2 for x2 in np.arange(x1_min, x2_max, 0.1)]
    X1, X2 = np.meshgrid(X1, X2)

    # 真の関数
    y = - X1**2 - 3*X2**2
    surf = ax.plot_surface(X1, X2, y, alpha=0.2, cmap=plt.cm.coolwarm,
                      linewidth=0, antialiased=False)
    #ax.plot(X1, X2, y, label='true')

    # カラーバーの表示
    fig.colorbar(surf, shrink=0.3, aspect=10)


    # サンプル点
    print(bo.res[0])
    xs1 = [bo.res[p]['params']['x1'] for p in range(len(bo.res))]
    xs2 = [bo.res[p]['params']['x2'] for p in range(len(bo.res))]
    ys = [bo.res[p]['target'] for p in range(len(bo.res))]
    ax.scatter(xs1,xs2, ys, c='black', s=20, zorder=10, label='sample')


    # 最大値
    print(bo.max)
    max_x1 = bo.max['params']['x1']
    max_x2 = bo.max['params']['x2']
    max_y = bo.max['target']
    print('max = ', max_x1, max_x2, max_y)
    for y in range(-200, 100):
        ax.scatter([max_x1],[max_x2], y , c='green', s=1)
    ax.scatter([max_x1],[max_x2],[max_y], c='green', s=100, zorder=10, label='pred_max')


def plot_ani(bo):
    # プロット範囲 (決め打ち)
    X1 = [x1 for x1 in np.arange(x1_min, x1_max, 0.1)]
    X2 = [x2 for x2 in np.arange(x1_min, x2_max, 0.1)]
    X1, X2 = np.meshgrid(X1, X2)
    # 真の関数
    y = - X1**2 - 3*X2**2
    surf = ax.plot_surface(X1, X2, y, alpha=0.2, cmap=plt.cm.coolwarm,
                      linewidth=0, antialiased=False)
    #最大値
    max_x1 = bo.max['params']['x1']
    max_x2 = bo.max['params']['x2']
    max_y = bo.max['target']
    print('x_max = ', max_x1, max_x2,'max_y = ', max_y)
    for y in range(-400, 200):
        ax.scatter([max_x1],[max_x2], y , c='green', s=1)
    # サンプル点
    print(bo.res[0])
    xs1 = [bo.res[p]['params']['x1'] for p in range(len(bo.res))]
    xs2 = [bo.res[p]['params']['x2'] for p in range(len(bo.res))]
    ys = [bo.res[p]['target'] for p in range(len(bo.res))]

    line, = ax.plot([], [], '.g', markersize=15)

    def init():
        line.set_data([], [])
        return line

    def animate(i):
        line.set_data(xs1[i],xs2[i])
        line.set_3d_properties(ys[i])
        ax.scatter([xs1[i]],[xs2[i]], ys[i] , c='black', s=15)
        return line

    ani = animation.FuncAnimation(fig, animate, len(bo.res), interval=500)

    plt.show()
    ani.save("bayes_ani.gif", writer='pillow')

def main():
    # 探索するパラメータと範囲を決める
    pbounds = {'x1': (x1_min, x1_max), 'x2': (x2_min, x2_max)}

    # 探索対象の関数と、探索するパラメータと範囲を渡す
    bo = BayesianOptimization(f=f, pbounds=pbounds)
    # 最大化する
    bo.maximize(init_points=3, n_iter=10)

    print('OPTIMIZE: END')

    # 結果をグラフに描画する
    plot_bo(bo)
    # グラフを表示する
    plt.legend()
    plt.grid()
    plt.show()

    #plt.cla()

    #アニメーション--------------
    #plot_ani(bo)


if __name__ == '__main__':
    main()

先ほどの一変数関数のときとほとんど同じなので,解説はいらないと思います.

これを実行すると,以下のような結果が得られます.

緑の棒が,ベイズ最適化によって得られた最適解です.黒の点はサンプル点(正解データ)です.

アニメーションを表示する場合は,

# 結果をグラフに描画する
    plot_bo(bo)
    # グラフを表示する
    plt.legend()
    plt.grid()
    plt.show()

の部分をコメントアウトして,

#plot_ani(bo)

のコメントアウトを外して(#を消して)実行してください.

その結果はこちらです.

アニメーションといっても,ただ探索してる場所を表示させているだけですが…(笑)

分散の表示は余裕があればやります.

コメント

タイトルとURLをコピーしました