機械学習のパラメータチューニングというと大なり小なり大変な部分があって、今年のエイプリルフール記事に皆さん引っかかって下さったところを見るにパラメータチューニングを簡単に済ませたい！と願う人々は世の中多いようです（笑）。
少し前のMXnetを使った記事でも取り上げましたが、そのパラメータチューニングを迅速に済ませようというアイデアの一つがベイズ最適化(Bayesian Optimization)です。
要は、グリッドサーチのように網羅的に最適なパラメータを探索しに行くのではなく、一つのパラメータで精度をチェックしたらその次は精度が上がりやすそうな方向にベイズ的に逐次改善を行いながら探索していく、という方法のことです。
世の中色々seminar paper的なものがあるようですが、arXivから@enakai00さんが見つけてきて下さったのがこれ。
日本語のブログ記事なら、id:olanleedさんのこちらの記事が良いと思います。
ということでこれらのreferencesを紹介してしまえば理論的なポイントは全ておしまいなのですが、リンク先を読むのも面倒という人のためにごくごく簡単な説明だけを書いておきます。ただし僕は数学が大の苦手なので数式は紹介する気もない上に、そもそも上記のreferencesからの受け売り以外何もありません（笑）。
ベイズ最適化とは
上記のarXiv論文のこのFig.1が最も分かりやすいかと思います（ちょうど最適化したい関数も1次元でここで想定している「機械学習のCV精度」の最適化というテーマとマッチするので）。
端的に言えば「機械学習のCV精度をパラメータを変えながら改善していくという試行錯誤のプロセスを、ベイズ統計の考え方を活用して効率化したもの」です。特にここで挙げられている方法は「機械学習のCV精度がある未知のガウス過程に従う」と仮定した上で、これに対して事前に決め打ちで選んだ「
の期待値を算出してくれる関数
」を用意し、交互に次の
次の
と反復繰り返し更新していくことで、
を最大にする
を探し出す、というアルゴリズムです。
その様子が上の図に分かりやすく書かれています。まず初期値を出したいのでランダムに2点選んできて黒線で表されるの値が2つ求まったところから始まります（上段）。これに対して
が緑線で表される関数としてこのステップで求まります。そこで
の中をある程度ランダムに探索した中で最も高い値を与えた点（赤い▼のところ）を次のサンプリング点として、
の値を求めます（中段）。この結果を受けて
がさらに更新されてそのまた次の赤い▼が決まるので、これをそのまた次のサンプリング点として
の値を求めます（下段）。これを繰り返すことで、徐々に
を最大化する
へと近付いていく、というわけです。
以上の手続きを擬似コードにしたのが論文中のAlgorithm 1です。
こう書くと簡単そうに見えますが、色々な理屈が裏側にはあります。まずガウス過程に共分散を与えるためにカーネルを事前に決めておく必要があります（細かい話は「黄色い本」ことPRMLの6.4節を参照のこと）。次にについては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
ランダムフォレストで試してみる
> 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系で言えばレプリカ交換とか、別のアイデアで精度を上げることはできるのかもしれませんが、これ以上は僕の理解の範囲を超えるのでやめておきます。。。
