ベイズ最適化 (その2)：獲得関数 - グラフ機械学習と強化学習について

前回の記事ではベイズ最適化で使用されるガウス過程回帰（Gaussian Process Regression）についてまとめていきました。今回の記事では、ガウス過程を用いたベイズ最適化について行っていきたいと思います。よく用いられるのは、ハイパーパラメータの探索や、実験点の探索などです。利用と探索を考慮した最適化手法になります。

github.com

スターが多かったのでBaysianOptimizationのパッケージを使って試してみたいと思います。ベイズ最適化ではないですが、facebook reserachの最適化ライブラリNervergradも気になっています。

理論について下記を参考にしています。

[1807.02811] A Tutorial on Bayesian Optimization

以下の記事も大変よくまとまっております。ベイズ最適化の例では、主にこちらのコードを使わせていただきました。

krasserm.github.io

ベイズ最適化

ベイズ最適化は、機械学習ベースの最適化法の一種です。数式で記述すると次の問題を解くことと同一です。

\begin{align} x^* = \arg \max_{x \in A} f(x) \end{align}

すなわち、ある目的関数 $f(x)$ を最大にするようなxを求めていくことになります。

ベイズ最適化を行う前提条件として、以下のことを想定しています。

入力xは $\mathcal{R}^d$ が多すぎない。 $d \leq 20$ が成功しやすい。
パラメータセットAはシンプルである。Aはhyper rectangle $\{x \in \mathcal{R}^d : a_i \leq x_i \leq b_i\}$ もしくはd次元のsimplex $\{x \in \mathcal{R}^d : \sum_i x_i = 1\}$ が良い。
目的関数 $f$ は連続である。ガウス過程回帰が良く用いられる。
$f$ は"評価することが困難"である。評価自体が難しい or 時間・コストがかかるなど。
$f$ はブラックボックスである。凸や線形性などもわからない。
$f$ を評価するとき、f(x)だけあり、1次導関数や2次導関数がない。導関数があるならニュートン法や準ニュートン法を用いられる。
ノイズは各評価に対して独立。
局所解ではなく大域的最適化を見つける。

ブラックボックスな微分不要な関数の最適化ならばベイズ最適化がよく用いられるようになりました。特にニューラルネットワークのハイパーパラメータの調整（Snoek et al, 2012）です。また、AlphaGoでもベイズ最適化が用いられています(Yutian et al., arxiv)。

アルゴリズム概要

f:id:udnp:20190420195143p:plain

アルゴリズムの考え方自体は、非常にわかりやすいです。ガウシアンプロセスなどを用いて値を観察します。すべてのデータを用いて事後確率分布を更新し、このときの分布から獲得関数が最大値となる $x_n$ を求めた後に、 $y_n = f(x_n)$ を求めます。この操作を繰り返していきます。 $f(x)$ は平均 $\mu_n (x)$ や分散 $\sigma^{2}_n (x)$ の分布が得られます。

ガウス過程回帰

ガウス過程の理論は前回の記事

ベイズ最適化：ガウシアンプロセス - 創薬・材料探索のための機械学習

となりますが、簡単に復習すると

\begin{align} f(x_{k+1} | x_{1:k}) \sim \mathcal{N}(x_{k+1}|μ_0 (x_{1:k}), Σ_0 (x_{1:k}, x_{1:k})) \end{align}

となるモデルです。この平均と分散は各点 $x_i, x_j$ におけるカーネル関数（共分散関数）により得られます。カーネルの選択はデータの関係性を考慮し、決定します。

Mean Function とカーネル選択

カーネルは入力空間中で点同士が強い相関を示すという特性があります。近い点はより近く、遠くの点はより遠くに写像することができます。また、カーネルは半正定値であることが必要条件です。よく使用されるカーネルをいかに示します

Gaussian kernel (power exponential)

\begin{align} k(x, x') = \alpha \exp ( -||x- x'||^2) \end{align}

ここで、 $||x- x'||^2 = \sum_{i=1}^{d} \alpha_{i} (x-x')^{2}$ です。αはカーネルのパラメータとなります。

Matern kernel

\begin{align} k(x, x') = \alpha_{0} \frac{2^{1-\nu}}{\Gamma(\nu)} (\sqrt{2\nu}||x-x'||)^{\nu} K_{\nu} (\sqrt{2\nu}||x-x'||) \end{align}

ここで、 $K_{\nu}$ は修正ベッセル関数です。いずれもscikit-learnで使用できます。

良く使用される平均関数は $\mu_{0} = \mu$ ですが、何らかのトレンドを有するなどの事前知識を有しているとき、パラメトリックで表した方が良い場合もあります。

\begin{align} \mu_{0} = \mu + \sum_{i=1}^{p} \beta_{i} \Psi_{i}(x) \end{align}

ここで、 $\Psi$ は基底関数です。よく低次元の多項式が使用されます。多項式の場合では正規方程式を解けばよいので次数が小さければ簡単に求まります。

ハイパーパラメータの選択

平均関数とカーネルにはパラメータがありますが、これをハイパーパラメータ（超パラメータ）と呼びます。これを求めるための方法の1つに最尤推定（maximum likelihood estimate, MLE）があります。

\begin{align} \hat{η} = \arg \max_{η} P(f(x) | η) \end{align}

こちらは前回の記事で行いました。scipyのminimizeを使えば簡単にMLEを行えます。もう一つの方法は、事後分布最大化（maximum a posterior, MAP）推定です。事後分布を最大化させるパラメータを決定する方法です。

\begin{align} \hat{η} = \arg \max_{η} P( η | f(x) ) = \arg \max_{η} P(f(x) | η ) P(η) \end{align}

ベイズの定理を用いているがポイントです。分母の正規化定数（エビデンス）はパラメータに依存しないため省略することができます。MLEはMAPの特別な場合です。点推定か分布推定かの違いです。MCMCなどのサンプリング法で求めます。個人的にはHMCが良いと思います。

\begin{align} P(f(x)=y | f(x_{1:n})) &= \int P(f(x)=y| f(x_{1:n}), η) P(η | f(x_{1:n}) dη \\ & \approx \frac{1}{J} \sum_{j=1}^{J} P(f(x)=y| f(x_{1:n}), η=\hat{η_{j}}) \end{align}

MAP推定は完全なベイズ推論とみなせます。

獲得関数 (Acquisition Functions)

効用関数とも呼ばれています。次の点を探すための関数です。ノイズの評価、パラレルな評価、導関数を用いる方法、"exotic " extensionがあります。最もよく用いられるのが、Expected improvementです。簡単なわりに、精度も良いです。

Upper Confidence Bounds, UCB

まずは、UCBについて記載します。平均値と信頼区間の上限値を考慮した獲得関数です。

\begin{align} a_{UCB} = μ + \sqrt{\frac{\log N}{N}} σ \end{align}

データ密度の低い所を探していきます。

def ucb(X, X_sample, y_sample, gpr):
    mu, sigma = gpr.predict(X, return_std=True)
    mu_sample = gpr.predict(X_sample)

    n_sample = X_sample.shape[0]
    mu_sample_opt = np.max(mu_sample)
    z = mu_sample_opt + np.sqrt(np.log(n_sample) / n_sample) * sigma
    return z

Expected Improvement

思考実験由来のものです。ベイズ最適化で最もよく使用されています。

\begin{align} EI(x) = \mathbb{E} [\max (f(x) - f(x^{+}), 0)] \end{align}

ここで $f(x^{+})$ は最も良いサンプル値です。 $x^{+} = \arg \max_{x_{i} \in x_{1:t}} f(x_{i})$ です。GPのモデルは以下の式になります。

\begin{align} EI(x) = \begin{cases} (μ\mathbf{x}) - f(\mathbf{x}^{+}) - \xi) \Phi(Z) + σ(\mathbf{x})\phi(Z) &\text{if} \ σ(\mathbf{x}) > 0 \\ 0 & \text{if}\ σ(\mathbf{x}) = 0 \end{cases} \end{align}

ここで、Zは次のようになります。

\begin{align} Z = \begin{cases} \frac{μ(\mathbf{x}) - f(\mathbf{x}^{+}) - \xi}{σ(\mathbf{x})} &\text{if}\ σ(\mathbf{x}) > 0 \\ 0 & \text{if}\ σ(\mathbf{x}) = 0 \end{cases} \end{align}

また、 $\Phi, \phi$ はそれぞれ正規分布のCDFとPDFです。

def expected_improvement(X, X_sample, Y_sample, gpr, xi=0.01):
    mu, sigma = gpr.predict(X, return_std=True)
    mu_sample = gpr.predict(X_sample)
    sigma = sigma.reshape(-1, X_sample.shape[1])

    mu_sample_opt = np.max(mu_sample)

    with np.errstate(divide='warn'):
        imp = mu - mu_sample_opt - xi
        Z = imp / sigma
        ei = imp * norm.cdf(Z) + sigma * norm.pdf(Z)
        ei[sigma == 0.0] = 0.0

    return ei

Knowledge Gradient, KG

knowledge-gradientな獲得関数は、以前の評価点が最終的な解としてEIを最大化させようとします。この考え方は、評価がノイズフリーでリスクを嫌うときは良い方法となります。ノイズが含まれているときは提案点にノイズが含まれるため評価が困難となります。Frazier先生が提案している方法です。

\begin{align} KG (x) = \mathbb{E} [ μ_{n+1} - μ_{n} | x_{n+1} = x] \end{align}

こちらのアルゴリズムの最初はFranzier, 2012です。

f:id:udnp:20190420235207p:plain

Entropy Search

Entropy serach (ES)は情報理論に基づく方法です。

\begin{align} ES(x) = H(P(x)) - E[H(P(x|f(x)))] \end{align}

となっています。KG, ESも実装して試してみたいです。

Exotic Bayesian Optimization

上記の方法論はスタンダードなベイズ最適化です。超長方形やシンプレックス、勾配情報の不足、ノイズフリーな評価法といった前提を仮定しています。実際にはこれらの前提条件が崩れることがあり、それを考慮した方法が提案されています。

獲得関数の設計が大事になってきそうです。KGやESは下記の条件を満たせるように獲得関数を設定できます。

Noisy Evaluations
Parallel Evaluations
Constraints
Multi-Fidelity & Multi-Information Source Evaluations
Random Environmental Conditions and Multi-Task Bayesian Optimization
Derivative Observations

より適切にベイズ最適化を行うためには、こちらも詳しくサーベイしていく必要がありそうです。

ソフト

RとPythonとMatlabのものが紹介されています。わりとたくさんあります。

DiceKriging
laGP

Python

Metrics Optimization Engine
Spearmint
GpyOpt
GPFlow
GPyTorch

など。

UCBとEIの例

とりあえずスタンダードの例でテストしてみます。下記関数を考えます。

\begin{align} f(x) = e^{-(x - 2)^{2}} + e^{\frac{-(x - 6)^{2}}{10}} + \frac{1}{x^{2} + 1} + noise \end{align}

f:id:udnp:20190421134812g:plain

import numpy as np
import matplotlib.pyplot as plt
from bayes_opt import BayesianOptimization
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel
from scipy.stats import norm
from scipy.optimize import minimize

noise = 0
bounds = np.array([[-2.0, 10.0]])
def f(x, noise=noise):
    return np.exp(-(x - 2) ** 2) + np.exp(-(x - 6) ** 2 / 10) + 1 / (x ** 2 + 1) + noise * np.random.randn(*x.shape)

def plot_result(X, X_init, y_init, gpr, itr):
    y_mean, y_cov = gpr.predict(X, return_std=True)
    upr = y_mean - np.sqrt(y_cov.reshape(-1, 1))
    lwr = y_mean + np.sqrt(y_cov.reshape(-1, 1))

    plt.clf()
    plt.subplot(211)
    plt.title("Bayesian Optimization # {}".format(i+1))
    plt.plot(X_init, y_init, "ro", label="Initial samples")
    plt.plot(X, y_mean, label="Surrogate model")
    plt.plot(X, f(X), label="Objective")
    plt.xlim([-2.0, 10.0])
    plt.ylim([-0.5, 1.5])
    plt.fill_between(X.ravel(), upr.ravel(), lwr.ravel(), alpha=0.5)
    plt.legend()
    plt.subplot(212)
    #plt.plot(X, ucb(X, X_init, y_init, gpr), label="UCB")
    plt.plot(X, expected_improvement(X, X_init, y_init, gpr), label="EI")
    plt.xlim([-2.0, 10.0])
    plt.legend()
    plt.savefig("./BayesOptEI{}.png".format(itr))

def next_X(acq, X_sample, Y_sample, gpr):
    dim = X_sample.shape[1]
    # 獲得関数の最大値を求めるために負にしています。
    def min_obj(X):
        return -acq(X.reshape(-1, dim), X_sample, Y_sample, gpr)

    # 局所解に陥るので、初期値を変えます。PSOやGAなどの方が良いかもしれません。
    for x0 in np.random.uniform(bounds[:, 0], bounds[:, 1], size=(25, dim)):
        res = minimize(min_obj, x0=x0, bounds=bounds, method='L-BFGS-B')
        if res.fun < 1:
            min_val = res.fun[0]
            min_x = res.x
    X_next = min_x.reshape(-1, 1)
    return X_next

def main():
    X_init = np.array([[-0.9], [1.1]])
    y_init = f(X_init)
    kernel = ConstantKernel(1.0) + RBF(length_scale=2.0, length_scale_bounds=(0, 10))
    gpr = GaussianProcessRegressor(kernel=kernel, random_state=42).fit(X_init, y_init)
    for i in range(10):
        gpr = GaussianProcessRegressor(kernel=kernel, random_state=42).fit(X_init, y_init)
        plot_result(X, X_init, y_init, gpr, i)
        X_next = next_X(expected_improvement, X_init, y_init, gpr)
        X_init = np.r_[X_init, X_next]
        y_init = np.r_[y_init, f(X_next)]

なおパッケージを使えば簡単です。

def plot_bo(f, bo):
    x = np.linspace(-2, 10, 10000)
    mean, sigma = bo._gp.predict(x.reshape(-1, 1), return_std=True)

    plt.figure(figsize=(16, 9))
    plt.plot(x, f(x))
    plt.plot(x, mean)
    plt.fill_between(x, mean + sigma, mean - sigma, alpha=0.1)
    plt.scatter(bo.space.params.flatten(), bo.space.target, c="red", s=50, zorder=10)
    plt.show()

def BayesRun():
    bo = BayesianOptimization(
        f=f,
        pbounds={"x": (-2, 10)},
        verbose=0,
        random_state=987234,
    )
    bo.maximize(n_iter=10, acq="ucb", kappa=0.1)
    plot_bo(f, bo)

結果

|   iter    |  target   |     x     |
-------------------------------------
|  6        |  1.305    |  1.651    |
|  9        |  1.401    |  1.967    |
|  12       |  1.402    |  2.019    |
|  13       |  1.402    |  2.0      |
|  14       |  1.402    |  2.0      |
=====================================

f:id:udnp:20190421141307p:plain