はじめに

ベイズ最適化の話について、興味が沸いたので簡単にまとめます。 (ただのメモです)

詳しくはこちらなどを見ると良いと思います。

cyberagent.ai

ベイズ最適化

ベイズ最適化とは，未知の目的関数𝑓(𝑥)を最大化する $X^\ast$ を求める大域最適化の手法です。

$X^\ast = argmax_{𝑋∈𝜒} f(X)$

ここで、 $𝜒$ は探索空間の変数集合であり、離散値でも連続値でも良い。

新しいデータの獲得する候補を評価する獲得関数 $\alpha$ を用いる。

現在までに得られているデータ $D_n$ から獲得関数を最大化する $X_{n+1}$ について調べる。

$X_{n+1} = argmax_{𝑋∈𝜒} \alpha (X; D_n)$

そこから観測値 $y_{ n+1 } = f ( X_{n+1} )$ を取得し、新たなデータとして $D_n$ に加える。

具体的には、次のように逐次的にデータに加えていく。

$D_{n+1} = \{ D_n, ( X_{ n+1 } , y_{ n+1 } ) \}$

その後，獲得関数 $\alpha$ を最大にする $X$ を選び $y$ を観測することを繰り返しで最適値の探索を行う。

アルゴリズムは次のように表現される。

f:id:saltcooky:20200826022622p:plain

ガウス過程の導入

ベイズ最適化の枠組みでは，統計モデル $f(x)$ としてガウス過程(Gaussian process)が用いられる。

ベイス推定を用いて一つの関数を推定するのではなく、関数の分布を推定することができる。

<あとで追記>

得られた観測値の各点 $x_i$ はガウス分布 $N(\mu_i, \sigma_i)$ に従っており、それが横に並んでいると考えられる。

$X _{ 1:N }$ は、今までに観測されている値である。

するとN点が観測された $X_{1:N}$ の事後分布(予測分布)が[ tex: x _{N+1} ] に対して定まる。

この事後分布は、次の平均 $\mu$ と分散 $\sigma$ からなる正規分布となる。

$\mu_{N+1}= k^{T}(K+\beta I)^{-1} X_{1:N}$

$\sigma_{N+1}= k(x_{N+1},x_{N+1})-k^{T}(K+\beta I)^{-1}k$

$k()$ はカーネル関数、Kはカーネルによるグラム行列である。

この値を用いて獲得関数の更新を行うことになる。

用いられるカーネル

カーネルには様々なものが存在する。

よく用いられるカーネルは次のようなものである。

Squared Exponential Function

$k({\bf x}, {\bf x}^{\prime}) = \exp \Bigl( -\frac{1}{2}|| {\bf x} – {\bf x}^{\prime} ||^2 \Bigr)$

Gaussian Kernel

$k({\bf x}, {\bf x}^{\prime}) = \exp \Bigl( - || {\bf x} – {\bf x}^{\prime} ||^2 / 2\sigma^2\Bigr)$

Ornstein-Uhlenbeck Process

$k({\bf x}, {\bf x}^{\prime}) = \exp \Bigl( -\theta | {\bf x} – {\bf x}^{\prime} | \Bigr)$

探索と活用

最適値の探索では獲得関数 $\alpha$ が最大の点を探索することになる。

この獲得関数による観測点を決定する戦略は、複数提案されている。

基本的には過去のデータに基づいて平均的に最大となるもの（活用）とデータの少ない未知の領域について調べるもの（探索）の二種類である。

Probability of Improvement(PI)

PIの獲得関数 $\alpha$ は次のようになる。

$Z = \frac {\mu(x)-f^{+}- \xi}{\sigma(x)} \\ a_{PI}(x)= \Phi (Z)$

ここで、 $f^{+}$ は今までに観測された値の最大値、 $\Phi (Z)$ は正規分布の累積密度関数、 $\xi$ は探索と活用のどちらを重視するかのトレードオフパラメータ。

現在の最大値を超える確率が高い点を次に観測することを意味する。

局所解に落ちいる可能性がある。

Expected Improvement(EI)

EIの獲得関数 $\alpha$ は次のようになる。

$a_{EI}(x)=(\mu (x)-f^{+}-\xi)\Phi (Z)+\sigma (x)\phi (Z)$

評価値と現在の最大値の差の期待値が最も⾼くなる点を次に観測することを意味する。

Upper-Confidence Bound(UCB)

UCBの獲得関数 $\alpha$ は次にようになる。

$a_{UCB}(x)=\mu (x)+k_{N} \sigma (x)$

ここで、 $k$ は探索と活用のトレードオフを設定するパラメータである。

評価値の信頼区間の上限が最も高い点を次に観測することを意味する。

最適解を求めることができる理論的保証が存在する。

Rでベイズ最適化

rBayesianOptimization packageとmlrMBO packageを用いた方法のそれぞれを確認してみた。

rBayesianOptimization

リポジトリはこちら。

github.com

ここで用いるのは、オークランドの火山地帯にある火山の1つであるVolcano火山の地形情報であるvolcanoです。

このデータを用いて最も標高が高い点を探索してみます。

なお、獲得関数にはUCBを用いた。

library(rBayesianOptimization)

x <- (1:nrow(volcano))
y <- (1:ncol(volcano))
image(x, y, volcano, col = terrain.colors(100), axes = FALSE)

# 探索したい関数の定義
volcano_func <- function(x, y){
  Pred <- volcano[x,y]
  tmp <- list(Score=Pred, Pred=Pred)
  return(tmp)
}

set.seed(123)
LOOP <- 50

BO_Res <-
  BayesianOptimization(volcano_func, 
                       bounds = list(x=c(1L,87L), y=c(1L,61L)), 
                       init_points = LOOP, 
                       n_iter = 1, 
                       acq = "ucb",
                       kappa = 1.5, 
                       eps = 0, 
                       verbose = TRUE)

推定された最適値の情報を確認する。

> BO_Res$Best_Par
 x  y 
21 35 
> BO_Res$Best_Value
[1] 190

最適ではないが、ある程度最適値に近い値を探索することができている。

(明るい色の方が標高が高い)

f:id:saltcooky:20200831225233p:plain

mlrMBO

ここで探索に利用するデータは、次のような関数により生成されたデータです。

$f(x) = x^2 + 8 sin(1.98 \pi x) * cos(1.2 \pi x) , ( -0.5 \le x \le 0.5)$

mlrMBO packageを用いた方法が次のような物である。

# 探索対象の関数を定義
fun = function(x) {
  x <- as.numeric(x)
  x^2 + 8*sin(1.98 * pi * x) * cos(1.2 * pi * x)
}

# 探索するパラメータ空間を定義
ps = makeParamSet(
  makeNumericParam("x", lower = -5, upper = 5)
)

# smoof関数でラップする
fn = makeSingleObjectiveFunction(
  fn = fun,
  par.set = ps,
  has.simple.signature = FALSE
)

# 初期探索点を取得
des = generateDesign(n = 3, par.set = ps)

# 対応する関数の出力を取得
des$y = apply(des, 1, fn)

# MBOの初期化
ctrl = makeMBOControl()
ctrl = setMBOControlInfill(ctrl, crit = crit.ei)
ctrl = setMBOControlTermination(ctrl, iters = 20)

# learnerオブジェクトの作成
surr.km = makeLearner("regr.km",
                      predict.type = "se", 
                      covtype = "matern3_2", 
                      control = list(trace = FALSE))

# 
configureMlr(show.info = FALSE, show.learner.output = FALSE)

# 最適化処理
run = mbo(fun = fn, 
          design = des,
          learner = surr.km,
          control = ctrl)

# 
run = exampleRun(fn, 
                 learner = surr.km, 
                 control = ctrl, 
                 show.info = FALSE)

# グラフ化
plotExampleRun(run, iters = c(1L:20L), pause = F)