Windows 10 で GPU（CUDA）を利用するLightGBM のR-packageをセットアップする
基本的には公式サイトのガイドに従ってセットアップするだけ。だけなのだが、ちょこちょこ入れ忘れとか落とし穴があったりして、次から手間取らないように備忘録にしておく。
全体の準備
32bit 版のRをインストールしない
わりと嵌まったポイント。LightGBMは64bit版しかサポートしないが、32bit 版のRが入っているとビルドの際に32bit版をビルドしようとして失敗するとのことで、解決方法は、Rのインストール時に32bit版を入れないようにする（ホントにそう書いてある）。
Install CMake
公式サイト からダウンロードアイコンをクリック → インストーラーのダウンロード → セットアップ
Install Rtools
公式サイト からRtools35.exeをダウンロードしてインストール。
環境変数Pathに、Rtools\binとRtools\mingw_64\binを設定しておく。
Install Visual Studio
公式サイト からダウンロード。
インストールは、とりあえずC++だけでよい。
Install MS-API
MS-MPI SDK と setup.exe を両方公式サイトからダウンロード
msmpisdk.msi と msmpisetup.exe を実行してインストール。
Install git
公式サイト からダウンロード。
設定は適当に。
Install OpenCL
インストールガイドの
- For running on Intel, get Intel SDK for OpenCL.
- For running on AMD, get AMD APP SDK.
- For running on NVIDIA, get CUDA Toolkit.
に従って、今回はCUDA Toolkitをインストール。
今回は、特にカスタマイズしなくてもよかった。
Install Boost Binary
sourceforgeからダウンロード。
Visual Studio 2017 をインストールしたので、 msvc-14.1-64.exeをインストールする。
Boostの環境設定
BOOST_ROOTとBOOST_LIBRALYDIRのパスの設定を行う
環境変数の確認
「コントロール パネル > システムとセキュリティ > システム（Win+Pause/Break）」から「システムの詳細設定→環境変数」をチェック
以下の項目を確認：
- CMake関連のパスがPath変数に追加されている
- Rtools関連のパスがPath変数に追加されている
- MSAPI関連の環境変数がが登録されている
- openCL（CUDA）関連の環境変数が追加されている
- boost関連の環境変数が追加されている
Windows を再起動する
超重要。忘れる。
LightGBMのインストール
lgbdlを使ってインストールする。
lgbdlのインストール
devtools::install_github("Laurae2/lgbdl")
## Downloading GitHub repo Laurae2/lgbdl@master ## from URL https://api.github.com/repos/Laurae2/lgbdl/zipball/master ## Installing lgbdl
## "C:/PROGRA~1/R/R-35~1.1/bin/x64/R" --no-site-file --no-environ --no-save \ ## --no-restore --quiet CMD INSTALL \ ## "C:/Users/MyName/AppData/Local/Temp/Rtmp0SHyaL/devtools43c0670829ab/Laurae2-lgbdl-afff103" \ ## --library="C:/Users/MyName/Documents/R/win-library/3.5" \ ## --install-tests --no-multiarch
インストール
use_gpu = TRUEを指定する。R 3.5以降のバージョンではR35 = TRUEをセットするとご利益があるらしい（ビルドの際の最適化に関する警告が違う）。
lgbdl::lgb.dl(commit = "master", compiler = "vs", # Remove this for MinGW + GPU installation repo = "https://github.com/Microsoft/LightGBM", cores = 4, R35 = TRUE, use_gpu = TRUE)
## Installing package into 'C:/Users/MyName/Documents/R/win-library/3.5' ## (as 'lib' is unspecified)
インストール中のメッセージが大量に表示される。
## [1] TRUE
当初、LoadLibrary failure: %1 is not a valid Win32 application.というエラーが出て、インストールに失敗したので調べたら、Rのインストールの際に32bit版を入れてしまうと、両方でビルドしようとしてこける（64bit版しかサポートしない）とのこと。
このほか、ビルド中にどこかで止まる場合、チェックポイントのどれかが不十分と思われるので、再確認。
動作確認
library(lightgbm)
## Loading required package: R6
data(agaricus.train, package='lightgbm') train <- agaricus.train dtrain <- lgb.Dataset(train$data, label=train$label)
## Loading required package: Matrix
params <- list(objective="regression", metric="l2") model <- lgb.cv(params, dtrain, 10, nfold=5, min_data=1, learning_rate=1, early_stopping_rounds=10)
## [1]: valid's l2:1.43524e-16+9.92589e-17 ## [2]: valid's l2:6.30955e-32+5.75747e-32 ## [3]: valid's l2:7.50672e-35+1.33422e-34 ## [4]: valid's l2:7.03798e-35+1.35156e-34 ## [5]: valid's l2:6.93275e-35+1.35631e-34 ## [6]: valid's l2:6.91914e-35+1.35695e-34 ## [7]: valid's l2:6.91914e-35+1.35695e-34 ## [8]: valid's l2:6.91914e-35+1.35695e-34 ## [9]: valid's l2:6.91914e-35+1.35695e-34 ## [10]: valid's l2:6.91914e-35+1.35695e-34
（おまけ）install_github()を使って、非GPU版のインストール
GPU使わない場合はOpenCLとBoostの準備が不要なので、こちらの方が楽。
library(devtools)
options(devtools.install.args = "--no-multiarch") # if you have 64-bit R only, you can skip this
install_github("Microsoft/LightGBM", subdir = "R-package")
## Downloading GitHub repo Microsoft/LightGBM@master ## from URL https://api.github.com/repos/Microsoft/LightGBM/zipball/master ## Warning: GitHub repo contains submodules, may not function as expected! ## Installing lightgbm ## "C:/PROGRA~1/R/R-35~1.1/bin/x64/R" --no-site-file --no-environ --no-save \ ## --no-restore --quiet CMD INSTALL \ ## "C:/Users/MyName/AppData/Local/Temp/Rtmp0SHyaL/devtools43c039486000/Microsoft-LightGBM-aac8e8f/R-package" \ ## --library="C:/Users/MyName/Documents/R/win-library/3.5" \ ## --install-tests --no-multiarch ##
本当は、ビルド中にシステムメッセージが大量に表示される 。
autoxgboost を使ってみる
前回のTokyoRで＠hoxo-mさんがつぶやいていたautoxgboostを使ってみる。
xgboost の自動パラメータ調整は autoxgboost というのが便利そうだった。#tokyorhttps://t.co/LvwY9U2zyx
— hoxo_m (@hoxo_m) 2018年7月15日
何？
autoxgboostは、mlr と mlrMBOを使ってxgboostのモデルをベイズ最適化でチューニングするためのラッパー。気合い入れて検討する前にベースラインとる用途だけでも、かなり楽になるのではないかと思う。
(参考)
準備
インストールは github から。
install.packages("devtools") # if you have not installed "devtools" package devtools::install_github("ja-thomas/autoxgboost")
set.seed(0) data(spam, package = "kernlab")
70%を教師データとして使う。このデータのラベルはtype列。
train.set <- sample(1:NROW(spam), NROW(spam) * 0.7) table(spam$type[train.set]) #> #> nonspam spam #> 1954 1266
残りをテストデータに使う。
test.set <- setdiff(1:NROW(spam), train.set) table(spam$type[test.set]) #> #> nonspam spam #> 834 547
autoxgboost のロード。必要な周辺パッケージも併せてロードされる。
require(autoxgboost) #> Loading required package: autoxgboost #> Loading required package: ParamHelpers #> Warning: replacing previous import 'BBmisc::isFALSE' by #> 'backports::isFALSE' when loading 'ParamHelpers' #> Loading required package: mlr #> Loading required package: mlrMBO #> Loading required package: smoof #> Loading required package: BBmisc #> Loading required package: checkmate #> Loading required package: mlrCPO
モデルのチューニング
一度、教師データでチューニングしてから、テストデータで評価することにする。ほぼマニュアル通りに書き下し。MBOcontrol()の細かな設定は、後日改めて確認しておきたい。
tune.task <- makeClassifTask(data = spam[train.set, ], target = "type") ctrl <- makeMBOControl() ctrl <- setMBOControlTermination(ctrl, iters = 1L) #Speed up Tuning by only doing 1 iteration
チューニングの実行。これもほぼマニュアル通り。mlrの分類タスクは デフォルトの評価指標がthe mean misclassification error (mmce) なので、これをaucに変更する。このあたりはタスクにあった指標をお好みで選べばよい。
design.sizeとかtune.thresholdとか、まだ理解していない引数もあるので、これもおいおい確認。
res <- autoxgboost(tune.task, control = ctrl, measure = auc, nthread = 4, tune.threshold = FALSE) #> Computing y column(s) for design. Not provided. #> [mbo] 0: eta=0.0992; gamma=18.4; max_depth=19; colsample_bytree=0.521; colsample_bylevel=0.529; lambda=5.39; alpha=35.8; subsample=0.733; scale_pos_weight=0.00616 : y = 0.968 : 0.7 secs : initdesign #> [mbo] 0: eta=0.0715; gamma=0.109; max_depth=9; colsample_bytree=0.641; colsample_bylevel=0.952; lambda=1.46; alpha=480; subsample=0.865; scale_pos_weight=0.173 : y = 0.5 : 0.3 secs : initdesign #> [mbo] 0: eta=0.181; gamma=25.6; max_depth=13; colsample_bytree=0.928; colsample_bylevel=0.618; lambda=10.6; alpha=7.48; subsample=0.569; scale_pos_weight=0.606 : y = 0.97 : 0.8 secs : initdesign #> [mbo] 0: eta=0.138; gamma=0.0292; max_depth=16; colsample_bytree=0.728; colsample_bylevel=0.596; lambda=452; alpha=0.343; subsample=0.618; scale_pos_weight=0.00176 : y = 0.969 : 1.2 secs : initdesign #> [mbo] 0: eta=0.0208; gamma=2.48; max_depth=4; colsample_bytree=0.744; colsample_bylevel=0.752; lambda=0.0143; alpha=0.859; subsample=0.635; scale_pos_weight=0.0286 : y = 0.965 : 0.4 secs : initdesign #> [mbo] 0: eta=0.0567; gamma=36.3; max_depth=15; colsample_bytree=0.542; colsample_bylevel=0.769; lambda=30.5; alpha=153; subsample=0.979; scale_pos_weight=396 : y = 0.944 : 0.6 secs : initdesign #> [mbo] 0: eta=0.13; gamma=5.39; max_depth=17; colsample_bytree=0.785; colsample_bylevel=0.809; lambda=0.0852; alpha=0.192; subsample=0.672; scale_pos_weight=5.51 : y = 0.981 : 0.7 secs : initdesign #> [mbo] 0: eta=0.17; gamma=1.11; max_depth=8; colsample_bytree=0.569; colsample_bylevel=0.985; lambda=2.84; alpha=0.0632; subsample=0.785; scale_pos_weight=3.19 : y = 0.986 : 0.6 secs : initdesign #> [mbo] 0: eta=0.158; gamma=0.0599; max_depth=19; colsample_bytree=0.671; colsample_bylevel=0.693; lambda=275; alpha=0.0142; subsample=0.525; scale_pos_weight=131 : y = 0.954 : 0.3 secs : initdesign #> [mbo] 0: eta=0.0861; gamma=0.441; max_depth=5; colsample_bytree=0.818; colsample_bylevel=0.905; lambda=0.00104; alpha=0.00136; subsample=0.935; scale_pos_weight=0.0513 : y = 0.962 : 0.3 secs : initdesign #> [mbo] 0: eta=0.0414; gamma=0.0134; max_depth=3; colsample_bytree=0.999; colsample_bylevel=0.65; lambda=0.0032; alpha=19.7; subsample=0.914; scale_pos_weight=614 : y = 0.969 : 0.5 secs : initdesign #> [mbo] 0: eta=0.0284; gamma=0.0255; max_depth=11; colsample_bytree=0.943; colsample_bylevel=0.724; lambda=0.0355; alpha=1.83; subsample=0.547; scale_pos_weight=1.06 : y = 0.969 : 0.8 secs : initdesign #> [mbo] 0: eta=0.0812; gamma=0.185; max_depth=14; colsample_bytree=0.878; colsample_bylevel=0.565; lambda=0.115; alpha=0.0313; subsample=0.828; scale_pos_weight=0.0112 : y = 0.985 : 1.1 secs : initdesign #> [mbo] 0: eta=0.191; gamma=8.07; max_depth=7; colsample_bytree=0.851; colsample_bylevel=0.894; lambda=134; alpha=174; subsample=0.877; scale_pos_weight=19.5 : y = 0.945 : 0.3 secs : initdesign #> [mbo] 0: eta=0.12; gamma=0.542; max_depth=11; colsample_bytree=0.615; colsample_bylevel=0.846; lambda=0.534; alpha=0.00287; subsample=0.764; scale_pos_weight=28 : y = 0.981 : 0.7 secs : initdesign #> Loading required package: rgenoud #> ## rgenoud (Version 5.8-2.0, Build Date: 2018-04-03) #> ## See http://sekhon.berkeley.edu/rgenoud for additional documentation. #> ## Please cite software as: #> ## Walter Mebane, Jr. and Jasjeet S. Sekhon. 2011. #> ## ``Genetic Optimization Using Derivatives: The rgenoud package for R.'' #> ## Journal of Statistical Software, 42(11): 1-26. #> ## #> [mbo] 1: eta=0.125; gamma=0.0108; max_depth=5; colsample_bytree=0.978; colsample_bylevel=0.618; lambda=0.00409; alpha=0.137; subsample=0.799; scale_pos_weight=0.625 : y = 0.983 : 0.5 secs : infill_cb
チューニング後の結果を表示すると、おススメのパラメータセットを喋ってくれる。
print(res) #> Autoxgboost tuning result #> #> Recommended parameters: #> eta: 0.170 #> gamma: 1.111 #> max_depth: 8 #> colsample_bytree: 0.569 #> colsample_bylevel: 0.985 #> lambda: 2.836 #> alpha: 0.063 #> subsample: 0.785 #> scale_pos_weight: 3.189 #> nrounds: 25 #> #> #> Preprocessing pipeline: #> dropconst(rel.tol = 1e-08, abs.tol = 1e-08, ignore.na = FALSE) #> #> With tuning result: auc = 0.986
結果オブジェクトの構成は以下の通り。
str(res, 1) #> List of 5 #> $ optim.result :List of 9 #> ..- attr(*, "class")= chr [1:2] "MBOSingleObjResult" "MBOResult" #> $ final.learner :List of 11 #> ..- attr(*, "class")= chr [1:3] "CPOLearner" "BaseWrapper" "Learner" #> $ final.model :List of 8 #> ..- attr(*, "class")= chr [1:3] "CPOModel" "BaseWrapperModel" "WrappedModel" #> $ measure :List of 10 #> ..- attr(*, "class")= chr "Measure" #> $ preproc.pipeline:List of 24 #> ..- attr(*, "class")= chr [1:2] "CPOPrimitive" "CPO" #> - attr(*, "class")= chr "AutoxgbResult"
チューニングしたモデルの評価
ベストのlearnerとmodelも取っておいてくれるので、これを使う。
res$final.learner #> Learner classif.xgboost.custom.dropconst from package xgboost #> Type: classif #> Name: ; Short name: #> Class: CPOLearner #> Properties: numerics,missings,twoclass,multiclass,prob,featimp #> Predict-Type: prob #> Hyperparameters: nrounds=25,verbose=0,objective=binary:logistic,eta=0.17,gamma=1.11,max_depth=8,colsample_bytree=0.569,colsample_bylevel=0.985,lambda=2.84,alpha=0.0632,subsample=0.785,scale_pos_weight=3.19 res$final.model #> Model for learner.id=classif.xgboost.custom.dropconst; learner.class=CPOLearner #> Trained on: task.id = spam[train.set, ]; obs = 3220; features = 57 #> Hyperparameters: nrounds=25,verbose=0,objective=binary:logistic,eta=0.17,gamma=1.11,max_depth=8,colsample_bytree=0.569,colsample_bylevel=0.985,lambda=2.84,alpha=0.0632,subsample=0.785,scale_pos_weight=3.19
取り分けておいたテストデータを使って予測性能を評価する。
test.task <- makeClassifTask(data = spam[test.set, ], target = "type") pred <- predict(res$final.model, task = test.task) performance(pred, measures = auc) #> auc #> 0.9851731
このデータではテストセットの性能も、チューニングのときのものとあまり変わらない。
（修正あり）arulesの結果をdata.frameで探索する
R標準のデータフレーム＋dplyr等で、探索的にルール抽出⇔眺めるのにパッケージを作った。本家のarulesだとちょっとやりにくいなあと思っている人むけ。
使いかた
arules::inspect() の代わりに inspectDF() するだけ。
arules::DATAFRAME() でデータフレームとして出力できるとのこと。
抽出されたルールをdataframe化する関数としてはすでに arules::DATAFRAME() があると思いますが、新たに inspectDF を作られた理由ってありますでしょうか？
— データポエマー[,-5] (@bob3bob3) September 16, 2018
自作する前によく探すの、大事。
require(dplyr) require(arules) require(inspectDF) data(Groceries) params <- list(confidence=0.001, support=0.001, maxlen=7, minlen=2) glo.apriori <- apriori(Groceries, parameter = params, control = list(verbose=FALSE))
glo.inspectDF <- inspectDF(glo.apriori) glo.inspectDF %>% str #> 'data.frame': 40943 obs. of 8 variables: #> $ rule : chr "Rule 1" "Rule 2" "Rule 3" "Rule 4" ... #> $ LHS : chr "whole milk" "other vegetables" "other vegetables" "rolls/buns" ... #> $ RHS : chr "sparkling wine" "artif. sweetener" "bathroom cleaner" "nuts/prunes" ... #> $ n.lhs : num 1 1 1 1 1 1 1 1 1 1 ... #> $ support : num 0.00102 0.00102 0.00102 0.00102 0.00102 ... #> $ confidence: num 0.00398 0.00525 0.00525 0.00553 0.00583 ... #> $ lift : num 0.712 1.615 1.914 1.647 1.024 ... #> $ count : num 10 10 10 10 10 10 10 10 10 10 ...
取り出したルールは、グラフ描画用のユーティリティ関数で眺める。
set.seed(0) glo.inspectDF %>% arrange(support, confidence) %>% head(60) %>% plotRuleGraph()
インストール
GitHubから.
install.packages("devtools") # if you have not installed "devtools" package devtools::install_github("katokohaku/inspectDF")
ソースコードも同じ場所。
モチベーション
本家のsubset(), sort()とかでなく、通常のdata.frameの操作としてやりたい。
本家の inspect() がコンソールに全部出力するのがツライ
print(glo.apriori) #> set of 40943 rules rules.lhs %>% summary() #> rule LHS RHS n.lhs #> Length:8035 Length:8035 Length:8035 Min. :2.000 #> Class :character Class :character Class :character 1st Qu.:2.000 #> Mode :character Mode :character Mode :character Median :3.000 #> Mean :2.739 #> 3rd Qu.:3.000 #> Max. :5.000 #> support confidence lift count #> Min. :0.001017 Min. :0.01815 Min. : 0.4379 Min. : 10.00 #> 1st Qu.:0.001118 1st Qu.:0.14961 1st Qu.: 2.0673 1st Qu.: 11.00 #> Median :0.001322 Median :0.25175 Median : 2.6852 Median : 13.00 #> Mean :0.001740 Mean :0.30231 Mean : 2.9595 Mean : 17.11 #> 3rd Qu.:0.001830 3rd Qu.:0.41935 3rd Qu.: 3.5531 3rd Qu.: 18.00 #> Max. :0.022267 Max. :1.00000 Max. :22.1394 Max. :219.00
arulesは非常によくできたパッケージですが、高度にwrapされ過ぎていて、気軽に弄ろうと思うとツライ。本家のprint() も summary()も、探索的なルール分析にはちょっと不親切だと思いませんか。探索的なルール分析にはsubset(), sort()などが用意されているものの、表示はことごとく上記のそれです。ツライ。
また、本家の inspect() はprint()で返すため、コンソールにすべての出力が強制表示される。いつ終わるとも知れぬコンソール出力を眺めるの、ツライ。あと、標準では2番目の列名が無名なので、この後そのまま探索的に弄ろうとすると困る、というか矢印要る？ ツライ。
glo.inspect <- glo.apriori %>% head(10) %>% inspect() #> lhs rhs support confidence lift #> [1] {honey} => {whole milk} 0.001118454 0.733333333 2.870009 #> [2] {whole milk} => {honey} 0.001118454 0.004377238 2.870009 #> [3] {soap} => {whole milk} 0.001118454 0.423076923 1.655775 #> [4] {whole milk} => {soap} 0.001118454 0.004377238 1.655775 #> [5] {tidbits} => {soda} 0.001016777 0.434782609 2.493345 #> [6] {soda} => {tidbits} 0.001016777 0.005830904 2.493345 #> [7] {tidbits} => {rolls/buns} 0.001220132 0.521739130 2.836542 #> [8] {rolls/buns} => {tidbits} 0.001220132 0.006633499 2.836542 #> [9] {cocoa drinks} => {whole milk} 0.001321810 0.590909091 2.312611 #> [10] {whole milk} => {cocoa drinks} 0.001321810 0.005173100 2.312611 #> count #> [1] 11 #> [2] 11 #> [3] 11 #> [4] 11 #> [5] 10 #> [6] 10 #> [7] 12 #> [8] 12 #> [9] 13 #> [10] 13 glo.inspect %>% str #> 'data.frame': 10 obs. of 7 variables: #> $ lhs : Factor w/ 7 levels "{cocoa drinks}",..: 2 7 4 7 6 5 6 3 1 7 #> $ : Factor w/ 1 level "=>": 1 1 1 1 1 1 1 1 1 1 #> $ rhs : Factor w/ 7 levels "{cocoa drinks}",..: 7 2 7 4 5 6 3 6 7 1 #> $ support : num 0.00112 0.00112 0.00112 0.00112 0.00102 ... #> $ confidence: num 0.73333 0.00438 0.42308 0.00438 0.43478 ... #> $ lift : num 2.87 2.87 1.66 1.66 2.49 ... #> $ count : num 11 11 11 11 10 10 12 12 13 13
使い道
dplyrあたりと組み合わせて使うと、tidyに探索できるのではないか。
require(stringr) rules.lhs <- glo.inspectDF %>% filter(str_detect(LHS, pattern = "yogurt|sausage")) %>% arrange(confidence, lift) %>% filter(n.lhs > 1) rules.lhs %>% head() #> rule LHS RHS n support #> 1 Rule 97 whole milk,yogurt UHT-milk 2 0.001016777 #> 2 Rule 98 whole milk,yogurt red/blush wine 2 0.001016777 #> 3 Rule 99 whole milk,yogurt house keeping products 2 0.001016777 #> 4 Rule 100 whole milk,yogurt liver loaf 2 0.001016777 #> 5 Rule 7175 whole milk,yogurt chewing gum 2 0.001118454 #> 6 Rule 7176 whole milk,yogurt cling film/bags 2 0.001118454 #> confidence lift count #> 1 0.01814882 0.5425339 10 #> 2 0.01814882 0.9444108 10 #> 3 0.01814882 2.1767518 10 #> 4 0.01814882 3.5698730 10 #> 5 0.01996370 0.9485170 11 #> 6 0.01996370 1.7530626 11
プロットは hideaki さんのqiitaの記事をそのまま流用。エッジ上の円の面積はルールのsupport値、色の濃さはルールのconfidence値を表す。
adujust.support.sizeオプションでヒューリスティックなサイズの調整ができる。
require(stringr) rules.lhs <- glo.inspectDF %>% filter(str_detect(RHS, pattern = "yogurt")) %>% arrange(confidence, lift) %>% filter(n.lhs > 1) %>% head(15) par(mfrow = c(1,2)) set.seed(0) rules.lhs %>% plotRuleGraph(label = "default") set.seed(0) rules.lhs %>% plotRuleGraph(label = "adjusted rule size", adujust.support.size = 4000)
今後
{support, confidence, lift}と{size, color}の対応は可変にして、サイズも自動調整したい。
既存のGithubのレポジトリをRStudioでcloneする
RStudioの操作だけで完結したい人向け。タイトル通りのことしたくて戸惑ったので備忘録。
(1) File > New Project... > Version Control > Git
(2) 既存のレポジトリのアドレスをコピー
(3) コピーしたアドレスを(1)に張り付け
(4) clone完了
あとは、ローカルの変更をpushしたり、別の端末からpushされた変更をpullしたり。(1)からVersion Controlを選べばいいというのに気が付くまで、時間がかかった。
元ネタ？はこちら：
qiita.com
random forestを使った解釈性にもとづく推薦システム
[1706.06691] Interpretable Predictions of Tree-based Ensembles via Actionable Feature Tweakingという論文が提案されている。KDD2017でYahoo Researchの研究者によって発表されたもの、とのこと。Ensemble treesによって、あるサンプルがラベルAと予測されたとき、別のラベルBと予測されたかったら元のサンプルのどの特徴量にどのような操作を加えればよいか「改善案」を算出する。手元にあると便利そうだったのでRで実装したので紹介したい。
どんなことができるかの概要は、先日のTokyoRで喋ってきたので、そちらのスライドも参照して欲しい。
大雑把な手順は、下記のとおり：
- データのスケーリング
- 予測モデル（random forest）の構築
- ルール抽出
- 各事例への推薦の計算
なお、日本語の解説は、接点QBさんによる丁寧な記事がすでにあるので、そちらを参照されたい。
setten-qb.hatenablog.com
前準備
ライブラリとソースのロード。この記事で使用するコードはgithubに置いてあるので、適当にローカルにおいて読み込む。そのうち時間を作ってパッケージ化したい。
github.com
require(tidyverse) require(magrittr) require(randomForest) devtools::install_github("hoxo-m/pforeach") require(pforeach) source("./R/tweak_feature.R")
pforeachに隅々まで依存しているので、別途インストールが必須。
訓練セットとテストセットの準備
例によってirisを使いたいところだが、今回はspamデータセットを使ってspamと分類されるサンプルをnonspamに分類されるように推薦を行う。irisは多クラスだしルールも最小限なんで良いんだけど、アヤメの花びらを切って貼ってもカキツバタにはならないので今回はスキップ。実用上、変更できない説明変数が予測モデルに含まれる場合には、禁止リストをかませるとか適切な層別化するといった工夫が別途必要になる。
# data preparation ------------------------------------------------------------- set.seed(777) data(spam, package = "kernlab") dataset <- sample_frac(spam) n.test <- floor(NROW(dataset) *0.1) dataset.train <- chop(dataset, n.test) dataset.test <- tail(dataset, n.test) dim(dataset);dim(dataset.train);dim(dataset.test)
chop()はデータを後ろから指定したサイズだけ切り落とすための自作ユーティリティ。
random forestの学習とモデル選択
訓練セットで学習して重要度のプロットと弱学習器の数の評価を行う。この辺りはルーチンなので変わったことはしていないが、ある程度予測精度が担保されていることが前提なので、ここで予測精度がダメな場合は再検討。
# bilding randomForest ----------------------------------------- X <- dataset.train[, 1:(ncol(dataset.train)-1)] true.y <- dataset.train[, ncol(dataset.train)] forest.all <- randomForest(X, true.y, ntree=500) forest.all par(mfrow=c(1,2)) varImpPlot(forest.all) # to view varImp, x3 & x4 should be removed. plot(forest.all) par(mfrow=c(1,1))
この方法は総当たり探索的なアプローチなので、変数（ルール）の数とツリーの数が増えると計算量が爆発的に増える。したがって変数選択とツリーの絞り込みのほか、max_depthを減らす工夫が重要。
今回は上位12変数と100ツリーくらいあれば良さそうなので、適当に切り落とす。
# model shrinkage based on importance ------------------------- top.importance <- forest.all$importance %>% data.frame %>% tibble::rownames_to_column(var = "var") %>% arrange(desc(MeanDecreaseGini)) %>% head(12) dataset.train.fs <- dataset.train %>% select(top.importance$var) dataset.test.fs <- dataset.test %>% select(top.importance$var)
データのスケーリングと縮約済みモデルの学習
この手法はハイパーパラメータを一つだけにするため、各変数を標準化する必要がある。実用上、テストセットは訓練セットのスケーリングパラメータで標準化される必要があるのと、推薦後の改善は標準化前の尺度に復元する必要があり、毎回attributeを取ってくるのは面倒なので、それぞれrescale()とdescale()として、ユーティリティ関数を用意した。
# scaling feature-selected data --------------------------- X.train <- scale( dataset.train.fs ) X.test <- rescale( dataset.test.fs, scaled = X.train ) dataset.test.fs[1:6, 1:6] descale(X.test, scaled = X.train)[1:6, 1:6] descale(X.test, scaled = X.test)[1:6, 1:6]
ちゃんと復元されていることを確認して欲しい。スケーリングした訓練データでモデルの学習を行う。
forest.rf <- randomForest(X.train, true.y, ntree=100) forest.all forest.rf plot(forest.rf)
変数を減らしたので予測精度が下がっているが、「予測ミス（擬陽性）を起こすサンプルは本質的には改善されて良いのでは」という立場をとって許容範囲と思うことにして目をつぶることにする。
feature tweaking
ここからが提案手法の処理。手順は、random forest予測モデルを学習したのち、（１）ルールを抽出する＋変更ルールを算出する（２）個別のサンプルの予測と推薦案を算出する（３）標準化の復元、の1+3ステップとなる。
# （０）random forest予測モデルを学習する forest.rf <- randomForest(X.train, true.y, ntree=100) # （１）ルールを抽出する＋変更ルールを算出する es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE) # （２）個別のサンプルの予測と推薦案を算出する tweaked <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = TRUE) # （３）標準化の復元 dt <- descale.tweakedFeature(tweaked, X.test)
（１）ルールを抽出する＋変更ルールを算出する
random forestではgetTree()で弱学習器を取ってこれる。これをつかって各ラベルに分類されるルールセットは、ラベルごとにdata.frameのlistとして整理しなおして抽出する。
# test sampling (for demo) ep <- getRules.randomForest(forest.rf, k=2, label.to = NULL) ep %>% str(2) ep[["spam"]]$path[[1]]
xgboostにも同様にルールセットを取り出すxgb.model.dt.tree()があるので、そのうち対応したい。
抽出されたルールセットと変更ルール（候補案）のセットは対応しているため、抽出と同時に計算してしまう。元論文のアルゴリズムは各サンプルの評価時にこれを毎回計算する構造になっており、素直に実装してしまうと計算量が大変になるので、今から実装しようと思う人は注意。
epsiron は変更の強さを指定するハイパーパラメタ。原理上、どうやっても予測を覆せないサンプルがあり、この量は変更可能な提案が可能なサンプル数に影響する。
es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3)
ntreeはルールセットを抽出する弱学習器の数で、省略するとすべてのツリーを使用する。デフォルトでは先頭からntreeで指定したツリーを解析対象にするが、データ不均衡に対処するために訓練セットをアンダーサンプリングして学習した幾つかのモデルをrandomForest::combine(...)で結合したい場合がある。この場合は、先頭から順番に解析してしまうと困るので、resampleオプションをTRUEに指定することでランダムに弱学習器を評価することができる。
es.rf <- set.eSatisfactory(forest.rf, ntree = 30, epsiron = 0.3, resample = TRUE)
（２）個別のサンプルの予測と推薦案を算出する
予測対象のデータセットが label.from と予測されたとき、label.to と予測されるための改善案を計算する。学習データと同じスケールで標準されていないとマトモな推薦が出来ないので注意。
tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = FALSE)
各サンプルについてルールを全部評価するので結構な計算時間がかかるため、並列化オプション.doparはデフォルトではTRUEになっている。
tweaked_ <- tweak(es.rf, newdata= X.test, label.from = "spam", label.to = "nonspam", .dopar = TRUE)
（３）標準化の復元
元のサンプルの操作方向とその量が知りたいのでスケーリングを解除する。元のサンプルの特徴量(original)、予測を変えるための変更後の特徴量(suggest)、予測の変更に必要な努力量=両者の差分(diff)のlistを返す。
>|r|
dt <- descale.tweakedFeature(tweaked, X.test)
dt %>% str(1)
|
個別のサンプルへの推薦パターンの可視化
個別の推薦結果を可視化する。操作量を相対的に比較したいため、標準化されたままのものを可視化している。サンプルによって違う推薦パターンが見えるので面白い。
# plot suggestion for each instance which(tweaked$predict == "spam") plot.suggest(tweaked, 4) plot.suggest(tweaked, 11)
操作量の大きい順に、0でない変数だけ可視化したければ、それぞれ .ordered = TRUEと、.nonzero.only = TRUEを指定する。
plot.suggest(tweaked, 15) plot.suggest(tweaked, 15, .ordered = TRUE, .nonzero.only = TRUE)
変数重要度の可視化
論文中ではユーザーによる納得感（アンケート調査?）をもとに評価しているが、それだと機械的には評価できないので、「予測サンプル集団において操作回数が多かった変数は重要なんじゃないか」という考えと、「操作された変数の操作の方向性が揃っている変数は、その意味合いを解釈できそう」と考えて、これらを各変数の重要度と思うことにして、それぞれ可視化する。
# Plot population importances pp <- plot.tweakedPopulation(tweaked, type = "absoluteSum") pp <- plot.tweakedPopulation(tweaked, type = "direction") pp <- plot.tweakedPopulation(tweaked, type = "frequency")
type = "frequency"は、サンプル集団において各変数が操作をされた集計値を可視化する。
type = "direction"は、各サンプルの推薦で各変数がどのような操作をされたかを集計して可視化する。相対的な重要性を比較したいので、標準化された状態の推薦量を使用する。
変更量の総和を可視化する"absoluteSum"オプションも作ってみたが、これはあまり有意義と感じなかった。
Lassoの非ゼロ変数を代替する候補を可視化する
最近、モデルの解釈性（interpretability）について色々と研究が出ており、興味もあってぼんやりと追いかけている。
“Finding Alternate Features in Lassos”という論文が提案されている：
Satoshi Hara and Takanori Maehara (2016): "Finding Alternate Features in Lasso", in NIPS 2016 workshop on Interpretable Machine Learning for Complex Systems, December 10th, 2016. （arXiv, github）
著者の日本語ページにもコンパクトな解説が掲載されている。
Lassoには変数間の相関が高いとき、それらの変数の中から一つだけを選択したがる性質がある(「見落としの問題」と呼ぶそうだ）。この問題に対して、この論文では、LASSOによって選択された変数がもし使われなかったとしたら、かわりにどの変数が選択されるかと考える。著者らによると、代替となる変数を列挙提示することで、納得感や結果の解釈性が得られるのではないか、と期待している。顧客に選択肢を残しておくというのは良いんじゃないだろうか。
具体的な手続きはシンプルで、Lassoで選択された係数の推定量が非ゼロの変数
と推定量ゼロの変数
の組を考え、
のもとで全ての組
について
、
とおいて目的関数を解きなおす。
と再推定されたら、「
は
の代替となる変数である」と提示され、やっぱり
なら「
は
の置き換えにはならない」とする。
また、代替となる変数同士の優劣は、目的関数の変化量を指標スコアとして相対的に順位づけされる。
で、ざっと論文を読んだ後、著者のPythonコードを見たらRでも簡単にできそうだったので、シンプルなOLS + soft-thresholdな罰則項という最小限だけをＲで実装した(github)。
ローカルに持ってきてDESCRIPTIONを用意してビルドすればパッケージになる（はず）。二分グラフを縦にプロットする方法で地味に手こずったが、この論文の可視化以外にもどこかで役に立つだろう。
面白いなと思った点は、目的関数をタスクに合わせてフレキシブルに設計できるところ。実際に、論文では損失関数をloglossに差し替えて二値分類タスクに適用・実験している。正則化項のほうを弄れば、例えばadaptive Lassoなどへの拡張もできるのではないかと思う。
また、各変数の相関行列と代替性の指標スコアとを比較すると、必ずしも大きさの順番がが一致しないのも興味深い。相関の高い変数が上位に来るのはわかりやすい一方、罰則の強さ(）の選び方によっては、割と相関の高くない変数でも代替候補として挙げられてくる。これはもしかしたら、前処理としてマルチコを除く作業をこの手法によってモデル構築後の二次工程に回せるかもしれないし、データによっては発見的な議論につながるかもしれない。
の選び方について、スライドの例ではcross-validationで求めた値を使った。しかし、変数が増えていくと見なければいけない情報が多すぎるので、強めの罰則をかけてエッセンシャルな変数だけに絞って分析するほうがシンプルかな、とは思った。この辺りは好みもあるかもしれないが、元の論文中では推奨値に関する議論は特に言及がなかったと思う。
代替候補の変数は必ずオリジナルの変数よりも目的関数の評価が下がることもあり、おそらく予測はオリジナルののモデルでやったほうがよい。そういう理由でprediction()は実装しなかった。この方法で予測したいとか2値分類について同様の分析がしたい方は本家のPython版を利用されたい。
推定値がズレる縮小推定の問題や、サンプルサイズよりも変数の数が大きいときに見過ごすといった、Lassoの癖というか性質に対処するアプローチの一例については、以前紹介したことがあるのでそちらも参照されたい。
mxnetパッケージをUbuntuで使うためのセットアップ
Ubuntu14.04でRからmxnet使いたくなってセットアップ。ガイド見ながらやっていけば終わり...と思ったらinstall_depsでいくつか引っかかった。Windows/Macばりに甘やかされたい。
次からコピペでやりたいので自分用にメモしておく。
http://mxnet.readthedocs.io/en/latest/how_to/build.html#r-package-installationmxnet.readthedocs.io
# install mxnet [~]$ sudo apt-get install build-essential git libatlas-base-dev libopencv-dev [~]$ git clone --recursive https://github.com/dmlc/mxnet [~]$ cd mxnet [~/mxnet]$ make
今回、makeは特に失敗せず完了。makeがおわったら以下を実行する。インストールガイドどおりにコマンドラインからやっても良い。
require(devtools); require(methods); require(httr) setwd("mxnet/R-package/") set_config(use_proxy(url="プロキシ使ってれば指定",port=8080)) options(repos=c(CRAN='https://cran.rstudio.com')); install_deps(dependencies = TRUE)
install_depsが動作するディレクトリに注意。今回はプロキシ指定とかが(なぜか)上手くいかなかったので、コマンドラインじゃなくスクリプトであれこれ試行錯誤した。
いろいろ入れてる途中でimagerでfftw3.hがないよと言われて失敗したので追加。
[~]$ sudo apt-get install fftw3 fft3-dev install_depsが沈黙したら下記を実行。 >|sh| [~]$ cd mxnet [~/mxnet]$ make rpkg
あとはインストール。
[~]$ R CMD INSTALL mxnet_0.7.tar.gz
以上をシェルスクリプトにまとめれば終わり。
