physboを使って多目的最適化問題をベイズ最適化で解いてみた

• 製造業出身のデータサイエンティストがお送りする記事
• 今回はphysboを使って多目的最適化問題をベイズ最適化で解いてみました。

はじめに

過去に多目的最適化問題については、記事を何個か書いておりますので参考にして頂けますと幸いです。

• 多目的最適化問題のNSGA-Ⅱを勉強したの整理しました
• 多目的最適化問題のNSGA-Ⅱを実装してみた
• 多目的最適化問題のNSGA-Ⅲを勉強したの整理しました
• 多目的最適化問題のNSGA-Ⅲを実装してみた
• pymooを使って多目的最適化のNSGA-Ⅱを実装してみた

使用するライブラリー（physbo）

今回は最適化ライブラリphysboを使って実装したいと思います。

• Officail Page
• MateriApps
• マニュアル

physboの実装

今回は実装サンプルを試してみました。

# 必要なライブラリーのインストール
import numpy as np
import matplotlib.pyplot as plt
import physbo
import itertools
%matplotlib inline

次にテスト関数を作成します。

def vlmop2_minus(x):
    n = x.shape[1]
    y1 = 1 - np.exp(-1 * np.sum((x - 1/np.sqrt(n)) ** 2, axis = 1))
    y2 = 1 - np.exp(-1 * np.sum((x + 1/np.sqrt(n)) ** 2, axis = 1))

    return np.c_[-y1, -y2]

探索範囲を指定します。

a = np.linspace(-2, 2, 101)
test_X = np.array(list(itertools.product(a, a)))

シミュレータを定義します。

class simulator(object):
    def __init__(self, X):
        self.t = vlmop2_minus(X)

    def __call__( self, action):
        return self.t[action]

simu = simulator(test_X)

今回の設定した関数は、1つは右上にピークがあり、1つは左下にピークがある関数です。
可視化すると以下のような感じです。

# plot objective 1
plt.figure()
plt.imshow(simu.t[:,0].reshape((101,101)),
           vmin=-1.0,
           vmax=0.0,
           origin="lower",
           extent=[-2.0, 2.0, -2.0, 2.0]
          )
plt.title("objective 1")
plt.colorbar()
plt.plot([1.0/np.sqrt(2.0)],
         [1.0/np.sqrt(2.0)],
         '*'
        )
plt.show()

# plot objective 2
plt.figure()
plt.imshow(simu.t[:,1].reshape((101,101)),
           vmin=-1.0,
           vmax=0.0,
           origin="lower",
           extent=[-2.0, 2.0, -2.0, 2.0]
          )
plt.title("objective 2")
plt.colorbar()
plt.plot([-1.0/np.sqrt(2.0)],
         [-1.0/np.sqrt(2.0)],
         '*'
        )
plt.show()

次に最適化を実行していきます。
* policyのセット
* 多目的最適化用のphysbo.search.discrete_multi.policyを利用
* num_objectivesに目的関数の数を指定

policy = physbo.search.discrete_multi.policy(test_X=test_X,
                                             num_objectives=2
                                            )
policy.set_seed(0)

ベイズ最適化は、下記3パターンの方法があるそうです。
* HVPI (Hypervolume-based Probability of Improvement)
* EHVI (Expected Hyper-Volume Improvement)
* TS (Thompson Sampling)

今回はHVPIを使用してみたいと思います。

HVPI (Hypervolume-based Probability of Improvement)は、多次元の目的関数空間における非劣解領域 (non-dominated region) の改善確率をscoreとして算出する手法だそうです。

policy = physbo.search.discrete_multi.policy(test_X=test_X,
                                             num_objectives=2
                                            )
policy.set_seed(0)

policy.random_search(max_num_probes=10,
                     simulator=simu
                    )
res_HVPI = policy.bayes_search(max_num_probes=40,
                               simulator=simu,
                               score='HVPI',
                               interval=10
                              )

パレート解の取得は以下のコードでできます。

front, front_num = res_HVPI.export_pareto_front()
front, front_num

結果はこんな感じで返ってきます。
各目的変数の値と対象のインデックス番号が得られます。

(array([[-9.82978321e-01, -3.32414921e-04],
        [-9.71056940e-01, -3.45637262e-02],
        [-9.50786244e-01, -6.76228653e-02],
        [-9.34458810e-01, -1.15668709e-01],
        [-9.13824628e-01, -1.71906452e-01],
        [-8.88492700e-01, -2.36876025e-01],
        [-8.68576210e-01, -2.82707709e-01],
        [-8.57549341e-01, -3.05692550e-01],
        [-8.33710080e-01, -3.53626320e-01],
        [-8.20907701e-01, -3.78330099e-01],
        [-7.92933832e-01, -4.26778524e-01],
        [-7.77704642e-01, -4.50447880e-01],
        [-7.62117827e-01, -4.74823143e-01],
        [-7.45438101e-01, -4.98117246e-01],
        [-7.28459156e-01, -5.21910482e-01],
        [-7.14031659e-01, -5.50368010e-01],
        [-6.92015440e-01, -5.67552511e-01],
        [-6.72522879e-01, -5.89370083e-01],
        [-6.52909089e-01, -6.11332656e-01],
        [-6.32120559e-01, -6.32120559e-01],
        [-6.11332656e-01, -6.52909089e-01],
        [-5.89370083e-01, -6.72522879e-01],
        [-5.70311337e-01, -6.93980247e-01],
        [-5.44575729e-01, -7.10347737e-01],
        [-5.21910482e-01, -7.28459156e-01],
        [-4.98117246e-01, -7.45438101e-01],
        [-4.74823143e-01, -7.62117827e-01],
        [-4.26778524e-01, -7.92933832e-01],
        [-4.02089723e-01, -8.07119688e-01],
        [-3.90152305e-01, -8.24313473e-01],
        [-3.53626320e-01, -8.33710080e-01],
        [-3.05692550e-01, -8.57549341e-01],
        [-2.58961961e-01, -8.78749507e-01],
        [-2.16642598e-01, -8.97780597e-01],
        [-1.71906452e-01, -9.13824628e-01],
        [-1.15668709e-01, -9.34458810e-01],
        [-6.76228653e-02, -9.50786244e-01],
        [-3.51124618e-02, -9.67608404e-01],
        [-1.46847591e-03, -9.78680443e-01]]),
 array([21, 15, 35, 26, 32, 18, 43, 27, 40, 16, 31, 49, 23, 47, 33, 12, 34,
        45, 22, 44, 28, 41, 13, 42, 29, 48, 19, 30, 46, 10, 37, 24, 38, 17,
        36, 25, 39, 11, 20]))

パレート解の可視化は以下のような感じでできます。

def plot_pareto_front(res):
    front, front_num = res.export_pareto_front()
    dominated = [i for i in range(res.num_runs) if i not in front_num]
    points = res.fx[dominated, :]

    plt.figure(figsize=(7, 7))
    plt.scatter(res.fx[dominated,0], res.fx[dominated,1], c = "blue")
    plt.scatter(front[:, 0], front[:, 1], c = "red")
    plt.title('Pareto front')
    plt.xlabel('Objective 1')
    plt.ylabel('Objective 2')
    plt.xlim([-1.0, 0.0])
    plt.ylim([-1.0, 0.0])

plot_pareto_front(res_HVPI)

ちなみに全探索もできます。

policy = physbo.search.discrete_multi.policy(test_X=test_X,
                                             num_objectives=2
                                            )
policy.set_seed(0)

N = test_X.shape[0]
res_all = policy.random_search(max_num_probes=N,
                               simulator=simu
                              )

plot_pareto_front(res_all)

さいごに

最後まで読んで頂き、ありがとうございました。
今回は、physboを使って多目的最適化問題をベイズ最適化を使って解きました。
今回は単純な問題だったので時間はかかりませんでしたが、複雑な問題となった場合、実行時間がどこまでかかるのかは気になりました。
時間がある時に試してみようと思います。

訂正要望がありましたら、ご連絡頂けますと幸いです。