この jupyter notebook ファイルは ISSP Data Repository (develop branch) から入手できます。
多目的最適化
最適化したい目的関数が複数 (\(p\) 個) ある場合は、多目的最適化をおこないます。
本チュートリアルでは、「解」は目的関数の組 \(y = (y_1(x), y_2(x), \dots, y_p(x))\) を意味することに注意してください。
解の大小関係 \(\prec\) を以下のように定義します。
\(y \prec y^{'}\Longleftrightarrow \forall \ i \le p, y_i \le y^{'}_i \land \exists \ j \le p, y_j < y^{'}_j\)
（最大化問題における） パレート解 とは、上記の大小関係の上で、自身よりも大きな解がないような解を指します。
すなわち、任意の目的関数の値を改善しようとした場合、他の目的関数のうちどれかひとつは悪化するような解です。
目的関数間にトレードオフが存在する場合は、パレート解は複数存在するため、それらを効率的に求めることが課題となります。
PHYSBOでは、パレート解を効率的に求めるためのベイズ最適化手法を実装しています。
[1]:
import time
import numpy as np
import matplotlib.pyplot as plt
import physbo
%matplotlib inline
seed = 12345
num_random_search = 10
num_bayes_search = 40
テスト関数
本チュートリアルでは、多目的最適化のベンチマーク関数である VLMOP2 を利用します。
目的関数の数は2つです。
\[\begin{split}\begin{split}
y_1(\vec{x}) &= 1 - \exp\left[-\sum_{i=1}^N\left( x_i - 1/\sqrt{N}\right)^2\right] \\
y_2(\vec{x}) &= 1 - \exp\left[-\sum_{i=1}^N\left( x_i + 1/\sqrt{N}\right)^2\right]
\end{split}\end{split}\]
\(y_1\) と \(y_2\) はそれぞれ \(x_1 = x_2 = \cdots x_N = 1/\sqrt{N}\) と \(x_1 = x_2 = \cdots x_N = -1/\sqrt{N}\) に最小値があり、その値はともに０です。また、上界は1 です。
PHYSBO は最大化問題を仮定するため、-1を掛けたものをあらためて目的関数とします。
参考文献
Van Veldhuizen, David A. Multiobjective evolutionary algorithms: classifications, analyses, and new innovations. No. AFIT/DS/ENG/99-01. AIR FORCE INST OF TECH WRIGHT-PATTERSONAFB OH SCHOOL OF ENGINEERING, 1999.
[2]:
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]
探索候補データの準備
入力空間 \(\vec{x}\) は 2次元とし、[-2, 2] × [-2, 2] の上で候補点をグリッド状に生成します。
[3]:
min_X = np.array([-2.0, -2.0])
max_X = np.array([2.0, 2.0])
test_X = physbo.search.utility.make_grid(min_X=min_X, max_X=max_X, num_X=101)
[4]:
test_X
[4]:
array([[-2. , -2. ],
[-2. , -1.96],
[-2. , -1.92],
...,
[ 2. , 1.92],
[ 2. , 1.96],
[ 2. , 2. ]], shape=(10201, 2))
[5]:
test_X.shape
[5]:
(10201, 2)
simulator の定義
[6]:
simu = physbo.search.utility.Simulator(test_X, vlmop2_minus)
関数のプロット
2つの目的関数をそれぞれプロットしてみましょう。
1つ目の目的関数は右上にピークがあり、2つ目の目的関数は左下にピークがあるようなトレードオフがある状態となっています。 （星はピークの位置です）
1つ目の目的関数
[7]:
t_all = vlmop2_minus(test_X)
plt.figure()
plt.imshow(t_all[:,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()
2つ目の目的関数
[8]:
# plot objective 2
plt.figure()
plt.imshow(t_all[:,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 に目的関数の数を指定してください。[9]:
policy = physbo.search.discrete_multi.Policy(test_X=test_X, num_objectives=2)
policy.set_seed(seed)
通常の
physbo.search.discrete.policy (目的関数が1つの場合）と同じく、random_search または bayes_search メソッドを呼ぶことで最適化を行います。基本的なAPIや利用方法は
discrete.policy とおおよそ共通しています。ランダムサーチ
今回のチュートリアルでは、10 (num_random_search) 回のランダムサンプリングと40 (num_bayes_search) 回のベイズ最適化を行います。 比較のために、50 (num_random_search + num_bayes_search) 回のランダムサンプリングも行っておきます。
[10]:
policy = physbo.search.discrete_multi.Policy(test_X=test_X, num_objectives=2)
policy.set_seed(seed)
res_random = policy.random_search(max_num_probes=num_random_search + num_bayes_search, simulator=simu, is_disp=True)
0001-th step: f(x) = [-0.99017418 -0.09267839] (action = 2560)
Pareto set updated.
the number of Pareto frontiers = 1
0002-th step: f(x) = [-0.63368874 -0.93288369] (action = 4415)
Pareto set updated.
the number of Pareto frontiers = 2
0003-th step: f(x) = [-0.89500547 -0.99942328] (action = 10046)
0004-th step: f(x) = [-0.52496048 -0.73019147] (action = 5205)
Pareto set updated.
the number of Pareto frontiers = 2
0005-th step: f(x) = [-0.06210869 -0.9918999 ] (action = 7542)
Pareto set updated.
the number of Pareto frontiers = 3
0006-th step: f(x) = [-0.78149934 -0.9714881 ] (action = 4018)
0007-th step: f(x) = [-0.99988434 -0.7158834 ] (action = 2631)
0008-th step: f(x) = [-0.46814364 -0.99957312] (action = 8063)
0009-th step: f(x) = [-0.58802709 -0.95199271] (action = 7519)
0010-th step: f(x) = [-0.64563383 -0.99949924] (action = 9458)
0011-th step: f(x) = [-0.85901111 -0.77832121] (action = 3396)
0012-th step: f(x) = [-0.6258165 -0.91403294] (action = 4413)
0013-th step: f(x) = [-0.96339156 -0.50605126] (action = 2577)
Pareto set updated.
the number of Pareto frontiers = 4
0014-th step: f(x) = [-0.99433124 -0.07837382] (action = 2955)
Pareto set updated.
the number of Pareto frontiers = 5
0015-th step: f(x) = [-0.77709331 -0.99991896] (action = 7370)
0016-th step: f(x) = [-0.44837296 -0.80073444] (action = 5909)
Pareto set updated.
the number of Pareto frontiers = 6
0017-th step: f(x) = [-0.86544619 -0.41433595] (action = 3787)
Pareto set updated.
the number of Pareto frontiers = 6
0018-th step: f(x) = [-0.56161131 -0.97933576] (action = 8027)
0019-th step: f(x) = [-0.99973799 -0.54162892] (action = 2034)
0020-th step: f(x) = [-0.99984173 -0.61123124] (action = 1231)
0021-th step: f(x) = [-0.97820847 -0.8330008 ] (action = 6382)
0022-th step: f(x) = [-0.99526434 -0.9637083 ] (action = 7282)
0023-th step: f(x) = [-0.97787893 -0.26212082] (action = 2569)
Pareto set updated.
the number of Pareto frontiers = 7
0024-th step: f(x) = [-0.24330703 -0.9118224 ] (action = 5619)
Pareto set updated.
the number of Pareto frontiers = 8
0025-th step: f(x) = [-0.24153475 -0.99178622] (action = 8040)
Pareto set updated.
the number of Pareto frontiers = 9
0026-th step: f(x) = [-0.84190968 -0.91981439] (action = 3506)
0027-th step: f(x) = [-0.99796446 -0.34753927] (action = 3349)
0028-th step: f(x) = [-0.99708276 -0.9919242 ] (action = 8291)
0029-th step: f(x) = [-0.99920018 -0.98099933] (action = 172)
0030-th step: f(x) = [-0.99665356 -0.14456179] (action = 2751)
0031-th step: f(x) = [-0.81366367 -0.99993226] (action = 7170)
0032-th step: f(x) = [-0.54951231 -0.9993634 ] (action = 9158)
0033-th step: f(x) = [-0.19615996 -0.94679724] (action = 6924)
Pareto set updated.
the number of Pareto frontiers = 9
0034-th step: f(x) = [-0.99950572 -0.99575833] (action = 81)
0035-th step: f(x) = [-0.79425482 -0.43043543] (action = 4491)
Pareto set updated.
the number of Pareto frontiers = 10
0036-th step: f(x) = [-0.94523115 -0.58027949] (action = 2782)
0037-th step: f(x) = [-0.9930067 -0.93998722] (action = 6981)
0038-th step: f(x) = [-0.99560251 -0.10351929] (action = 2653)
0039-th step: f(x) = [-0.97602834 -0.99134066] (action = 2409)
0040-th step: f(x) = [-0.05507465 -0.95545924] (action = 6527)
Pareto set updated.
the number of Pareto frontiers = 10
0041-th step: f(x) = [-0.54838423 -0.96652902] (action = 7723)
0042-th step: f(x) = [-0.97560369 -0.9981919 ] (action = 9623)
0043-th step: f(x) = [-0.8884927 -0.23687603] (action = 4183)
Pareto set updated.
the number of Pareto frontiers = 10
0044-th step: f(x) = [-0.97621647 -0.0052244 ] (action = 3367)
Pareto set updated.
the number of Pareto frontiers = 9
0045-th step: f(x) = [-0.5974142 -0.99778865] (action = 9046)
0046-th step: f(x) = [-0.69847999 -0.99969655] (action = 9661)
0047-th step: f(x) = [-0.99999523 -0.89945685] (action = 1112)
0048-th step: f(x) = [-0.99988103 -0.67274056] (action = 2230)
0049-th step: f(x) = [-0.67089612 -0.99973586] (action = 7063)
0050-th step: f(x) = [-0.44743129 -0.8578429 ] (action = 6412)
Pareto set updated.
the number of Pareto frontiers = 10
目的関数の評価値(の array) とそのときの action ID が表示されます。
また、パレート解集合 (Pareto set) が更新されたときにメッセージを表示します。
Pareto set が更新された際に中身を表示したい場合は、
disp_pareto_set=True と指定します。Pareto set は1つ目の目的関数値の昇順でソートされています。
[11]:
policy = physbo.search.discrete_multi.Policy(test_X=test_X, num_objectives=2)
policy.set_seed(0)
res_random = policy.random_search(max_num_probes=num_random_search + num_bayes_search, simulator=simu, disp_pareto_set=True)
0001-th step: f(x) = [-0.99973003 -0.62334035] (action = 836)
Pareto set updated.
current Pareto set = [[-0.99973003 -0.62334035]] (steps = [1])
0002-th step: f(x) = [-0.99789981 -0.99866427] (action = 9404)
Pareto set updated.
current Pareto set = [[-0.99973003 -0.62334035]
[-0.99789981 -0.99866427]] (steps = [1 2])
0003-th step: f(x) = [-0.99090897 -0.46609239] (action = 4664)
Pareto set updated.
current Pareto set = [[-0.99090897 -0.46609239]] (steps = [3])
0004-th step: f(x) = [-0.92633083 -0.29208351] (action = 4780)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]] (steps = [4])
0005-th step: f(x) = [-0.67969326 -0.99981691] (action = 9566)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.67969326 -0.99981691]] (steps = [4 5])
0006-th step: f(x) = [-0.45601619 -0.99848443] (action = 8852)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.45601619 -0.99848443]] (steps = [4 6])
0007-th step: f(x) = [-0.92670204 -0.71508873] (action = 6088)
0008-th step: f(x) = [-0.58233995 -0.99952931] (action = 7060)
0009-th step: f(x) = [-0.99848229 -0.96780195] (action = 473)
0010-th step: f(x) = [-0.80479332 -0.99994946] (action = 7573)
0011-th step: f(x) = [-0.99700074 -0.99847873] (action = 1406)
0012-th step: f(x) = [-0.9992592 -0.93891121] (action = 6061)
0013-th step: f(x) = [-0.19963873 -0.93357674] (action = 5722)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.19963873 -0.93357674]] (steps = [ 4 13])
0014-th step: f(x) = [-0.98046765 -0.99294428] (action = 2309)
0015-th step: f(x) = [-0.99602549 -0.98620358] (action = 7989)
0016-th step: f(x) = [-0.99957128 -0.9973798 ] (action = 8484)
0017-th step: f(x) = [-0.52191048 -0.72845916] (action = 5405)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.52191048 -0.72845916]
[-0.19963873 -0.93357674]] (steps = [ 4 17 13])
0018-th step: f(x) = [-0.99916441 -0.40869572] (action = 2742)
0019-th step: f(x) = [-0.99480122 -0.7565076 ] (action = 1266)
0020-th step: f(x) = [-0.63329589 -0.63329589] (action = 5200)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.63329589 -0.63329589]
[-0.52191048 -0.72845916]
[-0.19963873 -0.93357674]] (steps = [ 4 20 17 13])
0021-th step: f(x) = [-0.95437918 -0.80142908] (action = 2487)
0022-th step: f(x) = [-0.99899466 -0.96646532] (action = 269)
0023-th step: f(x) = [-0.19473522 -0.99445365] (action = 8044)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.63329589 -0.63329589]
[-0.52191048 -0.72845916]
[-0.19963873 -0.93357674]
[-0.19473522 -0.99445365]] (steps = [ 4 20 17 13 23])
0024-th step: f(x) = [-0.99969529 -0.52395588] (action = 2135)
0025-th step: f(x) = [-0.59106078 -0.79258035] (action = 6206)
0026-th step: f(x) = [-0.78231041 -0.99997141] (action = 9079)
0027-th step: f(x) = [-0.99955573 -0.99930147] (action = 396)
0028-th step: f(x) = [-0.99825097 -0.99875436] (action = 9403)
0029-th step: f(x) = [-0.65387719 -0.99938669] (action = 9456)
0030-th step: f(x) = [-0.26132949 -0.87913689] (action = 5816)
Pareto set updated.
current Pareto set = [[-0.92633083 -0.29208351]
[-0.63329589 -0.63329589]
[-0.52191048 -0.72845916]
[-0.26132949 -0.87913689]
[-0.19963873 -0.93357674]
[-0.19473522 -0.99445365]] (steps = [ 4 20 17 30 13 23])
0031-th step: f(x) = [-0.74523952 -0.7724917 ] (action = 6301)
0032-th step: f(x) = [-0.33644513 -0.97772424] (action = 5530)
0033-th step: f(x) = [-0.99777557 -0.93373833] (action = 6470)
0034-th step: f(x) = [-0.9725524 -0.6296562] (action = 2277)
0035-th step: f(x) = [-0.99924826 -0.84674905] (action = 353)
0036-th step: f(x) = [-0.99743705 -0.17848438] (action = 2549)
Pareto set updated.
current Pareto set = [[-0.99743705 -0.17848438]
[-0.92633083 -0.29208351]
[-0.63329589 -0.63329589]
[-0.52191048 -0.72845916]
[-0.26132949 -0.87913689]
[-0.19963873 -0.93357674]
[-0.19473522 -0.99445365]] (steps = [36 4 20 17 30 13 23])
0037-th step: f(x) = [-0.99994916 -0.78013199] (action = 326)
0038-th step: f(x) = [-0.99887174 -0.86936149] (action = 5458)
0039-th step: f(x) = [-0.99362699 -0.99769786] (action = 9309)
0040-th step: f(x) = [-0.84889217 -0.99973225] (action = 10156)
0041-th step: f(x) = [-0.95713719 -0.09067194] (action = 4073)
Pareto set updated.
current Pareto set = [[-0.95713719 -0.09067194]
[-0.92633083 -0.29208351]
[-0.63329589 -0.63329589]
[-0.52191048 -0.72845916]
[-0.26132949 -0.87913689]
[-0.19963873 -0.93357674]
[-0.19473522 -0.99445365]] (steps = [41 4 20 17 30 13 23])
0042-th step: f(x) = [-0.17190645 -0.91382463] (action = 6120)
Pareto set updated.
current Pareto set = [[-0.95713719 -0.09067194]
[-0.92633083 -0.29208351]
[-0.63329589 -0.63329589]
[-0.52191048 -0.72845916]
[-0.26132949 -0.87913689]
[-0.17190645 -0.91382463]] (steps = [41 4 20 17 30 42])
0043-th step: f(x) = [-0.98406208 -0.94467587] (action = 7289)
0044-th step: f(x) = [-0.93974303 -0.92444262] (action = 2698)
0045-th step: f(x) = [-0.98414342 -0.84762781] (action = 1780)
0046-th step: f(x) = [-0.99699986 -0.98166426] (action = 7684)
0047-th step: f(x) = [-0.3868143 -0.99930896] (action = 8060)
0048-th step: f(x) = [-0.97964062 -0.65554315] (action = 2075)
0049-th step: f(x) = [-0.99907307 -0.73466786] (action = 4450)
0050-th step: f(x) = [-0.99720717 -0.85352507] (action = 865)
結果の確認
#### 評価値の履歴
[12]:
res_random.fx[0:res_random.num_runs]
[12]:
array([[-0.99973003, -0.62334035],
[-0.99789981, -0.99866427],
[-0.99090897, -0.46609239],
[-0.92633083, -0.29208351],
[-0.67969326, -0.99981691],
[-0.45601619, -0.99848443],
[-0.92670204, -0.71508873],
[-0.58233995, -0.99952931],
[-0.99848229, -0.96780195],
[-0.80479332, -0.99994946],
[-0.99700074, -0.99847873],
[-0.9992592 , -0.93891121],
[-0.19963873, -0.93357674],
[-0.98046765, -0.99294428],
[-0.99602549, -0.98620358],
[-0.99957128, -0.9973798 ],
[-0.52191048, -0.72845916],
[-0.99916441, -0.40869572],
[-0.99480122, -0.7565076 ],
[-0.63329589, -0.63329589],
[-0.95437918, -0.80142908],
[-0.99899466, -0.96646532],
[-0.19473522, -0.99445365],
[-0.99969529, -0.52395588],
[-0.59106078, -0.79258035],
[-0.78231041, -0.99997141],
[-0.99955573, -0.99930147],
[-0.99825097, -0.99875436],
[-0.65387719, -0.99938669],
[-0.26132949, -0.87913689],
[-0.74523952, -0.7724917 ],
[-0.33644513, -0.97772424],
[-0.99777557, -0.93373833],
[-0.9725524 , -0.6296562 ],
[-0.99924826, -0.84674905],
[-0.99743705, -0.17848438],
[-0.99994916, -0.78013199],
[-0.99887174, -0.86936149],
[-0.99362699, -0.99769786],
[-0.84889217, -0.99973225],
[-0.95713719, -0.09067194],
[-0.17190645, -0.91382463],
[-0.98406208, -0.94467587],
[-0.93974303, -0.92444262],
[-0.98414342, -0.84762781],
[-0.99699986, -0.98166426],
[-0.3868143 , -0.99930896],
[-0.97964062, -0.65554315],
[-0.99907307, -0.73466786],
[-0.99720717, -0.85352507]])
パレート解の取得
[13]:
front, front_num = res_random.export_pareto_front()
front, front_num
[13]:
(array([[-0.95713719, -0.09067194],
[-0.92633083, -0.29208351],
[-0.63329589, -0.63329589],
[-0.52191048, -0.72845916],
[-0.26132949, -0.87913689],
[-0.17190645, -0.91382463]]),
array([40, 3, 19, 16, 29, 41]))
解（評価値）のプロット
これ以降、図示する空間が \(y = (y_1, y_2)\) であり \(x = (x_1, x_2)\) ではないことにあらためて注意してください。
赤のプロットがパレート解です。
[14]:
fig, ax = plt.subplots(figsize=(7, 7))
physbo.search.utility.plot_pareto_front_all(res_random, ax=ax)
ax.set_xlim(-1, 0)
ax.set_ylim(-1, 0)
ax.set_title("random search")
[14]:
Text(0.5, 1.0, 'random search')
劣解領域 (dominated region) の体積を計算
パレート解ではない解、すなわち、「自らよりも優れた解 \(y'\) が存在する解 \(y\) 」 を劣解と呼びます (\(\exists y' y \prec y'\))。
解空間（の部分空間）のうち、劣解の占める空間である劣解領域の体積は、多目的最適化の結果を示す指標のひとつです。
この値が大きいほど、より良いパレート解が多く求まっていることになります。
res_random.pareto.volume_in_dominance(ref_min, ref_max) は、 ref_min, ref_max で指定された矩形(hyper-rectangle) 中の劣解領域体積を計算します。
[15]:
VID_random = res_random.pareto.volume_in_dominance([-1,-1],[0,0])
VID_random
[15]:
np.float64(0.2376881844865093)
ベイズ最適化
多目的の場合の bayes_search では、score には以下のいずれかを指定します。
HVPI (HyperVolume-based Probability of Improvement)
EHVI (Expected Hyper-Volume Improvement)
TS (Thompson Sampling)
以下、score を変えてそれぞれ 50回 (ランダムサーチ10回 + ベイズ最適化 40回) 評価を行います。
HVPI (HyperVolume-based Probability of Improvement)
多次元の目的関数空間における非劣解領域 (non-dominated region) の改善確率を score として求めます。
参考文献
Couckuyt, Ivo, Dirk Deschrijver, and Tom Dhaene. 「Fast calculation of multiobjective probability of improvement and expected improvement criteria for Pareto optimization.」 Journal of Global Optimization 60.3 (2014): 575-594.
[16]:
policy = physbo.search.discrete_multi.Policy(test_X=test_X, num_objectives=2)
policy.set_seed(seed)
policy.random_search(max_num_probes=num_random_search, simulator=simu, is_disp=False)
time_start = time.time()
res_HVPI = policy.bayes_search(
max_num_probes=num_bayes_search,
simulator=simu,
score="HVPI",
interval=10,
is_disp=False,
)
time_HVPI = time.time() - time_start
print(f"HVPI の計算時間: {time_HVPI} 秒")
HVPI の計算時間: 4.140035152435303 秒
パレート解のプロット
physbo.search.utility.plot_pareto_front_all に steps_end と steps_begin を指定することで、プロットするステップ数を指定できます。 style_common を指定することで、プロットの見た目を変更できます。 また、 style_dominated と style_pareto_front で、劣解領域とパレート解のプロットの見た目をそれぞれ変更できます。 ここでは最初のランダムサンプリングの結果を marker="+" でプロットし、その後のベイズ最適化の結果を marker="o" でプロットしています。
[17]:
fig, ax = plt.subplots(figsize=(7, 7))
physbo.search.utility.plot_pareto_front_all(
res_HVPI,
ax=ax,
steps_end=num_random_search,
style_common={"marker": "+"},
)
physbo.search.utility.plot_pareto_front_all(
res_HVPI,
ax=ax,
steps_begin=num_random_search,
style_common={"marker": "o"},
)
ax.set_xlim(-1, 0)
ax.set_ylim(-1, 0)
ax.set_title("HVPI")
[17]:
Text(0.5, 1.0, 'HVPI')
劣解領域体積
[18]:
VID_HVPI = res_HVPI.pareto.volume_in_dominance([-1,-1],[0,0])
VID_HVPI
[18]:
np.float64(0.3285497837573861)
ランダムサンプリングを50回おこなった結果と比べて、多くのパレート解が見つかっています。
EHVI (Expected Hyper-Volume Improvement)
多次元の目的関数空間における非劣解領域 (non-dominated region) の改善期待値を score として求めます。
参考文献
Couckuyt, Ivo, Dirk Deschrijver, and Tom Dhaene. 「Fast calculation of multiobjective probability of improvement and expected improvement criteria for Pareto optimization.」 Journal of Global Optimization 60.3 (2014): 575-594.
[19]:
policy = physbo.search.discrete_multi.Policy(test_X=test_X, num_objectives=2)
policy.set_seed(seed)
policy.random_search(max_num_probes=num_random_search, simulator=simu, is_disp=False)
time_start = time.time()
res_EHVI = policy.bayes_search(
max_num_probes=num_bayes_search,
simulator=simu,
score="EHVI",
interval=10,
is_disp=False,
)
time_EHVI = time.time() - time_start
print(f"EHVI の計算時間: {time_EHVI} 秒")
EHVI の計算時間: 34.10754203796387 秒
パレート解のプロット
[20]:
fig, ax = plt.subplots(figsize=(7, 7))
physbo.search.utility.plot_pareto_front_all(
res_EHVI,
ax=ax,
steps_end=num_random_search,
style_common={"marker": "+"},
)
physbo.search.utility.plot_pareto_front_all(
res_EHVI,
ax=ax,
steps_begin=num_random_search,
style_common={"marker": "o"},
)
ax.set_xlim(-1, 0)
ax.set_ylim(-1, 0)
ax.set_title("EHVI")
[20]:
Text(0.5, 1.0, 'EHVI')
劣解領域体積
[21]:
VID_EHVI = res_EHVI.pareto.volume_in_dominance([-1,-1],[0,0])
VID_EHVI
[21]:
np.float64(0.32220853509132785)
TS (Thompson Sampling)
単目的の場合の Thompson Sampling では、各候補(test_X)について、目的関数の事後分布からサンプリングを行い、値が最大となる候補を次の探索点として推薦します。
多目的の場合は、サンプリングした値についてパレートルールの上で最大となる候補、つまりパレート最適な候補の中からランダムに1つ選択して次の探索点とします。
参考文献
Yahyaa, Saba Q., and Bernard Manderick. 「Thompson sampling for multi-objective multi-armed bandits problem.」 Proc. Eur. Symp. Artif. Neural Netw., Comput. Intell. Mach. Learn.. 2015.
[22]:
policy = physbo.search.discrete_multi.Policy(test_X=test_X, num_objectives=2)
policy.set_seed(seed)
policy.random_search(max_num_probes=num_random_search, simulator=simu, is_disp=False)
time_start = time.time()
res_TS = policy.bayes_search(
max_num_probes=num_bayes_search,
simulator=simu,
score="TS",
interval=10,
num_rand_basis=500,
is_disp=False,
)
time_TS = time.time() - time_start
print(f"TS の計算時間: {time_TS} 秒")
TS の計算時間: 13.531005859375 秒
パレート解のプロット
[23]:
fig, ax = plt.subplots(figsize=(7, 7))
physbo.search.utility.plot_pareto_front_all(
res_TS,
ax=ax,
steps_end=num_random_search,
style_common={"marker": "+"},
)
physbo.search.utility.plot_pareto_front_all(
res_TS,
ax=ax,
steps_begin=num_random_search,
style_common={"marker": "o"},
)
ax.set_xlim(-1, 0)
ax.set_ylim(-1, 0)
ax.set_title("TS")
[23]:
Text(0.5, 1.0, 'TS')
劣解領域体積
[24]:
VID_TS = res_TS.pareto.volume_in_dominance([-1,-1],[0,0])
VID_TS
[24]:
np.float64(0.3061697389817548)
比較
これまでの結果を比較してみましょう。
[25]:
import pandas as pd
pd.DataFrame({
"アルゴリズム": ["HVPI", "EHVI", "TS"],
"計算時間": [time_HVPI, time_EHVI, time_TS],
"劣解領域体積": [VID_HVPI, VID_EHVI, VID_TS]
})
[25]:
| アルゴリズム | 計算時間 | 劣解領域体積 | |
|---|---|---|---|
| 0 | HVPI | 4.140035 | 0.328550 |
| 1 | EHVI | 34.107542 | 0.322209 |
| 2 | TS | 13.531006 | 0.306170 |
付録：全探索
random_search で max_num_probes に全データ数 (N = test_X.shape[0]) を渡すと手軽に全探索できます。全データの評価には時間がかかるため、あらかじめデータ数を減らしておきます。
[26]:
test_X_sparse = physbo.search.utility.make_grid(min_X=min_X, max_X=max_X, num_X=21)
simu_sparse = physbo.search.utility.Simulator(test_X_sparse, vlmop2_minus)
policy = physbo.search.discrete_multi.Policy(test_X=test_X_sparse, num_objectives=2)
policy.set_seed(0)
N = test_X_sparse.shape[0]
res_all = policy.random_search(max_num_probes=N, simulator=simu_sparse, is_disp=False)
パレート解のプロット
[27]:
fig, ax = plt.subplots(figsize=(7, 7))
physbo.search.utility.plot_pareto_front_all(res_all, ax=ax)
[27]:
array([[<Axes: xlabel='Objective 1', ylabel='Objective 2'>]], dtype=object)
