Python ベイズ最適化による単一目的最適化「GPyOpt」

　本記事では、目的関数Yの最小値または最大値と、説明変数Xを求めるための雛形コードを載せました。ベイズ最適化のPythonライブラリは複数あります。本記事ではGPyOptについてで、インストールはcondaで次のようにします。

conda install -c conda-forge gpyopt

 　はじめに、最適化する出力指標を目的関数と呼び、入力因子のことを説明変数と言います。例えば、下図のようなz=x(exp(-x/5)^2-(y/5)^2)の場合は、zが目的関数でxとyが説明変数です。このzの最大値を求めたいとします。期待する結果は(x, y)=(3.5前後, 0)で最大値は2前後です。（参考までに説明変数は3つ以上でも出来ます。また制約条件を設けることも可能です。本コード中にコメントアウト文で例を記載しています）

本コードを実行すると、(x, y)の探索履歴を下図のようにグラフで出力します。探索回数が増えるにつれて、〇の色の濃淡が濃くなるように表示しています。下図3つの内の一番上図は、9回探索した結果で(x, y)=(3.5, 0)近傍を多く探索している様子です。その次の中央の図は、19回探索した結果で(x, y)=(3.5, 0)付近にプロットがあり最大値を見つけたようです。一番下の図は、一応念のため他の箇所（x, yの各々の探索範囲の上下限）も探索している様子です。ここまですると普通にDOEしているのと変わらないような状況です。

そして、下図左のように最適結果、下図右のように探索履歴をcsvファイルに保存します。

参考までに、Jupyter Lab等での使用環境の場合、探索結果は下図のようにたった一行で可視化できます。

■本プログラム
　本コードは、Jupyter labで作成後に.pyファイルに出力して一部編集したものです。Anaconda PromptのようなCUIで実行できます。

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np 
import GPy
import GPyOpt

# 目的関数
def f(x):
    x0, x1 = x[:,0], x[:,1]
    return x0*(np.exp(-1 * pow(x0/5, 2) - pow(x1/5, 2)))

# 説明変数xの上下限の設定
bounds = [
# 連続値で指定したい場合は、typeをcontinuousとして数値は（下限, 上限）で設定する
{'name': 'x0', 'type': 'continuous', 'domain': (-10, 10)},
# 非連続値で指定したい場合は、typeをdiscreteとして数値を直打ちで設定する
{'name': 'x1', 'type': 'discrete', 'domain': (-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10)} 
]

''' 
# 制約条件を設定したい場合
constraints = [
{
"name": "constr_1",
"constraint": "(x[:,0] + x[:,1]) - 10 - 0.001" # x0 + x1 -10 <= 0.001 の制約設定したい場合
},
{
"name": "constr_2",
"constraint": "10 - (x[:,0] + x[:,1]) - 0.001" # x0 + x1 -10 >= -0.001 の制約設定したい場合
}
]
'''
flag_maximize = True # 目的関数を最大化したい場合にTrue。最小化したい場合はFalse
num = 30 # 探索回数（イタレーション）
myBopt = GPyOpt.methods.BayesianOptimization(
    f=f,
    domain=bounds,
    #constraints = constraints, # 制約条件を設定したい場合
    maximize = flag_maximize, # when True -f maximization of f is done by minimizing -f (default, False).
    model_type = 'GP' # default 'GP'
)
# 最適化を実行する
myBopt.run_optimization(max_iter = num)
# 最適値（最小値or最大値）の抽出
if flag_maximize:
    opt_fx = myBopt.fx_opt * (-1)
else:
    opt_fx = myBopt.fx_opt
# 最適値（最小値or最大値）の時の説明変数xの値の抽出
opt_x = myBopt.x_opt
print(opt_fx)
print(opt_x)


# In[2]:


# 最適化結果をcsvファイルに保存する
columns_name_list = ['x0', 'x1'] # 列名を作成
columns_name_list_str = ','.join(columns_name_list) # 要素間にカンマを入れて文字列へ
columns_name_list_str = columns_name_list_str + ',Y' + '\n' # ターゲット名と改行コードを入れる
with open('opt_data.csv', 'w') as f:
    f.write(columns_name_list_str)
    buf_list = []
    for X in opt_x.tolist():
        buf_list.append(X)
    buf_list.append(opt_fx)
    opt_list = list(map(str, buf_list))
    opt_str = ','.join(opt_list)
    f.write(opt_str + '\n')


# In[3]:


# 説明変数xの探索履歴をグラフで可視化する
import matplotlib.pyplot as plt
import japanize_matplotlib

step = int(num/3)
for i in range(step-1, num, step):
    history_X = myBopt.X[0:i]
    plt.rcParams['font.size'] = 13
    plt.figure()
    plt.title('開始から' + str(i) + '番目までの探索履歴')
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.xlim(-11, 11)
    plt.ylim(-11, 11)
    sc = plt.scatter(
        history_X[:,0], # X
        history_X[:,1], # Y
        s = 100,
        c = range(len(history_X)),
        cmap = "Greys")
    plt.colorbar(sc)
    plt.plot(history_X[:,0], history_X[:,1], linestyle=":", color='k')
    plt.grid()
    #plt.show()
    plt.savefig('history_' + str(i) + '.png')
    plt.close()


# In[4]:


# 探索履歴をcsvファイルに保存する
if flag_maximize:
    history_Y = myBopt.Y[0:num] * (-1)
else:
    history_Y = myBopt.Y[0:num]
history_Y_list = history_Y.tolist()
history_X_list = history_X.tolist()
with open('history_data.csv', 'w') as f:
    f.write(columns_name_list_str)
    for X, Y in zip(history_X_list, history_Y_list):
        X.extend(Y)
        buf_list = list(map(str, X))
        buf_str = ','.join(buf_list)
        f.write(buf_str + '\n')


# In[5]:


# 説明変数が1次元、2次元の場合の探索履歴は、次のように1行で書くこともできる。
#myBopt.plot_acquisition()
# 但し、最大化（maximize=True）に設定して解析した場合には、脳内で-1を掛ける必要がある。


# In[6]:


# 探索履歴は、次のように1行で書くこともできる。
#myBopt.plot_convergence()
# 但し、最大化（maximize=True）に設定して解析した場合には、右図は脳内で-1を掛ける必要がある。


# In[ ]:

以上

<広告>

リンク