最近、Kagglerの間で話題となっていたDeep Learningのライブラリ「MXNet」のRパッケージ{mxnet}についてTJOさんよりご紹介がありました。
tjo.hatenablog.com
多くの方がこのブログで初めて{mxnet}を知った、もしくはあまり触れてこなかったが久々に使ってみようという段階かもしれませんが、私をはじめ多くのKagglerはすでに{mxnet}を多少なり使用していることもあり、下記のように解説記事の書くようにと無茶ぶりを頂きましたので、私からもご紹介させて頂きたいと思います。
@Keiku いえいえ、Keikuさんの記事楽しみにしております！（Kaggle masterはデータ分析者の憧れですよ〜）
— TJO (@TJO_datasci) 2016年3月29日
TJOさんの解説記事では、Deep Learningの代表的な事例でもある「Convolutional Neural Network」について解説されていますが、多くのデータサイエンスに関わる方、皆がコンピュータビジョン（画像認識等）を業務としているわけではないと思います。逆にコンピュータビジョンを業務とされている方々は、最近ではPythonから利用できるChainer、KerasといったDeep Learningライブラリをすでに実践しているというのが現状です。そこでより多くの方が気軽にRの{mxnet}を利用できるようにKaggleでの事例をご紹介したいと思います。
Kaggle「Otto Group Product Classification Challenge」
今回は、Otto Challengeで実践することにしました。 www.kaggle.com
このコンペティションは、Otto Group（20カ国以上で子会社を持つ世界でも最大級のECカンパニーのひとつ）から提供されたデータセットを題材としており、20万以上の商品を特徴付けた93変数（これら変数名は匿名化された数値変数）とその商品についての9種類の商品カテゴリ（この目的変数も匿名化されており9種類のカテゴリカル変数）のラベルを持ったデータセットとなっています。本コンペ参加者は、未知のラベルのデータセットに対してこの9種類の商品カテゴリを予測し、その精度を競いました。
このOtto Challengeは過去のKaggleの歴史上、最も多くの参加者を集めたコンペティションであり、Kaggleの中でも代表的なコンペのひとつです。マーケテイングにおける商品カテゴリ分類はセグメンテーションの下地にもなる重要なマーケテイング課題のひとつであると思います。その問題設定の取っ付きやすさもあり多くの参加者を集めた理由のひとつと考えられます。
また、最近では、Kaggle等のコンペティションにてGBDTのライブラリであるXGBoostが猛威をふるっているとたびたび話題にあがりますが、Kagglerから正確に現状を伝えると、XGBoostとDeep Learningのアンサンブルが猛威をふるっているというのが今の状況です。そしてこのOtto Challengeは、XGBoostだけでなく、Deep Learningの威力をまざまざと見せつけられた代表的な事例であると言えるでしょう。このコンペティションの上位者は絶対と言っていいほど、XGBoostとDeep Learningのアンサンブルを行っています。XGBoostでほとんどの参加者が良い精度を出すのは当たり前となっているため、ゆえにDeep Learningがいかにできるかが勝敗を分ける鍵であったと思います。
今回は、その勝負の鍵であったDeep LearningをRの{mxnet}をもちいて実践してみることとしました。
Multi Layer Perceptron
今回の問題設定は、9種類の商品カテゴリ分類となるため、下記の図のようなニューラルネットワークを組むことが基本となります。入力層はデータセット93変数を入力とし、隠れ層ではユニットを512、活性化関数としてRectified Linear Unit (ReLU)、さらにBatch NormalizationとDropoutを設定しています。この隠れ層をさらに同じアーキテクチャーの層を追加し、最終的に出力層としてSoftmax関数で9つの出力をしています。
このアーキテクチャーは、Kerasの開発者で現在Googleにお勤めのfcholletさんから提供されたサンプルコードを参考にしています。詳しくは下記リンクを参照してください。
またMulti Layer PerceptronやReLU、Batch Normalizationの教科書的な解説は書籍に任せたいと思います。
Customized mlogloss metricとEarly Stopping Ruleの設定
Otto Challengeで競う評価指標(Metric)は、multi-class logarithmic lossという指標であり、下記のように表されます。
Rの{mxnet}では評価指標としてAccuracyやRMSEは標準で実装されていますが、このmulti-class loglossという指標は、標準でサポートされていなかったため、mx.metric.customを用いて自作する必要があります。以下のようになります。
# Customized mlogloss eval_metric mx.metric.mlogloss <- mx.metric.custom("mlogloss", function(label, pred){ label_mat <- data.frame(i = 1:length(label), j = label + 1, x = rep(1, length(label))) label_mat <- sparseMatrix(i = label_mat$i, j = label_mat$j, x = label_mat$x) label_mat <- as.matrix(label_mat) return(MultiLogLoss(label_mat, t(pred))) })
また、XGBoostには標準でサポートされているEarly Stopping Ruleの設定は{mxnet}ではmx.callback functionを使って自分で実装する必要があります。以下のようになります。
# Customized mx.callback function for early stopping mx.callback.early.stop <- function (period, nrounds, logger = NULL, maximize=TRUE, verbose=FALSE) { function(iteration, nbatch, env, verbose) { if (nbatch%%period == 0 && !is.null(env$metric)) { result <- env$metric$get(env$train.metric) if (nbatch != 0) cat(paste0("Batch [", nbatch, "] Train-", result$name, "=", result$value, "\n")) if (!is.null(logger)) { if (class(logger) != "mx.metric.logger") { stop("Invalid mx.metric.logger.") } logger$train <- c(logger$train, result$value) if (!is.null(env$eval.metric)) { result <- env$metric$get(env$eval.metric) if (nbatch != 0) cat(paste0("Batch [", nbatch, "] Validation-", result$name, "=", result$value, "\n")) logger$eval <- c(logger$eval, result$value) if (iteration > nrounds) { if(maximize) { if(sum(result$value > tail(logger$eval,nrounds)) == 0 ) { bst.round <- which.max(logger$eval) if (verbose) { print(paste0("Best round: ", bst.round, "; Best metric: ", max(logger$eval))) print(paste0("This round: ", result$value)) } return(FALSE) } } else if(sum(result$value < tail(logger$eval,nrounds)) == 0 ) { bst.round <- which.min(logger$eval) if (verbose) { print(paste0("Best round: ", bst.round, "; Best metric: ", min(logger$eval))) print(paste0("This round: ", result$value)) } return(FALSE) } } } } return(TRUE) } } }
実際にKaggleにSubmitしてみた
以下のコードを実行することで、submissionファイルが生成されます。
options(scipen = 100) options(dplyr.width = Inf) options(dplyr.print_max = Inf) require(mxnet) require(caret) require(dplyr) require(data.table) require(MLmetrics) # to evaluate multi-class logarithmic loss require(stringr) require(Matrix) setwd("kaggle/otto-group-product-classification-challenge") train <- fread("data/train.csv", data.table = F) test <- fread("data/test.csv", data.table = F) sampleSubmission <- fread("data/sampleSubmission.csv", data.table = F) y <- as.numeric(as.factor(train[["target"]])) - 1 # drop id and target X <- train[, -c(1, ncol(train))] dim(X) # [1] 61878 93 X_test <- test[, -c(1)] dim(X_test) # [1] 144368 93 All <- bind_rows(X, X_test) # Scaling All <- All %>% dplyr::mutate_each(funs(scale)) X <- data.matrix(All[1:nrow(X),]) X_test <- data.matrix(All[(nrow(X)+1):nrow(All),]) # ----------------------------------------- # Evaluate validation dataset # ----------------------------------------- # Extract validation dataset set.seed(123) cv_list <- createFolds(y, k = 5, list = TRUE, returnTrain = FALSE) X_train <- X[-cv_list[[1]], ] dim(X_train) # [1] 49501 93 X_valid <- X[cv_list[[1]], ] dim(X_valid) # [1] 12377 93 y_train <- y[-cv_list[[1]]] y_valid <- y[cv_list[[1]]] # Customized mlogloss eval_metric mx.metric.mlogloss <- mx.metric.custom("mlogloss", function(label, pred){ label_mat <- data.frame(i = 1:length(label), j = label + 1, x = rep(1, length(label))) label_mat <- sparseMatrix(i = label_mat$i, j = label_mat$j, x = label_mat$x) label_mat <- as.matrix(label_mat) return(MultiLogLoss(label_mat, t(pred))) }) # Customized mx.callback function for early stopping mx.callback.early.stop <- function (period, nrounds, logger = NULL, maximize=TRUE, verbose=FALSE) { function(iteration, nbatch, env, verbose) { if (nbatch%%period == 0 && !is.null(env$metric)) { result <- env$metric$get(env$train.metric) if (nbatch != 0) cat(paste0("Batch [", nbatch, "] Train-", result$name, "=", result$value, "\n")) if (!is.null(logger)) { if (class(logger) != "mx.metric.logger") { stop("Invalid mx.metric.logger.") } logger$train <- c(logger$train, result$value) if (!is.null(env$eval.metric)) { result <- env$metric$get(env$eval.metric) if (nbatch != 0) cat(paste0("Batch [", nbatch, "] Validation-", result$name, "=", result$value, "\n")) logger$eval <- c(logger$eval, result$value) if (iteration > nrounds) { if(maximize) { if(sum(result$value > tail(logger$eval,nrounds)) == 0 ) { bst.round <- which.max(logger$eval) if (verbose) { print(paste0("Best round: ", bst.round, "; Best metric: ", max(logger$eval))) print(paste0("This round: ", result$value)) } return(FALSE) } } else if(sum(result$value < tail(logger$eval,nrounds)) == 0 ) { bst.round <- which.min(logger$eval) if (verbose) { print(paste0("Best round: ", bst.round, "; Best metric: ", min(logger$eval))) print(paste0("This round: ", result$value)) } return(FALSE) } } } } return(TRUE) } } } # 2 layers MLP (Multi layer perceptron) data <- mx.symbol.Variable("data") fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=512) act1 <- mx.symbol.Activation(fc1, name="relu1", act_type="relu") act1 <- mx.symbol.BatchNorm(act1, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn1") dp1 <- mx.symbol.Dropout(act1, p = 0.5, name="dp1") fc2 <- mx.symbol.FullyConnected(dp1, name="fc2", num_hidden=512) act2 <- mx.symbol.Activation(fc2, name="relu2", act_type="relu") act2 <- mx.symbol.BatchNorm(act2, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn2") dp2 <- mx.symbol.Dropout(act2, p = 0.5, name="dp2") fc3 <- mx.symbol.FullyConnected(dp2, name="fc3", num_hidden=9) softmax <- mx.symbol.SoftmaxOutput(fc3, name="sm") devices <- mx.cpu() mx.set.seed(123) logger <- mx.metric.logger$new() model <- mx.model.FeedForward.create( softmax, X=X_train, y=y_train, eval.data=list(data=X_valid, label=y_valid), ctx=mx.cpu(), num.round=100, array.batch.size=128, optimizer = "sgd", learning.rate=0.01, momentum=0.9, initializer=mx.init.uniform(0.01), eval.metric=mx.metric.mlogloss, array.layout="rowmajor", epoch.end.callback = mx.callback.early.stop(1, 5, logger=logger, maximize=FALSE, verbose = TRUE) ) # With Batch Normalization # [1] "Best round: 19; Best metric: 0.499036037381019" # [1] "This round: 0.505431101369955" # Without Batch Normalization # [1] "Best round: 51; Best metric: 0.480038855711653" # [1] "This round: 0.481640159668552" # Check network architecture graph.viz(model$symbol$as.json()) # Predict validation dataset preds = predict(model, X_valid) preds = t(preds) y_valid_mat <- data.frame(i = 1:length(y_valid), j = y_valid + 1, x = rep(1, length(y_valid))) y_valid_mat <- sparseMatrix(i = y_valid_mat$i, j = y_valid_mat$j, x = y_valid_mat$x) y_valid_mat <- as.matrix(y_valid_mat) # Evaluate multi-class logarithmic loss print(MultiLogLoss(y_true=y_valid_mat, y_pred=preds)) # With Batch Normalization # [1] 0.5047899 # Without Batch Normalization # [1] 0.4805768 # ----------------------------------------- # Predict test dataset # ----------------------------------------- # 2 layers MLP (Multi layer perceptron) data <- mx.symbol.Variable("data") fc1 <- mx.symbol.FullyConnected(data, name="fc1", num_hidden=512) act1 <- mx.symbol.Activation(fc1, name="relu1", act_type="relu") act1 <- mx.symbol.BatchNorm(act1, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn1") dp1 <- mx.symbol.Dropout(act1, p = 0.5, name="dp1") fc2 <- mx.symbol.FullyConnected(dp1, name="fc2", num_hidden=512) act2 <- mx.symbol.Activation(fc2, name="relu2", act_type="relu") act2 <- mx.symbol.BatchNorm(act2, eps=1e-06, momentum=0.9, fix.gamma=FALSE, name="bn2") dp2 <- mx.symbol.Dropout(act2, p = 0.5, name="dp2") fc3 <- mx.symbol.FullyConnected(dp2, name="fc3", num_hidden=9) softmax <- mx.symbol.SoftmaxOutput(fc3, name="sm") devices <- mx.cpu() mx.set.seed(123) model <- mx.model.FeedForward.create( softmax, X=X, y=y, ctx=mx.cpu(), num.round=19, array.batch.size=128, optimizer = "sgd", learning.rate=0.01, momentum=0.9, initializer=mx.init.uniform(0.01), eval.metric=mx.metric.mlogloss ) # Check network architecture graph.viz(model$symbol$as.json()) # Predict test dataset preds = predict(model, X_test) preds = t(preds) preds = data.frame(preds) names(preds) <- str_c("Class_", 1:9) # Make a submission submission <- bind_cols(data_frame(id = sampleSubmission$id), preds) write.csv(submission, "submit/submission.csv", quote = F, row.names = F) # Private Leaderboard Score # With Batch Normalization # 0.50705 # Without Batch Normalization # 0.47950
モデリングに関するメモ
Batch Normalizationの有無での比較
上記のコードで、Batch Normalizationありなしで比較も行ってみましたが、Batch Normalizationを入れた場合の方が精度を落としてしまいました。元の変数が頻度を表すような変数であり、バッチごとに正規化されたことが悪さをしたのか？実のところよくわかっていません。
Test dataset評価時におけるnum.roundの設定
上記のコードでは、Early Stopping RuleによりValidation datasetにおいて最も精度の高いnum.roundを見つけてそれをTest datasetでの評価時に設定していますが、Validation datasetは80%のデータで作成したモデルであるため、Test datasetを予測するためのTraining dataset(100%)とは少し異なります。num.roundは少し多めに設定しても問題ないと思いますし、最終的には適当にsubmitしてLeaderboard Scoreが最も高くなるnum.roundにすべきとしか言いようがないですが、このあたりどのように設定することがいいのか常に悩むところになります。
今後どのDeep Learningフレームワークを使うべきか？
Rの{mxnet}についてほんの少しご紹介しましたが、実際これだけDeep Learningフレームワークがあるなかでどれを選べばいいのか困ると思います。私のTwitterのタイムラインでは多くのKagglerが「Keras」をあげていました。私も同じくKerasを使っていることがほとんどですが、整理をすると以下のように考えています。
Webサービスに組み込む、研究開発などの用途、またKaggleコンペティションにおいて魔改造が必要となる場合：Theano、Tensorflow、Caffe、Torch7、Chainer、MXNet(C++)など
通常のKaggleコンペティションやアドホックな分析用途：Keras、Lasagne、MXNet(Python/R)、Neonなど
この2番目の用途において、Rの{mxnet}によって分析可能となったことは多くのデータ分析者にとって喜ばしいことです。また、私としてはまだDeep Learningのモデリングは非常に難しく、さまざまなフレームワークでのDemoやExampleの実装を見て勉強する必要があると考えているため、Kerasに絞らずできるだけ多くのフレームワークに慣れ親しんでおこうと思っています。