はじめよう実験計画

実験を早く終わらせるための技術

ベイズ最適化と実験計画法

f:id:Sturgeon:20210619121509p:plain

はじめに

最近、ベイズ最適化による実験の短縮が流行っているようで。ベイズ最適化について調べてみると、関数の最適化であったり、機械学習におけるハイパーパラメータの調整に使われるケースが多いイメージがあります。

そこで、本記事では「ベイズ最適化を実験に使う」ことに重点を置いて話します。

 

実験におけるベイズ最適化の利用

ベイズ最適化を実験で使用する場合、図1に示すフローが一般的です。

 

f:id:Sturgeon:20210619123906p:plain

図1. ベイズ最適化を用いた実験のフロー

 

ベイズ最適化では、「ガウス過程回帰」という回帰の親玉のような手法を用いて、実験データをモデル化(入力データに対して予測値を返す関数を作る)します。ガウス過程回帰では、線形回帰では表現できないような複雑な応答を表現することが出来ます。

ガウス過程回帰ではある実験条件における予測値µとデータの予測分散σを計算することができ、このµとσを考慮した「獲得関数」という評価値に基づいて、最適値の提案を行います。図2にこのイメージを示しました。

 

f:id:Sturgeon:20210626104113p:plain

図2. 1変数のベイズ最適化における獲得関数と実験点の提案

ガウス過程回帰モデルによる予測値µと予測分散σを用いた獲得関数を計算し、獲得関数の最大点を次の実験点とする。

 

獲得関数には、UCB、EI、PIという代表的な3種類があります。最も単純な獲得関数はUCBであり、

UCB = µ + k×σ  (kは実験回数に応じて変化させる係数)

という形をしています。したがって、最大化問題であれば、実験データから予測される値µが大きく、かつ予測分散σが大きな点で獲得関数が大きくなります。

獲得関数はよく「活用」と「探索」という言葉で説明されます。すなわち、予測値µが大きくなるようにするのが「活用」であり、分からなさσが大きな領域を調べるのが「探索」です。このように、ベイズ最適化では、予測値だけなく、分からなさも考慮した実験条件の提案が行われます。

ベイズ最適化の実装自体は難しくないですが、ガウス過程回帰の理論は抑えたうえでやるべきだと思います。ガウス過程回帰は結構奥深くて、私の説明では間違いがあるかもしれませんので、ぜひ参考書をお読みください。『ガウス過程回帰と機械学習』という参考書がおすすめです。

 

 

デモ

よくあるベイズ最適化のデモは、図3のような1変数のものです。

 

f:id:Sturgeon:20210619111703p:plain

図3. ベイズ最適化のデモ(1変数)

 

図3でやっていることを説明します。まず、1回目の左図は初期の実験データと真の関数(実際は分からない)を示しています。右図では実験データに基づいてガウス過程回帰モデルを作成しています(黄色の線)、灰色の帯は予測分散σを表しており、この灰色の大きさσと黄色の線による予測値µを

獲得関数 = µ+3σ

で計算したものが獲得関数です(赤色の線)。したがって、獲得関数の最大点が予測値と予測分散を考慮した、次の実験点となります。

図3の3回目になると、この獲得関数が威力を発揮してきます。実験されていない右側の領域では、モデル(黄色の線)による予測値があまり大きくありません。しかし、灰色の帯の大きさ(分からなさ)が大きいため、獲得関数(赤色)が大きくなるので、実際の最大点に近い条件が提案されています(図3右下)。

このように、ベイズ最適化では局所的な最適値にとどまらず、未知の点を「探索」する仕組みがあり、自動的な実験条件の最適化ができます。

ただ、これだとイマイチありがたみが無い気がするんですね。というのも、1変数なら最初から数点実験して傾向取ればいいじゃん!ってなります。

そこで、本記事では2変数のベイズ最適化のデモもやってみました。

 

f:id:Sturgeon:20210619112252p:plain

図4. ベイズ最適化のデモ(2変数)

 

こっちの方が実験っぽいかなと思います。図4のすべてのグラフにおいて、青色の曲面が真の関数で(実際は分からない)、黄色の曲面がガウス過程回帰モデルです。

このモデルによる予測値と予測分散を考慮して、先述の獲得関数を考慮して、次の実験点(赤丸)が提案されています。

図4では、5回目でモデル(黄色)の形が真の関数(青色)にぐっと近づいて、提案点も最適値に迫っているのがわかります。

 

実装

図4のデモをするためのPythonコードをこちらに書きました。ガウス過程回帰はSckiit-learnのページを参考にしました。

 

 

import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel, DotProduct, ConstantKernel
from mpl_toolkits.mplot3d import Axes3D

#獲得関数
def aq(mu, sigma, N):
    a = mu+3*sigma # a = mu + k * sigmaでkを実験回数で変更してもよい
    i = np.argmax(a)
    return a, i

#真の関数(実際は分からない)
def g(X,Y):
    u = 2*np.exp(-7*((X+0.5)**2+(Y+0.5)**2-(X+0.5)*(Y+0.5)))
    v = 3*np.exp(-5*((X-0.3)**2+(Y-0.5)**2-(X-0.3)*(Y-0.5)))
    w = np.exp(-3*((X-0.6)**2+(Y+0.6)**2-(X-0.6)*(Y+0.6)))
    return u+v+w

#初期の実験点
D_measured = np.array([[-0.7, -0.7],
                    [-0.7, 0.7],
                    [0.7, -0.7],
                    [0.7, 0.7],
                    [0,0]])
NP = np.array([[0,0]])
NP = np.delete(NP, 0, 0)
print("Reset NP:", NP)
print(D_measured)


print(NP.shape[0], "Newly measured points: ")
if NP.shape[0] > 0:
    print("The latest data measured at x=", NP[-1])


"""以下コードの繰返し"""
#  観測  
D_train = np.vstack((D_measured, NP))
np.random.seed(3)
Z_train = np.zeros(len(D_train))
for i in range(0,len(D_train)):
    Z_train[i] = g(D_train[i][0], D_train[i][1])+np.random.standard_normal(1)*0.1

# ガウス過程回帰モデルの構築
kernel = ConstantKernel()*RBF(length_scale_bounds=[0.01, 0.5])+WhiteKernel(noise_level_bounds=[1e-5, 0.5])
GP = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=20)
GP.fit(D_train, Z_train)
print("Kernel: ", GP.kernel_)

# ガウス過程回帰モデルによる予測値の計算
X1 = np.arange(-1, 1, 0.05)
X2 = np.arange(-1, 1, 0.05)
# Meshgrid for 3D plot
X1X1, X2X2 = np.meshgrid(X1, X2, indexing='ij')
#X1X2 is argument for GP
X1X2 = np.c_[X1X1.flatten(), X2X2.flatten()]
G, SD = GP.predict(X1X2, return_std=True)
G_mesh = G.reshape((len(X1), len(X2)))#mesh for 3D plot

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X1X1, X2X2, G_mesh, rstride=1, cstride=1, alpha=0.5, color='y')#prediction
for i in range(0, len(D_train)):
    ax.scatter(D_train[i][0], D_train[i][1], Z_train[i], color='b')

# 真の関数(実際は分からない)の表示
Z_true = np.zeros(len(X1X2))
for i in range(0,len(X1X2)):
    Z_true[i] = g(X1X2[i][0], X1X2[i][1])
Z_true_mesh = Z_true.reshape((len(X1), len(X2)))
ax.plot_surface(X1X1, X2X2, Z_true_mesh, rstride=1, cstride=1, alpha=0.3, color='b')#true function

# 次の実験点の提案
# 獲得関数の計算
A, index = aq(G, SD, len(X1X2))
A_mesh = A.reshape((len(X1), len(X2)))
indices = np.argwhere(A_mesh==A[index])#the location of maximum Aq
NP_x1 = X1[indices[0][0]]
NP_x2 = X2[indices[0][1]]
ax.scatter(NP_x1, NP_x2, G[index], color='r', s=50)
plt.show()

print("Suggested point NP is (%f,%f)" %(NP_x1, NP_x2))
NP = np.vstack((NP, [NP_x1, NP_x2]) )#add the suggested point to NP
"""以上コードの繰返し"""

 

上のコードには反映していないのですが、過学習(今あるデータの当てはまりが良すぎて、将来の予測に使えないこと)を防ぐために、クロスバリデーションによってガウス過程回帰モデルのカーネルを定める必要があります。

実験にベイズ最適化を用いる場合、実験点は数点~数十点でしょうから、テストデータは数点が良いでしょう。

 

ベイズ最適化による実験の短縮についての疑問

前述の図4に示したように、ガウス過程回帰は非常に柔軟で、複雑な応答を滑らかな曲線で表現することが出来ます。ベイズ最適化では、そのガウス過程回帰モデルを利用して、効率的な実験点を自動的に提案してくれます。

そのため、「実験の自動化」という点では間違いなく利用価値があります。しかし、ベイズ最適化それ自身には実験を短縮する効果はないのでは?と個人的には思います。

実際にガウス過程をやってみればわかるのですが、ガウス過程回帰は「観測した点の滑らかに近似するだけ」の手法です。別の言い方をすれば、

観測していない点に関しては何もわからない。

事実、未知の領域について何も分からないから、獲得関数に予測分散を反映して実験していない点を提案するようにしているわけです(図5)。

 

f:id:Sturgeon:20210626105953p:plain

図5, ベイズ最適化では、予測分散を反映した獲得関数が大きくなる未知の領域を探索する。なぜなら実験していない領域については何もわからないから。

 

これは人間と同じどころか、人間のように現象の物理を考えたりしているわけではないので、行き当たりばったりのモグラたたき探索になる危険性があります。

 

ベイズ最適化によって実験回数が大幅に短縮する場合、それは線形モデルのような簡単なモデルでも表せる現象である可能性が高いです。この場合、数式としての解釈ができる線形モデルの方が適していると思います。

複雑なモデルでしか表せない現象であれば、先ほど言ったように、ベイズ最適化は自動的な実験条件の探索として利用する価値があります。しかし複雑な現象では「実験しなければわからない」ため、実験回数の短縮はできないでしょう。

本当の意味で、実験回数を削減するためには、実験していない点の予測(外挿)ができる必要があります。外挿は、様々な物理を考慮して総合的に判断できる「人間」にしか今のところできないんじゃないかな?と思います。それが出来るスーパー人工知能が現れたら恐ろしいですが...

 

ベイズ最適化と要因の数と実験計画法の利用

多くの変数(要因)を調整して、最適な実験条件を見つけたいとき、すべての変数を用いてベイズ最適化を行うことは効率的とは言えません。なぜなら、変数が増えれば増えるほど、実験していない領域が多く(変数空間の外側の領域は変数の数のべき乗で増える)、探索を重視する傾向のあるベイズ最適化の初期段階では、実験領域の端ばかり提案するからです。数個の要因を扱う場合は、それでもいいのですが、要因が10, 20となると、実験がいつまでたっても終わりません。

そのため、多くの要因を用いた最適化を行う場合、スクリーニングの目的で実験計画法を用いて、効果の小さい要因を除外し、効果の大きな要因だけでベイズ最適化を行うべきです。

他の分野についてはよく分かりませんが、少なくとも工学の分野に関して言えば、オッカムの剃刀(ある事柄を説明するためには、必要以上に多くを仮定するべきでない)が正しいと思います。

まとめると、

  1. 実験計画法により多くの要因(例えば10~20個)から効果の大きな要因(例えば3~6個)に絞る。
  2. ベイズ最適化の利用

このような、実験計画法+ベイズ最適化という2ステップであれば、多くの要因を扱う実験の回数削減につながるでしょう。スクリーニング目的の実験計画としては、Plackett-Burmman計画や決定的スクリーニング計画をお勧めします。

 

www.doe-get-started.com

 

実験計画法の参考書としては、以下をお勧めしています。

 

Design and Analysis of Experiments, Wiley; 第10版 (2020/6/23)

待望の新版で、比較的新しい「決定的スクリーニング計画」の説明が追加されました!