mlrMBOを用いたベイズ最適化（Bayesian Optimization）

1.はじめに

Black-box関数最適化問題というのは，目的関数の代数的な表現，つまり，数式自体が与えられず，解に対する目的関数値のみを利用することができる問題で、これを解決する方法の一つにベイズ最適化（Bayesian Optimization）があります。Black-box関数最適化では解の1回の評価の計算コストが大きいことを前提しています。
例えば農作物の栽培試験では、栽培してから収穫まで時間がかかるので、実験回数が限られてしまい、1回の評価コストが大きくなります。
そこで、できるだけ少ない評価回数でベストなパラメータを推定する方法にベイズ最適化が用いられ、Rでは"rBayesianOptimization"パッケージが提供されています。
しかし、mlrで提供されている同様のmlrMBOパッケージについては、日本語による情報があまりありませんので、今回はmlrMBOパッケージによるベイズ最適化を試してみたいと思います。

2.ベイズ最適化とは

ベイズ最適化では評価に大きな計算コストのかかる目的関数の代わりに、ガウス過程回帰と獲得関数と呼ばれる代理関数を最適化することで，次の探索点を選択します。

• ガウス過程回帰によって予測した平均，分散を使って計算した獲得関数を最適化して，次の探索点を得る
• 得られた探索点を評価し，評価値を得る
• 探索点と評価値の組を保存し，ガウス過程回帰を更新する

ガウス過程回帰では、関数の事後分布を予測する（平均と分散）。上記図では青色の部分。
図下の緑色の部分は事後分布から計算された獲得関数(acquisition function)。
獲得関数が最大値をとる値が次の観測値の候補となる。

3.mlrMBOを使ったベイズ最適化

3.1　パラメータと最大値の推定

目的となる未知の1次関数を作成します。
この関数の形はわかりませんが、パラメータがとりうる範囲内で最大値を推定したいと思います。

library(mlrMBO)
target = makeSingleObjectiveFunction(
  name = "未知の1次関数",
  fn = function(x) {
  x = x * 10
  x * sin(x) 
  },
  par.set = makeParamSet(
    makeNumericParam("X", lower = 0, upper = 1.5)
  ),
  minimize = FALSE
)
plot(target)

Xがパラメータでｙを最大にする条件を求めます。
Xは0~1.5の範囲をとります。

主に使われることが多い獲得関数は次の２つですが、mlrMBOでは次のとおり用意されています。

まず獲得関数（LCB)とイテレーションを6回、定義します。

control = makeMBOControl()
control = setMBOControlTermination(control, iters = 6)
control = setMBOControlInfill(control, crit = makeMBOInfillCritCB())

主に使われることが多い獲得関数は次の２つです（と思っている）。

• Expected Improvement (EI)
• Lower Confidence Bound (LCB)

詳しくは下記を参考にしてください。

シンプルなベイズ最適化について

mlrMBOではその他にも次のとおり獲得関数が用意されています。

• makeMBOInfillCritMeanResponse()
• makeMBOInfillCritStandardError()
• makeMBOInfillCritEI(se.threshold = 1e-06) 　　　Expected Improvement (EI)
• makeMBOInfillCritCB(cb.lambda = NULL) 　　　 Lower Confidence Bound (LCB)
• makeMBOInfillCritAEI(aei.use.nugget = FALSE, se.threshold = 1e-06)
• makeMBOInfillCritEQI(eqi.beta = 0.75, se.threshold = 1e-06)
• makeMBOInfillCritDIB(cb.lambda = 1, sms.eps = NULL)
• makeMBOInfillCritAdaCB(cb.lambda.start = NULL, cb.lambda.end = NULL)

ガウス過程回帰を定義します。
この関数の細かいところはよくわからないので、デフォルトで使っています。

surr.km = makeLearner("regr.km", predict.type = "se",)

それでは、走らせてみます。

set.seed(12)
run1 = mbo(target, learner = surr.km, control = control, show.info = TRUE )

$>Computing y column(s) for design. Not provided.
$>[mbo] 0: X=1.08 : y = -10.7 : 0.0 secs : initdesign
$>[mbo] 0: X=0.207 : y = 1.82 : 0.0 secs : initdesign
$>[mbo] 0: X=1.3 : y = 5 : 0.0 secs : initdesign
$>[mbo] 0: X=0.734 : y = 6.38 : 0.0 secs : initdesign
$>
$>optimisation start
$>------------------
$>* estimation method   : MLE 
$>* optimisation method : BFGS 
$>* analytical gradient : used
$>* trend model : ~1
$>* covariance model : 
$>  - type :  matern5_2 
$>  - nugget : NO
$>  - parameters lower bounds :  1e-10 
$>  - parameters upper bounds :  2.178885 
$>  - best initial criterion value(s) :  -13.42582 
$>　~ 以下省略 ~
$>iterations 6
$>function evaluations 10
$>segments explored during Cauchy searches 6
$>BFGS updates skipped 0
$>active bounds at final generalized Cauchy point 0
$>norm of the final projected gradient 1.42266e-05
$>final function value 26.7622
$>
$>F = 26.7622
$>final  value 26.762208 
$>converged
$>[mbo] 6: X=1.43 : y = 14.2 : 0.0 secs : infill_cb
$> ~ 以下省略 ~

イテレーション6回目でX=1.43のときにy=14.2を推定しました。
逐次的に推定する過程を図で示します。

set.seed(123)
run2 = exampleRun(target, learner = surr.km, 
                 control = control, 
                 show.info = FALSE)
plotExampleRun(run2, iters = c(1,3,6), pause = FALSE)

6回のイテレーションで未知の関数に対してほぼ正確に最大値を推定できています。
Initial Designを指定していないため、デフォルトで４点が初期ポイントとして自動的にランダムに指定されます（赤い点）。
初期ポイントの配置によっては、真ん中の極大値にはまってしまうことがあります。

3.2　関数を定義しない場合のパラメータの提案

ここからは、農作物の架空のデータを使います。
4年間の栽培データです。NがチッソでKがカリの肥料です。yが反当たりの収量です。
5年目の栽培をするに当たって、収量を最大化するにはNとKをどれぐらい投入すればよいのかを推定します。
全く目的関数を定義することなしに全て結果を手入力します。

N	K	y
10	2	20
15	4	25
30	5	30
25	7	32

ps = makeParamSet(
  makeNumericParam("N", lower = 1, upper = 50),
  makeNumericParam("K", lower = 1, upper = 20)
)
des = data.frame("N"=c(10,15,30,25),"K"=c(2,4,5,7))
des$y = c(20,25,30,32)

パラメータがとりうる範囲を指定し、desに表データを手入力しました。

獲得関数にCBを指定します。

ctrl = makeMBOControl()
ctrl = setMBOControlInfill(ctrl, crit = crit.cb)

次の肥料の投入量の推定を行います。

set.seed(123)
opt.state = initSMBO(par.set = ps, design = des, control = ctrl, minimize = F, noisy = F)
plot(opt.state)

左側のcbの青色が濃い部分のところのパラメータの組み合わせが示唆されています。

proposePoints(opt.state)

$>$prop.points
$>           N        K
$>347 43.56749 9.218127
$>
$>$propose.time
$>[1] 0.61
$>
$>$prop.type
$>[1] "infill_cb"
$>
$>$crit.vals
$>          [,1]
$>[1,] -35.60235
$>
$>$crit.components
$>        se     mean lambda
$>1 4.089923 31.51242      1
$>
$>$errors.model
$>[1] NA
$>
$>attr(,"class")
$>[1] "Proposal" "list"    

提案された値のNが43.56749、Kが9.218127kg/10aで収穫は最大化するそうです。
5年目の肥料の投入量は、この辺りを試すとよさそうです。
ベイズ最適化は機械学習のパラメータ推定以外にも使えそうですね。

最後に、何か間違っていることがあるかもしれませんので、ご指摘いただけたら大変ありがたいです。

4.参考

シンプルなベイズ最適化について
mlrMBO
A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning