渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

{rBayesianOptimization}パッケージによるベイズ最適化で機械学習パラメータチューニングをお手軽に

機械学習のパラメータチューニングというと大なり小なり大変な部分があって、今年のエイプリルフール記事に皆さん引っかかって下さったところを見るにパラメータチューニングを簡単に済ませたい!と願う人々は世の中多いようです(笑)。


少し前のMXnetを使った記事でも取り上げましたが、そのパラメータチューニングを迅速に済ませようというアイデアの一つがベイズ最適化(Bayesian Optimization)です。

要は、グリッドサーチのように網羅的に最適なパラメータを探索しに行くのではなく、一つのパラメータで精度をチェックしたらその次は精度が上がりやすそうな方向にベイズ的に逐次改善を行いながら探索していく、という方法のことです。


世の中色々seminar paper的なものがあるようですが、arXivから@さんが見つけてきて下さったのがこれ。

日本語のブログ記事なら、id:olanleedさんのこちらの記事が良いと思います。


ということでこれらのreferencesを紹介してしまえば理論的なポイントは全ておしまいなのですが、リンク先を読むのも面倒という人のためにごくごく簡単な説明だけを書いておきます。ただし僕は数学が大の苦手なので数式は紹介する気もない上に、そもそも上記のreferencesからの受け売り以外何もありません(笑)。


ベイズ最適化とは


f:id:TJO:20160815120951p:plain

上記のarXiv論文のこのFig.1が最も分かりやすいかと思います(ちょうど最適化したい関数も1次元でここで想定している「機械学習のCV精度」の最適化というテーマとマッチするので)。


端的に言えば「機械学習のCV精度をパラメータを変えながら改善していくという試行錯誤のプロセスを、ベイズ統計の考え方を活用して効率化したもの」です。特にここで挙げられている方法は「機械学習のCV精度がある未知のガウス過程f(x)に従う」と仮定した上で、これに対して事前に決め打ちで選んだ「f(x)の期待値を算出してくれる関数u(\cdot)」を用意し、交互に次のx次のxと反復繰り返し更新していくことで、f(x)を最大にするxを探し出す、というアルゴリズムです。


その様子が上の図に分かりやすく書かれています。まず初期値を出したいのでランダムに2点選んできて黒線で表されるf(x)の値が2つ求まったところから始まります(上段)。これに対してu(\cdot)が緑線で表される関数としてこのステップで求まります。そこでu(\cdot)の中をある程度ランダムに探索した中で最も高い値を与えた点(赤い▼のところ)を次のサンプリング点として、f(x)の値を求めます(中段)。この結果を受けてu(\cdot)がさらに更新されてそのまた次の赤い▼が決まるので、これをそのまた次のサンプリング点としてf(x)の値を求めます(下段)。これを繰り返すことで、徐々にf(x)を最大化するxへと近付いていく、というわけです。


以上の手続きを擬似コードにしたのが論文中のAlgorithm 1です。
f:id:TJO:20160815123627p:plain


こう書くと簡単そうに見えますが、色々な理屈が裏側にはあります。まずガウス過程に共分散を与えるためにカーネルを事前に決めておく必要があります(細かい話は「黄色い本」ことPRMLの6.4節を参照のこと)。次にu(\cdot)についてはProbaiblity of Improvement (PI), Expected Improvement (EI), GP Upper Confidence Bound (GP-UCB)という3つのアルゴリズムが提案されており、いずれかを選ぶ必要があります。ベイズ最適化で機械学習パラメータチューニングを行う場合は、この辺のベイズ最適化自体のチューニングにも依存するという点に注意が必要です。


とりあえず、これ以上は普通に上記の関連資料をご覧いただいた方が遥かに手っ取り早い&正確なので割愛します。機械学習パラメータチューニングに話を絞らなければ色々と広がりのあるテーマなので個人的には面白そうだと思うのですが、数学大の苦手人間としては深追いはしないでおきます(笑)。


{rBayesianOptimization}をインストールする



普通にCRANに上がっているので、stable versionならinstall.packagesで入ります。β版を試したければdevtools::install_githubで入れればOKです。使い方も簡単です。

BayesianOptimization(Test_Fun, bounds, init_points, n_iter, acq = "ucb",
  kappa = 2.576, eps = 0, verbose = TRUE, ...)

「最大化」したい目的関数をR functionとして書き下しておいて、これを第1引数Test_Funに入れます。この目的関数は必ず"Score", "Pred"という2つの返値をリストとして持っている必要があります。一応"Score"が最大化するターゲットなので、ここにターゲットにしたいパフォーマンス(例えばテストデータに対するパフォーマンス)を代入しておきます。"Pred"にはCVに対するパフォーマンスを代入することになっていますが、いまいち使い道がわからないので僕はScoreと同じものを入れています。そしてBayesianOptimization自体が標準出力にレポートを吐くので、Test_Funの方はverbose = FALSEにしておくことをお薦めします。


これだけなんですが、このベイズ最適化関数は「最大化」するものなので、例えばRMSEなど「最小化」したいパフォーマンスをターゲットにする場合は返値を定義する際に符号を負にしておきましょう。


SVMで試してみる


では、早速試してみましょう。まずはSVM。例えばですが、MNISTに対する多項式カーネルの次数を決める課題でやってみようと思います。

> train <- read.csv('mnist/short_prac_train.csv')
> test <- read.csv('mnist/short_prac_test.csv')
> train$label <- as.factor(train$label)
> test$label <- as.factor(test$label)

> library(e1071)
> library(rBayesianOptimization)

> svm_holdout <- function(deg, c){
+     model <- svm(label~., train, kernel='polynomial', degree=deg, coef0=c)
+     t.pred <- predict(model, newdata=test[,-1])
+     Pred <- sum(diag(table(test$label, t.pred)))/nrow(test)
+     list(Score=Pred, Pred=Pred)
+ }

> opt_svm <- BayesianOptimization(svm_holdout,
+                                 bounds=list(deg=c(3L,10L),
+                                             c=c(0,10)),
+                                 init_points=10, n_iter=1, acq='ei', kappa=2.576,
+                                 eps=0.0, verbose=TRUE)
elapsed = 30.02	Round = 1	deg = 8.0000	c = 7.6517	Value = 0.8210 
elapsed = 34.01	Round = 2	deg = 9.0000	c = 3.1119	Value = 0.7880 
elapsed = 32.12	Round = 3	deg = 8.0000	c = 5.1357	Value = 0.8210 
elapsed = 29.32	Round = 4	deg = 6.0000	c = 9.5419	Value = 0.8740 
elapsed = 32.69	Round = 5	deg = 8.0000	c = 0.8852	Value = 0.8220 
elapsed = 23.57	Round = 6	deg = 4.0000	c = 5.0165	Value = 0.9270 
elapsed = 28.99	Round = 7	deg = 7.0000	c = 3.8474	Value = 0.8480 
elapsed = 23.23	Round = 8	deg = 4.0000	c = 1.4057	Value = 0.9270 
elapsed = 28.99	Round = 9	deg = 7.0000	c = 2.1059	Value = 0.8480 
elapsed = 33.52	Round = 10	deg = 9.0000	c = 3.6292	Value = 0.7890 
elapsed = 20.74	Round = 11	deg = 3.0000	c = 0.0000	Value = 0.9510 

 Best Parameters Found: 
Round = 11	deg = 3.0000	c = 0.0000	Value = 0.9510 

3次の多項式カーネルでcoef0 = 0.0000として、95.1%のtest ACCが得られました。


ランダムフォレストで試してみる

> library(randomForest)
> rf_holdout <- function(mnum){
+     model <- randomForest(label~., train, mtry=mnum)
+     t.pred <- predict(model, newdata=test[,-1])
+     Pred <- sum(diag(table(test$label, t.pred)))/nrow(test)
+     list(Score=Pred, Pred=Pred)
+ }
> opt_rf <- BayesianOptimization(rf_holdout,
+                                 bounds=list(mnum=c(3L,40L)),
+                                 init_points=10, n_iter=1, acq='ei', kappa=2.576,
+                                 eps=0.0, verbose=TRUE)
elapsed = 152.45	Round = 1	mnum = 33.0000	Value = 0.9460 
elapsed = 142.13	Round = 2	mnum = 22.0000	Value = 0.9510 
elapsed = 137.53	Round = 3	mnum = 11.0000	Value = 0.9500 
elapsed = 130.25	Round = 4	mnum = 4.0000	Value = 0.9340 
elapsed = 151.26	Round = 5	mnum = 38.0000	Value = 0.9500 
elapsed = 138.97	Round = 6	mnum = 16.0000	Value = 0.9500 
elapsed = 145.66	Round = 7	mnum = 29.0000	Value = 0.9510 
elapsed = 144.63	Round = 8	mnum = 11.0000	Value = 0.9530 
elapsed = 151.59	Round = 9	mnum = 23.0000	Value = 0.9460 
elapsed = 141.68	Round = 10	mnum = 21.0000	Value = 0.9510 
elapsed = 167.45	Round = 11	mnum = 28.0000	Value = 0.9530 

 Best Parameters Found: 
Round = 8	mnum = 11.0000	Value = 0.9530

デフォルト値のmtry = 28(特徴次元数の平方根)ではなく、mtry = 11でtest ACC 95.3%という結果になりました。


Xgboostで試してみる

> library(xgboost)
> library(xgboost)
> library(Matrix) # データの前処理に必要
> train.mx<-sparse.model.matrix(label~., train)
> test.mx<-sparse.model.matrix(label~., test)
> # データセットをxgboostで扱える形式に直す
> dtrain<-xgb.DMatrix(train.mx, label=as.integer(train$label)-1)
> dtest<-xgb.DMatrix(test.mx, label=as.integer(test$label)-1)
> xgb_holdout <- function(ex, mx, nr){
+     model <- xgb.train(params=list(objective="multi:softmax", num_class=10,
+                                     eval_metric="mlogloss",
+                                     eta=ex/10, max_depth=mx),
+                                     data=dtrain, nrounds=nr)
+     t.pred <- predict(model, newdata=dtest)
+     Pred <- sum(diag(table(test$label, t.pred)))/nrow(test)
+     list(Score=Pred, Pred=Pred)
+ }
> opt_xgb <- BayesianOptimization(xgb_holdout,
+                                 bounds=list(ex=c(2L,5L), mx=c(4L,5L), nr=c(70L,160L)),
+                                 init_points=20, n_iter=1, acq='ei', kappa=2.576,
+                                 eps=0.0, verbose=TRUE)
elapsed = 26.36	Round = 1	ex = 5.0000	mx = 5.0000	nr = 96.0000	Value = 0.9390 
elapsed = 36.38	Round = 2	ex = 3.0000	mx = 5.0000	nr = 117.0000	Value = 0.9550 
elapsed = 41.42	Round = 3	ex = 3.0000	mx = 5.0000	nr = 142.0000	Value = 0.9530 
elapsed = 34.00	Round = 4	ex = 4.0000	mx = 5.0000	nr = 125.0000	Value = 0.9560 
elapsed = 27.23	Round = 5	ex = 3.0000	mx = 4.0000	nr = 94.0000	Value = 0.9540 
elapsed = 43.66	Round = 6	ex = 2.0000	mx = 4.0000	nr = 152.0000	Value = 0.9570 
elapsed = 24.46	Round = 7	ex = 4.0000	mx = 5.0000	nr = 75.0000	Value = 0.9410 
elapsed = 29.90	Round = 8	ex = 3.0000	mx = 5.0000	nr = 89.0000	Value = 0.9460 
elapsed = 36.41	Round = 9	ex = 3.0000	mx = 4.0000	nr = 127.0000	Value = 0.9510 
elapsed = 32.22	Round = 10	ex = 2.0000	mx = 5.0000	nr = 83.0000	Value = 0.9520 
elapsed = 39.67	Round = 11	ex = 4.0000	mx = 5.0000	nr = 149.0000	Value = 0.9460 
elapsed = 32.22	Round = 12	ex = 4.0000	mx = 5.0000	nr = 103.0000	Value = 0.9510 
elapsed = 30.76	Round = 13	ex = 4.0000	mx = 5.0000	nr = 97.0000	Value = 0.9440 
elapsed = 38.54	Round = 14	ex = 4.0000	mx = 5.0000	nr = 147.0000	Value = 0.9460 
elapsed = 46.11	Round = 15	ex = 2.0000	mx = 5.0000	nr = 136.0000	Value = 0.9500 
elapsed = 30.57	Round = 16	ex = 4.0000	mx = 5.0000	nr = 105.0000	Value = 0.9480 
elapsed = 24.91	Round = 17	ex = 3.0000	mx = 4.0000	nr = 85.0000	Value = 0.9470 
elapsed = 39.63	Round = 18	ex = 4.0000	mx = 5.0000	nr = 159.0000	Value = 0.9410 
elapsed = 23.48	Round = 19	ex = 2.0000	mx = 4.0000	nr = 77.0000	Value = 0.9410 
elapsed = 30.29	Round = 20	ex = 4.0000	mx = 4.0000	nr = 116.0000	Value = 0.9480 
elapsed = 33.45	Round = 21	ex = 4.0000	mx = 5.0000	nr = 122.0000	Value = 0.9420 

 Best Parameters Found: 
Round = 6	ex = 2.0000	mx = 4.0000	nr = 152.0000	Value = 0.9570 

実は以前出したハイスコア(95.8%)には及ばないんですが、test ACC 95.7%に達しました。


Deep Learning (DNN)で試してみる


これは先日の記事でやった通りですので、ここではそれを再掲するに留めておきます。MXnetでオンラインニュース人気度のデータを回帰したものですね。ただし一点大きな修正点がありまして、Test_FunのところでScore返値に思いっ切り学習誤差を入れてしまっていたので、テスト誤差を入れるようにしてありますorz なので結果も変わっております。。。

> d <- read.csv('OnlineNewsPopularity.csv')
> d <- d[,-c(1,2)]
> idx <- sample(nrow(d),5000,replace=F)
> d_train <- d[-idx,]
> d_test <- d[idx,]
> library(mxnet)
> train <- data.matrix(d_train)
> test <- data.matrix(d_test)
> train.x <- train[,-59]
> train.y <- train[,59]
> test.x <- test[,-59]
> test.y <- test[,59]
> train_means <- apply(train.x, 2, mean)
> train_stds <- apply(train.x, 2, sd)
> test_means <- apply(test.x, 2, mean)
> test_stds <- apply(test.x, 2, sd)
> train.x <- t((t(train.x)-train_means)/train_stds)
> test.x <- t((t(test.x)-test_means)/test_stds)
> train.y <- log(train.y)
> test.y <- log(test.y)
> mxnet_holdout_bayes <- function(unit1, unit2, unit3, num_r, learn_r){
+     data <- mx.symbol.Variable("data")
+     fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=unit1)
+     act1 <- mx.symbol.Activation(fc1, name="tanh1", act_type="tanh")
+     drop1 <- mx.symbol.Dropout(act1, p=0.2)
+     fc2 <- mx.symbol.FullyConnected(drop1, name="fc2", num_hidden=unit2)
+     act2 <- mx.symbol.Activation(fc2, name="tanh2", act_type="tanh")
+     drop2 <- mx.symbol.Dropout(act2, p=0.2)
+     fc3 <- mx.symbol.FullyConnected(drop2, name="fc3", num_hidden=unit3)
+     act3 <- mx.symbol.Activation(fc3, name="tanh3", act_type="tanh")
+     drop3 <- mx.symbol.Dropout(act3, p=0.2)
+     fc4 <- mx.symbol.FullyConnected(drop3, name="fc4", num_hidden=1)
+     output <- mx.symbol.LinearRegressionOutput(fc4, name="linreg")
+     devices <- mx.cpu()
+     mx.set.seed(71)
+     model <- mx.model.FeedForward.create(output, X=train.x, y=train.y, 
+                                          ctx=devices, num.round=num_r, array.batch.size=200,
+                                          learning.rate=learn_r, momentum=0.99,
+                                          eval.metric=mx.metric.rmse,
+                                          initializer=mx.init.uniform(0.5),
+                                          array.layout = "rowmajor",
+                                          epoch.end.callback=mx.callback.log.train.metric(20),
+                                          verbose=FALSE)
+     preds <- predict(model, test.x, array.layout='rowmajor')
+     holdout_score <- rmse(preds, test.y)
+     list(Score=-holdout_score, Pred=-holdout_score)
+     # 「最大化」されてしまうので、最小化したい場合はメトリクスを負にする
+ }
> library(Metrics)
> library(rBayesianOptimization)
> opt_res <- BayesianOptimization(mxnet_holdout_bayes,
+                                 bounds=list(unit1=c(200L,500L),
+                                             unit2=c(10L,80L),
+                                             unit3=c(10L,50L),
+                                             num_r=c(80L,150L),
+                                             learn_r=c(1e-5,1e-3)),
+                                 init_points=20, n_iter=1, acq='ei',
+                                 kappa=2.576, eps=0.0, verbose=TRUE)
elapsed = 112.90	Round = 1	unit1 = 417.0000	unit2 = 59.0000	unit3 = 34.0000	num_r = 139.0000	learn_r = 0.0002	Value = -0.8849 
elapsed = 43.19	Round = 2	unit1 = 252.0000	unit2 = 63.0000	unit3 = 27.0000	num_r = 87.0000	learn_r = 0.0005	Value = -0.8860 
elapsed = 60.58	Round = 3	unit1 = 232.0000	unit2 = 19.0000	unit3 = 18.0000	num_r = 134.0000	learn_r = 0.0001	Value = -0.8905 
elapsed = 73.56	Round = 4	unit1 = 200.0000	unit2 = 50.0000	unit3 = 44.0000	num_r = 147.0000	learn_r = 0.0007	Value = -0.8769 
elapsed = 89.40	Round = 5	unit1 = 334.0000	unit2 = 27.0000	unit3 = 46.0000	num_r = 144.0000	learn_r = 0.0006	Value = -0.8825 
elapsed = 65.20	Round = 6	unit1 = 362.0000	unit2 = 11.0000	unit3 = 18.0000	num_r = 114.0000	learn_r = 0.0006	Value = -0.8828 
elapsed = 56.21	Round = 7	unit1 = 207.0000	unit2 = 18.0000	unit3 = 14.0000	num_r = 142.0000	learn_r = 0.0006	Value = -0.8815 
elapsed = 76.21	Round = 8	unit1 = 373.0000	unit2 = 16.0000	unit3 = 25.0000	num_r = 135.0000	learn_r = 0.0001	Value = -0.8958 
elapsed = 108.85	Round = 9	unit1 = 434.0000	unit2 = 79.0000	unit3 = 14.0000	num_r = 139.0000	learn_r = 0.0009	Value = -0.8798 
elapsed = 68.35	Round = 10	unit1 = 274.0000	unit2 = 46.0000	unit3 = 28.0000	num_r = 135.0000	learn_r = 0.0009	Value = -0.8803 
elapsed = 76.58	Round = 11	unit1 = 385.0000	unit2 = 12.0000	unit3 = 34.0000	num_r = 131.0000	learn_r = 0.0006	Value = -0.8846 
elapsed = 65.12	Round = 12	unit1 = 306.0000	unit2 = 80.0000	unit3 = 29.0000	num_r = 95.0000	learn_r = 0.0006	Value = -0.8818 
elapsed = 91.64	Round = 13	unit1 = 408.0000	unit2 = 59.0000	unit3 = 28.0000	num_r = 130.0000	learn_r = 0.0009	Value = -0.8803 
elapsed = 33.73	Round = 14	unit1 = 216.0000	unit2 = 30.0000	unit3 = 10.0000	num_r = 89.0000	learn_r = 0.0001	Value = -0.8921 
elapsed = 112.30	Round = 15	unit1 = 479.0000	unit2 = 57.0000	unit3 = 12.0000	num_r = 147.0000	learn_r = 0.0004	Value = -0.8849 
elapsed = 65.10	Round = 16	unit1 = 425.0000	unit2 = 14.0000	unit3 = 21.0000	num_r = 91.0000	learn_r = 0.0002	Value = -0.8858 
elapsed = 59.71	Round = 17	unit1 = 291.0000	unit2 = 33.0000	unit3 = 45.0000	num_r = 115.0000	learn_r = 0.0007	Value = -0.8837 
elapsed = 108.94	Round = 18	unit1 = 454.0000	unit2 = 52.0000	unit3 = 28.0000	num_r = 148.0000	learn_r = 0.0002	Value = -0.8871 
elapsed = 56.60	Round = 19	unit1 = 246.0000	unit2 = 71.0000	unit3 = 17.0000	num_r = 111.0000	learn_r = 0.0004	Value = -0.8847 
elapsed = 57.43	Round = 20	unit1 = 254.0000	unit2 = 18.0000	unit3 = 39.0000	num_r = 128.0000	learn_r = 0.0009	Value = -0.8788 
elapsed = 68.05	Round = 21	unit1 = 200.0000	unit2 = 61.0000	unit3 = 50.0000	num_r = 150.0000	learn_r = 0.0008	Value = -0.8828 

 Best Parameters Found: 
Round = 4	unit1 = 200.0000	unit2 = 50.0000	unit3 = 44.0000	num_r = 147.0000	learn_r = 0.0007	Value = -0.8769 

前の記事でやったランダムフォレストはおろかL1正則化線形回帰すら上回れませんでしたorz Deep Learningは難しい。。。多分パラメータの探索範囲が適切でない(もっとユニット数とエポック数の上限を増やす必要がある)のだろうとは思いますが、マシン負荷を考えるとあんまり大々的に増やすことも出来ないのでデモとしてはこんなものですよ、ということで。


ランダムサーチの方が良いのでは説


前回の記事でも引用したように、ベイズ最適化より普通にランダムサーチした方がパフォーマンスが良いのでは?というツッコミがあるようです。

個人的にはMCMCアルゴリズムの比較みたいな話に近いのかなという印象もありますが、現状ではMCMCとは違って機械学習の結果としての分類・回帰パフォーマンスというブラックボックスを陽に表せるわけではない(だからこそガウス過程で代用している)ので、例えばGibbs-samplingみたいなアイデアに行けるわけではないのでしょう。


もしかしたら(例えばの話ですが)最適化計画系で言えば量子アニーリングとか、MCMC系で言えばレプリカ交換とか、別のアイデアで精度を上げることはできるのかもしれませんが、これ以上は僕の理解の範囲を超えるのでやめておきます。。。