Think more, try less

Owen Zhangの言葉を大切にして,日々Kaggleに取り組んでいます.

Deep Learningライブラリ「MXNet」のR版をKaggle Otto Challengeで実践してみた

最近、Kagglerの間で話題となっていたDeep Learningのライブラリ「MXNet」のRパッケージ{mxnet}についてTJOさんよりご紹介がありました。 tjo.hatenablog.com

多くの方がこのブログで初めて{mxnet}を知った、もしくはあまり触れてこなかったが久々に使ってみようという段階かもしれませんが、私をはじめ多くのKagglerはすでに{mxnet}を多少なり使用していることもあり、下記のように解説記事の書くようにと無茶ぶりを頂きましたので、私からもご紹介させて頂きたいと思います。

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さんから提供されたサンプルコードを参考にしています。詳しくは下記リンクを参照してください。

Achieve 0.48 in 5 min with a deep net (feat. BatchNorm, PReLU) - Otto Group Product Classification Challenge | Kaggle

またMulti Layer PerceptronやReLU、Batch Normalizationの教科書的な解説は書籍に任せたいと思います。

f:id:Keiku:20160330184900p:plain

Customized mlogloss metricとEarly Stopping Ruleの設定

Otto Challengeで競う評価指標(Metric)は、multi-class logarithmic lossという指標であり、下記のように表されます。

{ \displaystyle
log loss = -\frac{1}{N}\sum_{i=1}^N\sum_{j=1}^My_{ij}\log(p_{ij})
}

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ファイルが生成されます。次の手順でモデリングを行っています。

  1. データを読み込み、読み込んだデータをscaleを使って標準化する。
  2. caret::createFoldsを使ってValidation datasetを用意する。
  3. Validation datasetに対して予測を行う。
  4. Validation datasetの結果、最適なnum.roundを見つけて、Test datasetに対して予測を行う。

結果は、Batch Normalizationありのモデルで0.50705、Batch Normalizationなしのモデルで0.47950となりました。

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を使っていることがほとんどですが、整理をすると以下のように考えています。

  1. Webサービスに組み込む、研究開発などの用途、またKaggleコンペティションにおいて魔改造が必要となる場合:Theano、Tensorflow、Caffe、Torch7、Chainer、MXNet(C++)など

  2. 通常のKaggleコンペティションやアドホックな分析用途:Keras、Lasagne、MXNet(Python/R)、Neonなど

この2番目の用途において、Rの{mxnet}によって分析可能となったことは多くのデータ分析者にとって喜ばしいことです。また、私としてはまだDeep Learningのモデリングは非常に難しく、さまざまなフレームワークでのDemoやExampleの実装を見て勉強する必要があると考えているため、Kerasに絞らずできるだけ多くのフレームワークに慣れ親しんでおこうと思っています。

追記 (2016/04/01)

TOP Kagglerのid:ultraistさんから貴重なコメント頂きました。