ベイズ最適化で関数の最適解を探索【python実装】

2020年11月27日

今回は、ベイズ最適化によってAckley関数の最適解を探索し、任意のベイズ最適化の繰り返しの回数ごとに、ガウス過程回帰によってどのような関数推定が成されているのかを可視化します。

細かい説明はまた別日にするとして、今回はGPyでの実装コードと結果を述べる程度の記事です。

僕が書いたコードよりももっと良い方法がある。or 間違いがある等あれば、どんどん指摘してください!お願いします。

ベイズ最適化に用いる関数を定義

まず、ベイズ最適化で用いるメソッドを関数化しておきます。

環境はPython 3.5。

ライブラリは、

・GPy

・numpy

・scipy

"""
ベイズ最適化に用いるメソッドをまとめたファイル。

"""

import GPy
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm


def GP_Model(dim, data_x, data_y):
    """
    ガウス過程回帰のモデルを構築する関数

    param dim: 入力パラメータの次元数
    param data_x: 入力パラメータ
    param data_y: y

    return:
    """
    kern = GPy.kern.RBF(dim) + GPy.kern.White(dim)
    Gpy_model = GPy.models.GPRegression(data_x, data_y, kern, None)
    Gpy_model.optimize()
    return Gpy_model


def Produce_Best(purpose, value):
    """
    工程ごとの実験結果を入れれば、最適値を返してくれる関数

    param purpose: 最適化の目的。 最大化 or 最小化を想定。
    param value: 実験値。(最適化したい出力) 。 一次元を想定。

    return: 現状の実験値の中で最も最適な実験値を一つ返す。
    """
    if purpose == 'max':
        Best_value1 = np.max(value)
        return Best_value1
    if purpose == 'min':
        Best_value2 = np.min(value)
        return Best_value2


def EI_forMax(Best_value, model, x, purpose='max'):
    """
    最大値、または最小値を探索する時に用いる獲得関数 EI

    param Best_value: 現状で一番最適な値
    param model: GPモデル
    param x: 入力パラメータ。 ココにいれた入力値における獲得関数を返す。

    return: 指定した入力値での獲得関数EIの値。
    """
    def GP_predict(GPmodel, x_pre):
        mean_cal, var_cal = GPmodel.predict(x_pre)
        return mean_cal, var_cal

    mean, var = GP_predict(model, x)

    sigma = var ** .5
    imp = 0

    if purpose == 'max':
        imp = np.max(mean - Best_value, 0)
    elif purpose == 'min':
        imp = np.max(Best_value - mean, 0)

    if sigma != 0:
        Z = imp / sigma
        EI = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
    else:
        EI = 0.0

    return EI


def EI_forApproach(Best_value, model, x):
    """
    出力をある値に近づけるためのEI.

    param Best_value: 近づけたい値。(固定値)
    param model: GPモデル
    param x: 入力パラメータ。 ココにいれた入力値における獲得関数を返す。

    return: 指定した入力値での獲得関数EIの値。
    """
    def GP_predict(GPmodel, x_pre):
        mean_cal, var_cal = GPmodel.predict(x_pre)
        return mean_cal, var_cal

    mean, var = GP_predict(model, x)

    sigma = var ** .5
    Bad = np.max(abs(Best_value - mean))    # 求めたい値に近い程、不良度が下がる。

    if sigma != 0:
        Z = Bad / sigma
        EI = Bad * norm.cdf(Z) + sigma * norm.pdf(Z)
    else:
        EI = 0.0 

    return EI


def Optimization(acquisition, Best_value, model, seed, n_restart, dim, bounds_of_input, purpose='max'):
    """
    入力された関数を準ニュートン法で最適化し、最適解の入力値を返す関数 (獲得関数の最適化 & 次の実験値の導出に使用).
    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)

    def MinObject_approach(xx):
        return acquisition(Best_value, model, xx.reshape(-1, dim))

    if acquisition == EI_forMax:
        min_val = 1
        min_x = None

        for x0 in np.round(seed * np.random.rand(n_restart, dim), 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

    elif acquisition == EI_forApproach:
        min_val = 1000000
        min_x = None
        for x0 in np.round(seed * np.random.rand(n_restart, dim), 0):
            res = minimize(MinObject_approach, 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

コードの中でも言及していますが、

・GP_Model : ガウス過程回帰のモデルを構築

・Produce_Best : 現状での最適値を返す関数(最大化、最小化のベイズ最適化で使用)

・EI_forMax : 最大化(または最小化)のための獲得関数EI

・EI_forApproach : 目的関数をある値に漸近させるための獲得関数EI(みたいなやつ)

・Optimization : 獲得関数を最適化する関数

これらの関数を使って、関数の最適解を探索します。

Ackley 関数について

推定する関数は、Ackley functionと呼ばれる、有名なベンチマーク関数。

可視化してみると、以下のような感じになります。

これは、プロットのメッシュ間隔を 1 にした時。

これを見た感じ、かなり単純な形に見えます。

しかし実際はかなり細かくギザギザとしている関数であり、局所解が複数存在します。

メッシュを二倍に細かくした図が以下。

ベイズ最適化による最適解の探索

ベイズ最適化によってAckley関数の最小値(最適解)を探索していきます。

また、その過程でガウス過程回帰がどのような予測分布を出力しているのかを可視化します。

コードは以下。

可視化にはmatplotlibを使用。

一番初めに紹介したファイルを「BO_Method」という名前にし、冒頭で中身の関数を取り出しています。

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 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)
        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)
# inputを保存したいときはここのコメントアウトを外す
# df_input = pd.DataFrame(x)
# df_input.to_csv('C:/卒研関係/Outputs/BO-3D-revise/all_inputs.csv', header=False, index=False)

・Optimization を、BO_Methodで書いたものから便宜的に書き直しています。
(Ackley関数を(-30 ≦ x ≦ 30)で最適化するため。)

・初期点は20個で、ベイズ最適化の繰り返しの回数は、40回としている。

・matplotlibのimshowで軸がちゃんと表示されないところで、結構ハマりました。
普通にmeanの値をmesh状にして可視化するだけだと、軸が[0, 60]になってしまうんですね。
軸の数値を書き換えたり色々試しましたが、今回のように、ax.imshow の extent=[]をイジるのがベターのようです。

・(修正)メッシュグリッドを生成する時、
x1_mesh, x2_mesh = np.meshgrid(x1, x2)
としていましたが、これだと描画したい図に対してX2が上下逆の関係性で出力されていました。
なので、
x1_mesh, x2_mesh_inversion = np.meshgrid(x1, x2)
x2_mesh = np.flipud(x2_mesh_inversion)
として、上下関係を反転させ、目的にあったメッシュに直しています。

出力結果

gp_heat_mapという関数を使って、ベイズ最適化繰り返しの2回ごとに推定した関数形のヒートマップを出力しました。

たとえば、繰り返し6回目の推定関数は以下のような感じ。

推定関数のヒートマップ

大まかな形は予測できています。

もっと探索を進めていくと、Ackley関数の「ギザギザ具合」を学習し始め、濃淡のブチがはっきりしていきます。

繰り返し40回目だと以下の感じ。

推定関数の評価関数

ベイズ最適化の繰り返しにおける、関数推定の様子を可視化できました。

まとめ

次回はヒートマップに探索点のプロットを追加し、最適解探索の様子を可視化したいと思います。

また、別の機会に、このベイズ最適化が私の専攻である材料工学でどのように活用されているのかも紹介したいと思います。

では。

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

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

Python

Posted by arabica