Bayesian Optimization / Gaussian Process

前回、matplotlibのアニメーションについては大体使えるようになったので、「ベイズ最適化・ガウス過程」を実験してみました。いろいろと専用のライブラリがあるのですが（GPy、GPFlow、Edward、Pyroなど）、今回はscikit-learn Gaussian Processesを使っています。

ベイズ最適化はハイパーパラメータチューニングに使われるようで、ブラックボックス関数における最大値（あるいは最小値）を全領域をスキャンしなくても探し出すことができる方法のようです。

今回の場合は、グラフ内の赤破線の最大値がどこにあるか探索しています。赤破線を見ると、左と右にピークが二箇所ありますが、微妙に右側のほうが高いという設定にしてあります。それを間違えないように探し出せるかという実験です。

右側の赤丸が正解（最大値）の地点です。この地点は、x座標-2〜2までの範囲を0.1刻みで200回スキャンし、その結果から最大値を求めています。

しかし、ベイズ最適化の場合は、200回スキャンしなくても約10回程度で、どこが最大値なのかを効率よく探してくれるようです。

今回の手順として：
・カーネル：Matern、獲得関数：UCBを選択

・xの範囲（-2〜2）からランダムにまず2点（青丸：start印）選んで初期観測値とする

・その観測値から、全域のμ（青破線）とσ（青塗）を求める

・μとσを獲得関数（UCB）に代入

・獲得関数の結果（緑実線）から最大値の箇所を見つける

・獲得関数の最大値から、それに相当するxの位置を見つける

・xの位置から、それに相当するyの位置を求める

・ここで求めたxとyの座標を次回の探索箇所にする

・次回探索箇所を観測値リストに加える

・あとは繰り返し

グラフ下部（緑実線）にあるのが獲得関数UCBの値です。緑実線の一番高いところが、次の探索地のx座標となります。

今回使用したUCBの中身は単純で、

UCB最大値 = mu + kappa * sigma

となっています。係数kappa=1.0の場合は、μ値（青破線）にσ値（青塗）を足した値がそのままUCB最大値となるので、グラフ内の青塗を含めた一番高い箇所が次の探索箇所となります。μ値が高いところほど選ばれがちですが、kappaを大きくするとσ値も大きくなるので、μ値が低い箇所でも探索される可能性があがります。

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
import matplotlib.animation as animation
from IPython.display import HTML

np.random.seed(111)

def f(x):
    return -np.sin(x*4.5)*x*0.3+0.02*x                  # 真の関数（赤破線表示）

x = np.linspace(-2,2,200)                               # x値範囲
y = f(x)                                                # 真のy値

fig = plt.figure(figsize=(12,6))                        # 画面サイズ設定
ax = plt.subplot()                                      # サブプロットを用意
ax.axis([x.min(),x.max(),y.min()*2,y.max()*2])          # グラフ範囲設定（x,yの値に依存）
ax.set_aspect(1.0)                                      # アスペクト比設定

ax.plot(x,y,'r--',alpha=0.8)                            # 真の関数曲線をプロット（赤破線）
ax.plot(x[np.argmax(y)],y.max(), 'ro',alpha=0.8)        # y正解値；x全域スキャンによるy最大値（赤丸）

init_pt = 2                                             # 初期観測値個数
xi = np.random.uniform(x.min(),x.max(),init_pt)         # 初期観測値：xの範囲から2箇所ランダム選択
yi = f(xi)                                              # 初期観測値：yの値
ax.plot(xi,yi,'bo')                                     # 初期観測値プロット（青丸）
for i in range(len(xi)):
    ax.text(xi[i],yi[i],'  start'+str(i),color='b')     # 初期観測値：'start'文字表示

def UCB(mu_x, sigma_x, kappa=1.0):                      # 今回使用する獲得関数：upper_confidence_bound
    return mu_x + kappa * sigma_x                       # kappa増大-->σ増大-->探索領域拡大

kernel = Matern(length_scale_bounds="fixed")            # 今回使用するカーネル：Matern
gp = GaussianProcessRegressor(kernel=kernel,alpha=0.01) # scikit-learnガウス過程回帰オブジェクト

ims = []                                                # アニメ用リスト
itr = 15                                                # 探索ループ回数

for i in range(itr):
    gp.fit(xi[:,None],f(xi))                                 # 観測値フィッティング（xiをリシェイプして代入）
    mu_x, sigma_x = gp.predict(x[:, None], return_std=True)  # μ(x)とσ(x)を獲得（事後分布）

    a_x = UCB(mu_x, sigma_x, kappa=2.0)                      # UCBでa(x)最大値を得る
    index = np.argmax(a_x)                                   # a(x)最大値のindexを得る
    x_next = x[index]                                        # a(x)最大値のindexから次のx探索値を決定
    y_next = f(x_next)                                       # y探索値決定
    
    P1 = ax.plot(x, mu_x, 'b--')                                               # μ値（青破線）
    P2 = ax.fill_between(x, mu_x-sigma_x, mu_x+sigma_x, color='b', alpha=0.2)  # σの範囲（青塗）
    P3 = ax.plot(xi,yi,'bx',markersize=8)                                      # 探索済み箇所（青クロス）
    P4 = ax.plot(x_next, y_next, 'mo')                                         # 現在探索箇所（マジェンタ丸）
    P5 = ax.plot([x_next,x_next],[y.min()*2, y_next] ,'m-', alpha=0.6)         # 現在探索箇所（マジェンタ縦線）
    P6 = ax.plot(x_next,(a_x[index]-a_x.min())/2+y.min()*2, 'go',alpha=0.8)    # a(x)最大値：σ（青丸）
    P7 = ax.plot(x, (a_x-a_x.min())/2+y.min()*2, 'g-')                         # 獲得関数値：a(x)
    P8 = ax.plot([xi,xi],[[y.min()*2]*len(yi),yi], 'b-', alpha=0.2)            # 探索済み箇所（青縦線）
       
    xi = np.r_[xi,x_next]                                  # 現在探索箇所をリストに加える（r_: concatenateの代用）
    yi = np.r_[yi,y_next]
    
    ims.append(P1 + [P2] + P3 + P4 + P5 + P6 + P7 + P8)    # アニメ用リストへ各プロットを入れる
    
print('x :',x[y.argmax()], '   y :',y.max())               # 正解値出力（y値）
print('xi:',xi[yi.argmax()], '   yi:',yi.max())            # 最終最大値出力（yi値）

ani = animation.ArtistAnimation(fig, ims, interval=1000)   # アニメ関数
#ani.save("gp_sim.gif", writer = "imagemagick")            # GIFアニメ保存
plt.close()
HTML(ani.to_html5_video())                                 # HTML5 Video 出力（mp4ダウンロード可）
#HTML(ani.to_jshtml())                                     # JavascriptHTML 出力

GaussianProcessRegressor()にはカーネルを指定する箇所があり、今回は、scikit-learnのMaternを使っています。ちなみにカーネルを指定せずデフォルトのままだと、以下のような感じになります。

σ（青塗）の領域が結構変わります。
scikit-learnのこのページにカーネルの違いについて書いてあります。

ハイパーパラメータチューニングと言っても、カーネルの選び方でもかなり違いがでてきます。
また、獲得関数もUCB以外にいくつかあるので、それも選ばなければいけないという感じです。


参考にしたサイト：
Scikit-learn Gaussian Processes
llSourcell/hyperparameter_optimization_strategies
ベイズ最適化入門
機械学習のハイパーパラメータ探索 : ベイズ最適化の活用

関連：
GPyOpt: Digital Recognizer(MNIST) CNN Keras ハイパーパラメータ最適化