GPyとScipyで実装するベイズ最適化 -最適化プロセスの可視化-

2020年1月7日

前回の記事で、

・ベイズ最適化に使用するメソッドの作成
・Ackley関数の最適解を探索していく中で、ガウス過程回帰によってなされている関数推定を可視化

を行いました。

今回は、前回出力した関数推定のヒートマップに、

・今までのデータ点
・ベイズ最適化によって提案された次のデータ点

をプロットし、Ackley関数の最適解の探索がどのように進んでいるのかを可視化します。

可視化には前回と同じくmatplotlibを用います。

予測平均のヒートマップに、探索点のプロットを追加。

from BO_Method import GP_Model, EI_forMax, Produce_Best
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import pandas as pd


def f(in_x):
    """
    ノイズ入り、Ackley_function
    param in_x: 入力パラメータ。二次元。
    return: 出力。
    """
    in_x = np.array(in_x).reshape(-1, 2)
    in_x1, in_x2 = np.split(in_x, 2, 1)
    # Ackley関数
    z = 20 - 20 * np.exp(-0.2 * np.sqrt((in_x1 ** 2 + in_x2 ** 2) / 2)) + np.exp(1) - np.exp((np.cos(2 * np.pi * in_x1)
                                                                                              + np.cos(
                2 * np.pi * in_x2)) / 2) + np.random.normal(0, .01)
    z = np.array(z).reshape(-1, 1)
    return z


def produce_mesh_grid(input_bound_low, input_bound_high, interval):
    """
    input_bound_low : grid bounds(low)
    input_bound_high : grid bounds(high)
    return: mesh grid value.
    """
    x1 = np.arange(input_bound_low, input_bound_high, interval)
    x2 = np.arange(input_bound_low, input_bound_high, interval)
    x1_mesh, x2_mesh_inversion = np.meshgrid(x1, x2)   # 描画したい図に対して、X2のグリッドが上下逆。
    x2_mesh = np.flipud(x2_mesh_inversion)
    return x1_mesh, x2_mesh


def gp_heat_map(model_of_GP, iteration_num, Produce_mesh):
    """
    ガウス過程回帰の予測平均のヒートマップを生成。
    param model_of_GP: model of GP that is already learned.
    param iteration_num: iteration number of Bayesian Optimization.
    return: Single graph with colormap.
    """

    X1, X2 = Produce_mesh(input_bound_low=-30, input_bound_high=30, interval=1)

    X1_1d = X1.ravel().reshape(1, -1)
    X2_1d = X2.ravel().reshape(1, -1)
    X_grid_values = np.vstack((X1_1d, X2_1d)).transpose()
    mean, var = model_of_GP.predict(X_grid_values)

    pred = np.array(mean).reshape(X1.shape)

    figure, axes = plt.subplots(1, 1, figsize=(6, 6))
    pos = axes.imshow(pred, cmap='Blues', interpolation='none', extent=[-30, 30, -30, 30])
    axes.set_title("{}th_GP_prediction's_colormap".format(iteration_num))
    axes.set_xlim(-30, 30)
    axes.set_ylim(-30, 30)
    cb = figure.colorbar(pos)
    cb.set_label('prediction_of_Output')

    return figure, axes


def add_scatter(ax_exist, opt_x, prior_x):
    """
    すでに存在するfigure(GPの予測平均分布を想定)中の axes に、プロットを追加する関数。

    param ax_exist: すでに存在しているaxes
    param opt_x: 次の提案値
    param prior_x: これまでの実験値
    return: 点を追加したaxesを返す
    """
    x1_opt = opt_x[0]
    x2_opt = opt_x[1]
    scale_opt = 50
    scale_prior = 30

    if prior_x is None:
        ax_exist.scatter(x1_opt, x2_opt, c='white', s=scale_opt, label='opt_x', edgecolors='black')
        return ax_exist

    else:
        x1_prior, x2_prior = np.hsplit(prior_x, 2)
        ax_exist.scatter(x1_prior, x2_prior, c='blue', s=scale_prior, label='prior_x', edgecolors='white')

        ax_exist.scatter(x1_opt, x2_opt, c='red', s=scale_opt, label='opt_x', edgecolors='none')
        ax_exist.legend()
        return ax_exist


def Optimization(acquisition, Best_value, model, seed, Translation, n_restart, dim, bounds_of_input, purpose='max'):
    """
    BO_methodにあるOptimizationに、Translationの引数を新たに加えた。
    (準ニュートン法の初期値を、今回のAckley関数の最適化問題にfitさせるため。)

    入力された関数を準ニュートン法で最適化し、最適解の入力値を返す関数 (獲得関数を最適化し、提案値を導出する).
    bounds_of_input must be tuple objects.

    param acquisition:
          最適化したい値を生成する関数
    param Best_value:
          現状の最適値(最大化 or 最小化 の場合)。 または近づけたい値。
    param model:
          GP model
    param seed:
          整数。 np.random.randで、0.0以上、1.0未満の最適化における初期点を生成するが、seedをかけることで、そのスケールを制御。
    param n_restart:
          準ニュートン法では局所解にたどり着きにくいように、最適化計算を何度か試行する。 その繰り返しの回数を指定。
    param dim:
          入力パラメータの次元数
    param bounds_of_input:
          入力パラメータの定義域を、タプルで指定。

    return:  最適解を出力するような入力値を返す。
    """

    def MinObject_max(xx):
        return (-1) * acquisition(Best_value, model, xx.reshape(-1, dim), purpose)

    min_val = 1
    min_x = None

    for x0 in np.round(seed * np.random.rand(n_restart, dim) - Translation, 0):
        res = minimize(MinObject_max, x0=x0, bounds=bounds_of_input, method='L-BFGS-B')
        if res.fun < min_val:
            min_val = res.fun
            min_x = res.x

    return min_x


input_dim = 2
output_dim = 1

# 初期点20個。 20 × 2 の行列. -30 < x < 30 の乱数。
x = np.round(60 * np.random.rand(20, 2), 0) - 30
print('init_x :', x)

y = f(x)

Iter = 40
# ベイズ最適化での最適値を補完するリスト
A = []

#  iteration of BO

for i in range(Iter):
    print("_______________roop %d th________________" % (i + 1))
    # それぞれのベストを計算
    Best = Produce_Best('min', y)
    Best = np.round(Best, 3)
    A.append(Best)
    print("Best", Best)
    # ___________________BO_____________________
    GPM = GP_Model(input_dim, x, y)

    # 全入力の範囲が同じなら、最適化の入力の数を入力するだけでOKになる
    bounds_generate = tuple([(-30.0, 30.0) for i in range(input_dim)])

    # Optimization EI
    propose_x = Optimization(EI_forMax, Best, GPM, seed=60, Translation=30, n_restart=1000, dim=input_dim,
                             bounds_of_input=bounds_generate, purpose='min')
    print('propose_x:', propose_x)

    draw_span = 1  # Iteration何回ごとに描画するか。

    # ============予測分布のヒートマップと実験値の散布図のプロット==========
    if (i + 1) % draw_span == 0:
         fig1, ax1 = gp_heat_map(GPM, i + 1, produce_mesh_grid)
         ax1 = add_scatter(ax1, propose_x, x)
         fig1.savefig('C:/Outputs/BO-3D-pred_plot/{}th_prediction'.format(i + 1), dpi=150)

    # next examination
    y_next = f(propose_x)

    # データを更新
    x = np.vstack((x, propose_x))
    y = np.vstack((y, y_next))

Best_final = np.min(y)
Best_final = np.round(Best_final, 3)
A.append(Best_final)

print('Best_final', Best_final)
print('A', A)

# df_input = pd.DataFrame(x)
# df_input.to_csv('C:/Outputs/BO-3D-revise/all_inputs.csv', header=False, index=False)

・初期点20個。ベイズ最適化のIterationは40回。

・add_scatterという関数を新たに追加。
それまでのデータ点を青色の点、次の提案値を赤色の点としてプロットしています。

出力結果

ベイズ最適化の繰り返しの数に対応した、ヒートマップとプロットが出力されます。

繰り返し二回目だと以下。

右下が小さい(右下に最適解がありそう)ような予測がされており、次の提案値として右下の赤丸が提案されています。

さらに探索を続けて、9回目だと以下。

大まかな関数形が予測されており、真ん中あたりに最適解があると予測しています。

そしてここらへん(9回目以降)から、真ん中あたりを重点的に探索し始めます。

探索の様子をgif画像にする

イテレーションごとの最適化の様子を出力できたので、最適解探索の様子をより見やすくするために、gifにしてやりましょう。

コードは以下。外部ライブラリのpillowを使用します。

from PIL import Image

IMAGE_DIR = "C:/Outputs/BO-3D-pred_plot/"

X = []

for i in range(40):
    img_name = IMAGE_DIR + '{}th_prediction.png'.format(i + 1)
    img = Image.open(img_name)
    X.append(img)

X[0].save('C:/Outputs/BO-3D-pred_plot/BO.gif', save_all=True,
          append_images=X[1:], optimize=False, duration=200, loop=0)

・IMAGE_DIRでは、画像を保存した一つ上のディレクトリのパスを保存しておきます。

・.saveで、一枚目の画像にそれ以降の画像を連番で保存しgif化しています。
durationの数値は一枚の画像の表示時間(ミリ秒)。
もっと大きくすれば、より「ゆっくりなgif」になります。

出力結果

・左上に表示されている数字がベイズ最適化の繰り返しの数です。
最初の方はまだ探索されていない端っこの方を調べ、関数形の予測がある程度定まってからは、最適解の存在する(0, 0)付近を重点的に探索しています。

・最後の方は一見同じところを探索しているように見えますが、小数点以下のオーダーで違う位置を探索しています。
gif上でも良く見ると、探索位置や予測分布が微妙に変化しているのが確認できます。

まとめ

今回は、ベイズ最適化による最適解の探索の様子を可視化しました。
次回は、「EI以外の獲得関数を用いた最適化」か「ベイズ最適化の各イテレーションにおける獲得関数の可視化」を行いたいと思います。

では。

(コードが動かない・コードがおかしい・コードのここが分からない等ありましたら、コメントお願いします。)

↓1日1タップで応援よろしくお願いします!

にほんブログ村 にほんブログ村へ

Python

Posted by arabica