tuneRF関数の挙動の検証
RのrandomForestパッケージのtuneRF関数が遅いというお話があったので，調べてみた．
tuneRF関数は，ランダムフォレストを用いて予測モデルを構築する際に使用する特徴量の個数を貪欲的な方法により求める．
tuneRF関数だけでなく，Rで機械学習のアルゴリズムをチューニングする方法については，id:TJOさんの記事「Rで機械学習するならチューニングもグリッドサーチ関数orオプションでお手軽に」が詳しいので，是非参照してほしい．
検証に使用するデータと環境
UCI Machine Learning RepositoryのBank Marketingデータセットを使用する．このデータセットは，ポルトガルの銀行で電話でダイレクトマーケティングを実施した際に収集したデータを用いて，予測モデルのアルゴリズムに投入できるように特徴量を構築したもの．
まずは，データを取得して解凍する．
$ wget http://archive.ics.uci.edu/ml/machine-learning-databases/00222/bank.zip $ unzip bank.zip
同じディレクトリに"bank-full.csv", "bank-names.txt", "bank.csv"という3つのファイルが解凍される．"bank-full.csv"がすべてのレコードを含んでおり，ここではこのデータを使用する．
データをRに読み込んで，行頭やサイズ，各列のデータ型を確認する．
> bank <- read.csv("bank-full.csv", sep=";") > head(bank, 3) age job marital education default balance housing loan contact day 1 58 management married tertiary no 2143 yes no unknown 5 2 44 technician single secondary no 29 yes no unknown 5 3 33 entrepreneur married secondary no 2 yes yes unknown 5 month duration campaign pdays previous poutcome y 1 may 261 1 -1 0 unknown no 2 may 151 1 -1 0 unknown no 3 may 76 1 -1 0 unknown no > dim(bank) [1] 45211 17 > sapply(bank, class) age job marital education default balance housing loan "integer" "factor" "factor" "factor" "factor" "integer" "factor" "factor" contact day month duration campaign pdays previous poutcome "factor" "integer" "factor" "integer" "integer" "integer" "integer" "factor" y "factor"
行数が45,211，列数が17であることを確認できる．最後の列(y)は，ダイレクトマーケティングに顧客が反応したかどうかを表すカテゴリ変数であり，この項目を目的変数として使用する(Yes: 反応した，No: 反応していない)．残りの16個の列を説明変数と使用する．上記の結果を見ると，16個のうち9項目がカテゴリ変数であることを確認できる．
なお，以下で検証に使用する環境は次の通り．
まずは単純にtuneRF関数を実行してみる
まずは特別な設定を行わずにtuneRF関数を実行してみよう．tuneRF関数の第1引数には説明変数，第2引数には目的変数を指定する．また，doBest引数をTRUEに指定すると，評価が最も良いモデルを返すようになる．
> set.seed(123) > system.time(tuneRF(bank %>% select(-y), bank$y, doBest=TRUE)) mtry = 4 OOB error = 9.46% Searching left ... mtry = 2 OOB error = 9.78% -0.03295162 0.05 Searching right ... mtry = 8 OOB error = 9.48% -0.001402197 0.05 ユーザ システム 経過 77.131 1.725 81.135
上記を見ると，おおよそ1分20秒程度実行に要したことを確認できる．この結果，特徴量の個数が
- 4個のときに，Out-of-Bag誤差(OOB error)は9.46%
- 2個のときに，Out-of-Bag誤差は9.78%
- 8個のときに，Out-of-Bag誤差は9.48%
となり，特徴量の個数が4個のときにOut-of-Bag誤差が最少となり，この個数に設定するのが良さそうであることがわかる*1．
構築する決定木の個数を増やしてみる
ntreeTry引数はデフォルトでは50となっており，50個の決定木を構築することがわかる．500個の決定木を構築するように指定してみよう．
> set.seed(123) > system.time(tuneRF(bank %>% select(-y), bank$y, ntreeTry=500, doBest=TRUE)) mtry = 4 OOB error = 9.1% Searching left ... mtry = 2 OOB error = 9.53% -0.04689018 0.05 Searching right ... mtry = 8 OOB error = 9.18% -0.008503401 0.05 ユーザ システム 経過 237.656 2.943 244.725
構築する決定木の個数が10倍になったので，実行時間も約10倍になるのかと思うところであるが，3倍程度に抑えられていることを確認できる．
tuneRF関数の引数と仕様
tuneRF関数にはいくつかの引数を指定できる．
> args(tuneRF) function (x, y, mtryStart = if (is.factor(y)) floor(sqrt(ncol(x))) else floor(ncol(x)/3), ntreeTry = 50, stepFactor = 2, improve = 0.05, trace = TRUE, plot = TRUE, doBest = FALSE, ...) NULL
| 引数 | 説明 | データ型 |
| x | 説明変数 | data.frame |
| y | 目的変数 | factor or numeric |
| mtryStart | 選択する特徴量の個数の初期値 | integer |
| ntreeTry | 構築する決定木の個数 | integer |
| stepFactor | 特徴量の個数を減少または増加させるファクター | integer |
| improve | モデルの評価指標の改善率 | numeric |
| trace | チューニングの過程をトレースして画面に表示するかどうかを指定するフラグ | logical |
| doBest | 評価が最も良いモデルを返すかどうかを指定 | logical |
また，tuneRF関数自体もそれほど長くないので，中身を読んでみる．
tuneRF <- function(x, y, mtryStart=if(is.factor(y)) floor(sqrt(ncol(x))) else floor(ncol(x)/3), ntreeTry=50, stepFactor=2, improve=0.05, trace=TRUE, plot=TRUE, doBest=FALSE, ...) { if (improve < 0) stop ("improve must be non-negative.") # yが因子かどうかを判定 classRF <- is.factor(y) errorOld <- if (classRF) { # yが因子の場合はクラス分類を実行し，誤分類率により評価 randomForest(x, y, mtry=mtryStart, ntree=ntreeTry, keep.forest=FALSE, ...)$err.rate[ntreeTry,1] } else { # yが因子でない場合は回帰を実行し，平均二乗平均誤差により評価 randomForest(x, y, mtry=mtryStart, ntree=ntreeTry, keep.forest=FALSE, ...)$mse[ntreeTry] } if (errorOld < 0) stop("Initial setting gave 0 error and no room for improvement.") # traceがTRUEに指定された場合は，チューニング過程を画面に出力 if (trace) { cat("mtry =", mtryStart, " OOB error =", if (classRF) paste(100*round(errorOld, 4), "%", sep="") else errorOld, "\n") } oobError <- list() oobError[[1]] <- errorOld names(oobError)[1] <- mtryStart # 以下，特徴量の個数を増減させて予測モデルを構築し結果を評価(left: 個数を削減，right: 個数を増加) for (direction in c("left", "right")) { if (trace) cat("Searching", direction, "...\n") Improve <- 1.1*improve mtryBest <- mtryStart mtryCur <- mtryStart while (Improve >= improve) { mtryOld <- mtryCur mtryCur <- if (direction == "left") { # 個数の削減は，現在の個数をstepFactor引数で割ったものを新しい個数とする max(1, ceiling(mtryCur / stepFactor)) } else { # 個数の増加は，現在の個数にstepFactor引数を掛けたものを新しい個数とする min(ncol(x), floor(mtryCur * stepFactor)) } if (mtryCur == mtryOld) break errorCur <- if (classRF) { # yが因子の場合はクラス分類を実行し，誤分類率により評価 randomForest(x, y, mtry=mtryCur, ntree=ntreeTry, keep.forest=FALSE, ...)$err.rate[ntreeTry,"OOB"] } else { # yが因子でない場合は回帰を実行し，平均二乗平均誤差により評価 randomForest(x, y, mtry=mtryCur, ntree=ntreeTry, keep.forest=FALSE, ...)$mse[ntreeTry] } if (trace) { cat("mtry =",mtryCur, "\tOOB error =", if (classRF) paste(100*round(errorCur, 4), "%", sep="") else errorCur, "\n") } oobError[[as.character(mtryCur)]] <- errorCur # 予測精度(誤差)の改善率の算出 Improve <- 1 - errorCur/errorOld cat(Improve, improve, "\n") if (Improve > improve) { errorOld <- errorCur mtryBest <- mtryCur } } } mtry <- sort(as.numeric(names(oobError))) res <- unlist(oobError[as.character(mtry)]) res <- cbind(mtry=mtry, OOBError=res) # plotがTRUEに指定された場合は，横軸を特徴量の個数，縦軸を評価値とする折れ線グラフをプロット if (plot) { plot(res, xlab=expression(m[try]), ylab="OOB Error", type="o", log="x", xaxt="n") axis(1, at=res[,"mtry"]) } # doBestがTRUEに指定された場合は，評価が最も良いモデルを返り値とする if (doBest) res <- randomForest(x, y, mtry=res[which.min(res[,2]), 1], ...) res }
以上のtuneRF関数の実装を読むと，特徴量の個数と評価結果を保持するノードを木構造で表現していくことがわかる．
そのため，この処理の高速化の方針として，
- 木の左側の計算と右側の計算を並列化する
- 個々のrandomForest関数で決定木を構築する処理を並列化する
の2つが考えられる．これ以外の方法は，randomForest関数，およびその内部のC++やFortranのコードを改変することになる．
なお，前者は，木が深い時に特に有効であり，後者は木の深さに関係なく，ある程度の個数の決定木を構築するならば有効な手段だと思われる．
*1:ランダムフォレストはOut-of-Bag誤差が求まるのでクロスバリデーション不要とも言われるが，実際はクロスバリデーションなども組み合わせて，もう少し精緻に検証することが望ましいと思われる．
Juliaによる機械学習の予測モデル構築・評価
これは，Julia Advent Calendar 2014 14日目の記事です．MLBaseパッケージを用いて機械学習の予測モデルを構築し，評価する方法について説明します．
以下では，Julia0.3.2，MLBase0.5.1，DecisionTree0.3.4，RDatasets0.1.1を使用しています．
Juliaで使用できる機械学習の手法
Juliaで使用できる機械学習の手法には，以下のようなものがある．
| 手法 | パッケージ |
|---|---|
| 決定木 | DecisionTree |
| ランダムフォレスト | DecisionTree, RandomForests(by @bicycle1885さん) |
| サポートベクタマシン | SVM, LIBSVM |
他の手法については，Awesome Machine Learningにまとまっている．
ランダムフォレストを試してみる
Juliaでランダムフォレストを実行するためには，DecisionTreeパッケージ，またはRandomForestsパッケージ(by @bicycle1885さん)を用いる．ここでは，DecisionTreeパッケージによる実行方法を示す．
Pima Indianデータセットに対して，ランダムフォレストの予測モデルを構築し，テストデータに対する予測結果を評価している．
using RDatasets using DecisionTree srand(123) # 訓練データの作成 Pima_tr = dataset("MASS", "Pima.tr") X_tr = array(Pima_tr[:, 1:7]) y_tr = array(Pima_tr[:, 8]) # テストデータの作成 Pima_te = dataset("MASS", "Pima.te") X_te = array(Pima_te[:, 1:7]) y_te = array(Pima_te[:, 8]) # 予測モデルの構築(特徴量の個数:4, 木の個数:2000, 訓練データの割合:70%) model = build_forest(y_tr, X_tr, 4, 2000, 0.7) # テストデータに対する予測 pred = apply_forest(model, array(Pima_te[:, 1:7])) # 分割表 confusion_matrix(y_te[:, 1], pred[:, 1])
以上を実行すると，下記の結果を得る．
Classes: {"No","Yes"}
Matrix:
2x2 Array{Int64,2}:
195 28
46 63
Accuracy: 0.7771084337349398
Kappa: 0.4723594347321851
他にも，DecisionTreeパッケージにはクロスバリデーションを実行するnfoldCV_forest関数が用意されている．次の例は，Pima Indianデータセットに対して，10-foldのクロスバリデーションを実行している．
# 10-foldクロスバリデーションの実行(特徴量の個数:4, 木の個数:200, fold数:10, 訓練データの割合:70%) accuracy = nfoldCV_forest(labels, features, 4, 200, 10, 0.7)
以上のように，DecisionTreeパッケージを使用すると，ランダムフォレストの予測モデルを構築し評価することが可能である．しかし，予測モデルに用いるアルゴリズムが異なると，別のパッケージを使用し異なる関数を用いて上記のタスクを実行する必要がある．アルゴリズムによらず，統一的なフレームワークにより予測モデルを構築し評価することができれば非常に便利である．
MLBaseパッケージによる予測モデル構築・評価の統一的なフレームワーク
MLBaseは，Rのcaretやmlr，Pythonのscikit-learnのように，機械学習の予測モデルを構築・評価する汎用的なパッケージである．このパッケージを使用することにより，クロスバリデーション，ハイパーパラメータのグリッドサーチなどの予測モデルの構築・評価に関わるタスクを統一的なフレームワークで実行することが可能になる．
クロスバリデーション
MLBaseでは，クロスバリデーションはcross_validate関数を用いることにより実行できる．予測モデルを構築する関数，予測モデルを評価する関数の2つを作成して，cross_validate関数に渡す必要がある．
以下の例では，Pima Indianデータセットに対して，10-foldのクロスバリデーションを実行している．
using MLBase srand(123) # 訓練データの作成 Pima_tr = dataset("MASS", "Pima.tr") X_tr = array(Pima_tr[:, 1:7]) y_tr = array(Pima_tr[:, 8]) # テストデータの作成 Pima_te = dataset("MASS", "Pima.te") X_te = array(Pima_te[:, 1:7]) y_te = array(Pima_te[:, 8]) # 予測モデルを構築する関数 function trainfun(inds) # ランダムフォレストによる予測モデルの構築(特徴量の個数:4, 木の個数:2000, 訓練に用いるデータの割合:70%) model = build_forest(y_tr[inds, 1], X_tr[inds, :], 4, 2000, 0.7) return model end # 予測モデルを評価する関数 function evalfun(model, inds) # 予測 pred = apply_forest(model, X_tr[inds, :]) # 分割表 conf_mat = confusion_matrix(y_tr[inds, 1], pred).matrix # Precisionの算出 prec = conf_mat[1, 1]/sum(conf_mat[:, 1]) return prec end # サンプル数 const n = size(X_tr, 1) # クロスバリデーションの実行 scores = cross_validate( inds -> trainfun(inds), (model, inds) -> evalfun(model, inds), n, Kfold(n, 10))
以上のコードを実行すると，下記の結果が表示される．
10-element Array{Float64,1}: 0.785714 0.705882 0.8 0.615385 0.692308 0.733333 0.642857 0.857143 0.846154 0.947368
cross_validate関数の第4引数に，各foldごとに用いるデータのインデクスを指定する．以上では，10-foldのクロスバリデーションによりデータを分割した．クロスバリデーションのスキームとしては，以下の5つの手法
- K-foldクロスバリーデション(Kfold関数)
- 層別K-foldクロスバリデーション(StratifiedKfold関数)
- 一つ抜きクロスバリデーション(LOOCV関数)
- ランダムサブサンプリング(RandomSub関数)
- 層別ランダムサブサンプリング(StratifiedRandomSum関数)
の5つが提供されている．
ハイパーパラメータのグリッドサーチ
MLBaseでハイパーパラメータのグリッドサーチを実行するためには，gridtune関数を使用する．
using MLBase srand(123) # 訓練データの作成 Pima_tr = dataset("MASS", "Pima.tr") X_tr = array(Pima_tr[:, 1:7]) y_tr = array(Pima_tr[:, 8]) # テストデータの作成 Pima_te = dataset("MASS", "Pima.te") X_te = array(Pima_te[:, 1:7]) y_te = array(Pima_te[:, 8]) # 予測モデルを構築する関数 function estfun(ntree, nfeature) # 予測モデルの構築(訓練データの割合:70%) model = build_forest(y_tr, X_tr, nfeature, ntree, 0.7) return model end # 予測精度を評価する関数 function evalfun(model) # テストデータに対する予測 pred = apply_forest(model, X_te) # 分割表 conf_mat = confusion_matrix(y_te[:, 1], pred[:, 1]).matrix # 評価指標(Precision)の算出 prec = conf_mat[1, 1]/sum(conf_mat[:, 1]) return prec end # グリッドサーチによるハイパーパラメータの評価(木の個数:{5,10,15}, 特徴量の個数:{4,5,6}) r = gridtune(estfun, evalfun, ("ntree", [5, 10, 15]), ("nfeature", [4, 5, 6]); verbose=true)
以上の例では，ランダムフォレストの木の個数を5,10,15個，特徴量の個数を4,5,6個として，これらの2つのハイパーパラメータのグリッドサーチを実行している．この結果，下記が表示される．
[ntree=5, nfeature=4] => 0.8034934497816594 [ntree=10, nfeature=4] => 0.8035714285714286 [ntree=15, nfeature=4] => 0.8122270742358079 [ntree=5, nfeature=5] => 0.8301886792452831 [ntree=10, nfeature=5] => 0.8285714285714286 [ntree=15, nfeature=5] => 0.7773279352226721 [ntree=5, nfeature=6] => 0.8109243697478992 [ntree=10, nfeature=6] => 0.8302752293577982 [ntree=15, nfeature=6] => 0.8060344827586207 (Ensemble of Decision Trees Trees: 10 Avg Leaves: 20.6 Avg Depth: 8.5,(10,6),0.8302752293577982)
この結果，木の個数を10個，特徴量の個数を6個としたときに，Precisionは0.830275と最大値となり，この値にハイパーパラメータを設定するのが良さそうであることが分かる．
ハイパーパラメータのグリッドサーチとクロスバリデーションを組み合わせて実行する関数は提供されていないようだが，以上で説明した事項を組み合わせると実行できる．
mlrパッケージによる予測モデルの構築・評価
これは，R Advent Calendar 2014 6日目の記事です．
本記事では，Rで機械学習の予測モデルの構築・評価を統一的なフレームワークで実行するmlrパッケージについて入門的な説明を行います．
mlrパッケージとは
mlrパッケージは，Michael Lang氏によって開発されており，機械学習の予測モデルを構築し評価する統一的なフレームワークを提供する．use!R 2014でも発表が行われている．
基本的な流れ
kernlabパッケージに収録されているspamデータセットに対して，メールがスパムかどうかを判定するための学習器を作成してみよう．ここでは，カーネルSVMを使用して10-foldのクロスバリデーションを実行してみる．
> install.packages("mlr") > library(mlr) > install.packages("kernlab") > library(kernlab) > data(spam) > # spamデータセットに対して目的変数をtypeとしたクラス分類タスクの作成 > task <- makeClassifTask(id="spam", data=spam, target="type") > # カーネルSVMの学習器の作成 > learner <- makeLearner("classif.ksvm") > # 10-foldクロスバリデーションを実行するためのデータセットの作成 > rdesc <- makeResampleDesc(method="CV", iters=10) > # 学習の実行 > r <- resample(learner=learner, task=task, rdesc)
以下，各機能について概要的な説明を行う．
実行可能なタスク
下表のように，5つのタスクが提供されている．
| 関数名 | タスク |
|---|---|
| makeClassifTask | クラス分類 |
| makeClusterTask | クラスタリング |
| makeCostSensTask | コスト考慮型学習 |
| makeRegrTask | 回帰 |
| makeSurvTask | 生存時間分析 |
リサンプリング方法
リサンプリング方法は，以下の6つが利用できる．カッコ内は，makeResampleDesc関数のmethod引数に指定する文字列を表している．
- クロスバリデーション("CV")
- 1つ抜きクロスバリデーション("LOOCV")
- 繰り返しクロスバリデーション("RepCV")
- Out-of-bagブートストラップ("Bootstrap")
- サプサンプリング("Subsample")
- ホールドアウト検定("Holdout")
ハイパーパラメータの最適化
ハイパーパラメータの最適化は，makeParamSet関数により探索するパラメータの範囲を指定することにより実行する．
> library(mlr) > library(kernlab) > data(spam) > task <- makeClassifTask(id="spam", data=spam, target="type") > learner <- makeLearner("classif.ksvm") > rdesc <- makeResampleDesc(method="CV", iters=10) > par.set <- makeParamSet( + makeDiscreteParam("C", values=2^(-2:2)), + makeDiscreteParam("sigma", values=2^(-2:2)) > ) > ctrl <- makeTuneControlGrid() > res <- tuneParams("classif.ksvm", task=task, resampling=rdesc, par.set=par.set, control=ctrl)
属性選択
一般的に，属性選択はフィルタ法，ラッパ法の2つに大別される．フィルタ法は，学習器を構築することなく各属性の重要度を定量化するアプローチである．一方で，ラッパ法は，学習器を構築しながら書く属性の重要度を定量化する．
mlrパッケージでは，ラッパ法については以下のように4つの属性選択のアルゴリズムが提供されている．
| 関数名 | タスク |
|---|---|
| makeFeatSelControlExhaustive | Exhaustive search |
| makeFeatSelControlGA | 遺伝的アルゴリズム |
| makeFeatSelControlRandom | ランダム探索 |
| makeFeatSelControlSequential | 決定論的な前方/後方探索 |
並列計算
クロスバリデーションやグリッドサーチによるハイパーパラメータの最適化などの処理は，並列計算が可能である．
Javaで分散処理
Javaで分散処理する必要が生じたので、調査のメモ。今回は、以下のページを参考にCORBA+RMIで分散処理をしてみることにする。OSはUbuntu14.04。
上記のリンク先の説明を読むと、Java SE 5以前は、クライアント側にスタブが、サーバ側にスケルトンが必要だったが、Java SE 5以降は不要になってJavaプログラムだけが必要な状況のようだ。
JBossToolsのインストール
JBoss EAPをダウンロードして、bin直下のrun.shを実行しようとすると、以下のエラーが出る。
$ ./run.sh ./run.sh: 3: ./run.sh: Bad substitution ======================================================================================== To start a JBoss Enterprise Application Platform 6 Standalone Server, a single server instance, use the command: /home/sfchaos/Downloads/jboss-eap-6.3/bin/standalone.sh To start a JBoss Enterprise Application Platform 6 Managed Domain, allowing control and management of multiple instances, use the command: /home/sfchaos/Downloads/jboss-eap-6.3/bin/domain.sh
そこで、ここではstandaloneモードで起動しよう。
$ ./standalone.sh
分散アプリケーションの実装
サーバ側は、次のように接続方法の定義し、その接続方法を用いた待ち受け状態にし、クライアントからのメソッド呼び出しの動作について設定する。
ServerClient.java
package remoting.server.test; import org.jboss.remoting.InvokerLocator; import org.jboss.remoting.transport.Connector; public class SimpleServer { public static void main(String[] args) throws Exception { InvokerLocator myLocator = new InvokerLocator("socket://127.0.0.1:8080"); Connector connector = new Connector(); connector.setInvokerLocator(myLocator.getLocatorURI()); connector.start(); connector.addInvocationHandler("simpleSystem",new SimpleServerInvocationHandler()); try { while (true) { Thread.sleep(1000); } } catch (Exception e) { e.printStackTrace(); } } }
SimpleServerInvocator.java
package remoting.server.test; import javax.management.MBeanServer; import org.jboss.remoting.InvocationRequest; import org.jboss.remoting.ServerInvocationHandler; import org.jboss.remoting.ServerInvoker; import org.jboss.remoting.callback.InvokerCallbackHandler; /** * メソッド呼び出しのためのハンドラ */ public class SimpleServerInvocationHandler { /** * 受信したメッセージを処理し、メッセージを返却します。 */ public Object invoke(InvocationRequest invocation) { System.out.println(" request：「" + invocation.getParameter() + "」"); return "ServerからClientへ"; } public void addListener(InvokerCallbackHandler arg0) { } public void removeListener(InvokerCallbackHandler arg0) { } public void setInvoker(ServerInvoker arg0) { } public void setMBeanServer(MBeanServer arg0) { } }
一方で、クライアント側は次のように、リモート呼び出しが実行されるinvokeメソッドを実装して、呼び出し元へ処理結果を返す。
SimpleClient.java
package remoting.server.test; import javax.management.MBeanServer; import org.jboss.remoting.InvocationRequest; import org.jboss.remoting.ServerInvocationHandler; import org.jboss.remoting.ServerInvoker; import org.jboss.remoting.callback.InvokerCallbackHandler; /** * メソッド呼び出しのためのハンドラ */ public class SimpleServerInvocationHandler { /** * 受信したメッセージを処理し、メッセージを返却します。 */ public Object invoke(InvocationRequest invocation) { System.out.println(" request：「" + invocation.getParameter() + "」"); return "ServerからClientへ"; } public void addListener(InvokerCallbackHandler arg0) { } public void removeListener(InvokerCallbackHandler arg0) { } public void setInvoker(ServerInvoker arg0) { } public void setMBeanServer(MBeanServer arg0) { } }
あとは、このプログラムを実行すればO.K.
doParallel関数に直接ワーカープロセス数を指定するとゾンビプロセスが残る件
次から次へと迫り来る原稿の嵐に追われている休みの昼下がり、何気なく目を向けたTLにこんなつぶやきが。
並列化した残骸のRScrpt.exeどうやって処分したらいいの…#メモリを圧迫し続けています
— Hadleyに憑依されてるテラモナギ (@teramonagi) 2014, 9月 22というわけで調べてみました。
状況の再現(Ubuntu)
まずは、Ubuntu-14.04での再現。
> library(foreach) # doParallelを読みこめばforeachも読み込まれるが念のため > library(doParallel) > registerDoParallel(4) > foreach (i=1:32) %dopar% sqrt(i) > system("ps") PID TTY TIME CMD 19531 pts/14 00:00:00 bash 19766 pts/14 00:00:00 R 19771 pts/14 00:00:00 R <defunct> 19772 pts/14 00:00:00 R <defunct> 19773 pts/14 00:00:00 R <defunct> 19774 pts/14 00:00:00 sh 19775 pts/14 00:00:00 ps > stopImplicitCluster() # クラスタを終了する > system("ps") PID TTY TIME CMD 19531 pts/14 00:00:00 bash 19766 pts/14 00:00:00 R 19771 pts/14 00:00:00 R <defunct> 19772 pts/14 00:00:00 R <defunct> 19773 pts/14 00:00:00 R <defunct> 19776 pts/14 00:00:00 sh 19777 pts/14 00:00:00 ps
ゾンビプロセスが残っていることが確認できる。このRのマスタープロセスを終了しない限り、ゾンビプロセスが残ることを確認。
続いて、
> library(foreach) # doParallelを読みこめばforeachも読み込まれるが念のため > library(doParallel) > cl <- makeCluster(4) > registerDoParallel(cl) > foreach (i=1:32) %dopar% sqrt(i) > stopCluster(cl) > system("ps") PID TTY TIME CMD 19531 pts/14 00:00:00 bash 19782 pts/14 00:00:00 R 19785 pts/14 00:00:00 R 19794 pts/14 00:00:00 R 19803 pts/14 00:00:00 R 19812 pts/14 00:00:00 R 19820 pts/14 00:00:00 sh 19821 pts/14 00:00:00 ps > stopCluster(cl) > system("ps") PID TTY TIME CMD 19531 pts/14 00:00:00 bash 19782 pts/14 00:00:00 R 19830 pts/14 00:00:00 sh 19831 pts/14 00:00:00 ps
無事にワーカープロセスを終了できていることが確認される。
状況の再現(Windows)
続いて、Windows。
> library(foreach) # doParallelを読みこめばforeachも読み込まれるが念のため > library(doParallel) > registerDoParallel(4) > foreach(i=1:32) %dopar% sqrt(i) > stopImplicitCluster()
ワーカープロセスは消えなかった。タスクマネージャ等は省略。
> library(foreach) # doParallelを読みこめばforeachも読み込まれるが念のため > library(doParallel) > cl <- makeCluster(4) > registerDoParallel(cl) > foreach (i=1:32) %dopar% sqrt(i) > stopCluster(cl)
どうしてこうなるのか
doParallelパッケージのregisterDoParallel関数で、上記の現象が関連する部分は23-45行目の以下の箇所。
23 if (missing(cl) || is.numeric(cl)) { 24 if (.Platform$OS.type == "windows") { 25 if (!missing(cl) && is.numeric(cl)) { 26 cl <- makeCluster(cl) 27 } 28 else { 29 if (!missing(cores) && is.numeric(cores)) { 30 cl <- makeCluster(cores) 31 } 32 else { 33 cl <- makeCluster(3) 34 } 35 } 36 assign(".revoDoParCluster", cl, pos = .options) 37 setDoPar(doParallelSNOW, cl, snowinfo) 38 } 39 else { 40 if (!missing(cl) && is.numeric(cl)) { 41 cores <- cl 42 } 43 setDoPar(doParallelMC, cores, mcinfo) 44 } 45 }
- 26行目でmakeCluster関数によりワーカープロセスを生成
- 36行目で、doParallelパッケージの.optionsオブジェクトが確保しているメモリアドレスに、revoDoParClusterオブジェクトにワーカープロセスの情報を保持したオブジェクトclを付値する。
という流れでワーカープロセスが生成され、情報が保持されていることが分かる。
一方で、Windows以外のOSだと、43行目でdoParallelMC関数を呼び出しているが、その中で同じような処理を行っている。
さて、こうして生成されたワーカープロセスは、並列計算が終了したら停止しなければならない。これは、doParallelパッケージのstopImplicitCluster関数を用いて行えるはずである。しかし、上記の再現結果を見れば分かるようにそうにはなっていない。そこで、stopImplicitCluster関数の実装を見てみると、次のようになっている。
1 if (exists(".revoDoParCluster", where = .options) && !is.null(.revoDoParCluster)) { 2 stopCluster(.revoDoParCluster) 3 remove(".revoDoParCluster", where = .options) 4 }
1行目が問題。if文の判定で.revoDoParClusterオブジェクトがNULLかどうかについて確かめているが、それが間違い。これは、本当は次のように.optionsのメモリアドレスに格納された.revoDoParClusterオブジェクトを取得して、そのオブジェクトがNULLかどうかを調べないといけない。さらに、それに呼応して、.optionsのメモリアドレスに格納されたオブジェクトを取得して、ワーカープロセスを停止する処理になっていなければならない。
1 if (exists(".revoDoParCluster", where = .options) && !is.null(get(".revoDoParCluster", envir = .options))) { 2 stopCluster(get(".revoDoParCluster", envir = .options)) 3 remov(".revoDoParCluster", where = .options) 4 }
私がregisterDoParallel関数とstopImplicitCluster関数をウロウロしている間に
@hoxo_m先生が原因を突き止めていらっしゃっていて、
@teramonagi @sfchaos @kos59125 これで止まるとこまではつきとめました。
stopCluster(get(".revoDoParCluster", pos = doParallel:::.options))
— hoxo_m (@hoxo_m) 2014, 9月 23さすがだな〜。
なお，以下の本の中でregisterDoParallel関数を使用している箇所は、すべてmakeCluster関数と合わせて使用しているため、上記の問題が顕在化しない．
Juliaで並列計算
Juliaでの並列計算に関する調査メモ。
Juliaでの並列計算の概要
Juliaでのマルチプロセッシング環境は、メッセージパッシングに基づいている。MPIなどの通常のメッセージパッシングは、プロセス間でデータや命令などを相互にやりとりする。しかし、Juliaのメッセージパッシングの実装はあるプロセスから他のプロセスへの一方通行であることが特色となっている。そのため、ユーザは片方のプロセスの管理だけを行えば良い。
あるプロセスから他のプロセスへのメッセージとして重要なのが、"remote call"(remotecall関数)と"remote reference"(fetch関数)である。これらはそれぞれ、「あるプロセスから他のプロセスへの処理の指示」、「指示された処理を行ったプロセスでの処理結果を指示を出したプロセスが参照」に対応する。
- remote call(remotecall関数)
あるプロセスから他のプロセスに引数を渡して関数を実行するための呼び出し。返り値は下記のremote referenceとなる。 - remote reference(fetch関数)
処理を行ったプロセスに保持されたオブジェクトを他のオブジェクトから参照すること。
基本的な実行方法
以下、実行環境はUbuntu14.04(VM Ware上で実行)、Juliaのバージョンは0.3.0。実行結果と次のプロンプトの間の空行は省略。
公式ドキュメントの"Parallel Map and Loops"までを参考にして、Juliaでの並列計算の基本的な実行方法を確認する。
次の例では、remotecall関数を用いて2つのワーカープロセスに一様乱数を生成させ、その値をマスタープロセスがfeatch関数を用いて参照している。
$ ./julia -p 2 julia> # 乱数種の設定 julia> srand(123) julia> # 2つのワーカープロセスで区間[0, 1]の一様乱数を2個ずつ生成 julia> r = remotecall(2, rand, 2, 2) RemoteRef(2,1,4) julia> # 生成した一様乱数の値を取得 julia> fetch(r) 2x2 Array{Float64,2}: 0.341326 0.0673966 0.766994 0.0680754
この「remotecall関数によるワーカープロセスへの処理の指示」+「fetch関数による処理結果の参照」という手間を軽減するために、remotecall_fetch関数が用意されている。
julia> remotecall_fetch(2, getindex, r, 1, 1) 0.3413258405054578
@spawnatマクロを使用すると、第2引数で指定した式を第1引数で指定したワーカープロセスで実行する。次の例では、先に生成した一様乱数に1を加算する処理を行って、その結果をマスタープロセスがfetch関数により参照している。
julia> s = @spawnat 2 1 + fetch(r) RemoteRef(2,1,6) julia> fetch(s) 2x2 Array{Float64,2}: 1.34133 1.0674 1.76699 1.06808
以上の例では、ワーカープロセスの個数を指定してJuliaを起動していた。一方で、Juliaを起動した後に、ワーカープロセスを起動することもできる。
julia> # 起動しているプロセス数(マスタープロセス+ワーカープロセス)の確認 julia> nprocs() 3 julia> nworkers() 2 julia> addprocs(2) 2-element Array{Any,1}: 4 5 julia> # ワーカープロセスが4個に、マスタープロセスとあわせてプロセス数が5個になったことの確認 julia> nworkers() 4 julia> nprocs() 5 julia> # 4番目と5番目のワーカープロセスを停止 julia> rmprocs(4, 5) :ok julia> nworkers() 2
remotecall関数を用いた表記は特に便利でもない。代替として、@spawnマクロを使用することもできる。
julia> r = @spawn rand(2,2) RemoteRef(3,1,10) julia> s = @spawn 1 .+ fetch(r) RemoteRef(3,1,11) julia> fetch(s) 2x2 Array{Float64,2}: 1.03 1.16674 1.99027 1.89547
2行目の処理で、"1 .+r"ではなく"1 .+fetch(r)"にしているのは、どのワーカープロセスでこの処理が実行されるかが分からないためである。
マスタープロセスで定義した関数はワーカープロセスで認識することができないので、単純に@spawnマクロで関数を呼び出すと以下のエラーが発生する。
julia> function rand2(dims...) return 2*rand(dims...) end rand2 (generic function with 1 method) julia> rand2(2,2) 2x2 Array{Float64,2}: 1.5369 1.34792 1.88103 0.790906 julia> @spawn rand2(2,2) RemoteRef(2,1,13) julia> exception on 2: ERROR: function rand2 not defined on process 2 in error at error.jl:21 in anonymous at serialize.jl:397 in anonymous at multi.jl:1279 in anonymous at multi.jl:848 in run_work_thunk at multi.jl:621 in run_work_thunk at multi.jl:630 in anonymous at task.jl:6 julia> @spawn rand2(2,2) RemoteRef(3,1,14) julia> exception on 3: ERROR: function rand2 not defined on process 3 in error at error.jl:21 in anonymous at serialize.jl:397 in anonymous at multi.jl:1279 in anonymous at multi.jl:848 in run_work_thunk at multi.jl:621 in run_work_thunk at multi.jl:630 in anonymous at task.jl:6
すべてのワーカープロセスで変数や関数を認識させて、上記のような事態を回避するには変数や関数の前に@everywhereマクロを付加する。
julia> @everywhere function rand2(dims...) 2 * rand(dims...) end julia> @spawn rand2(2,2) RemoteRef(2,1,20) julia> fetch(r) 2x2 Array{Float64,2}: 1.92811 0.535094 1.40141 0.0465075
同様に、Juliaのソースコードやパッケージをすべてのワーカープロセスに認識させるためには、include関数やusing関数の前に@everywhereを付加する。
@everywhere include("defs.jl") @everywhere using Clustering
Juliaを起動するときに、すべてのワーカープロセスにJuliaのソースコードの記述内容を認識させるためには、-Lオプションでファイル名を指定すれば良い。
$ julia -p <n> -L file1.jl -L file2.jl driver.jl
データの移動
並列計算においては、プロセス間でのデータや命令のやり取りを極力削減することが処理時間の短縮化につながる。
次の例では、の行列の2乗を計算する処理において、
という2つの方法を示している。
julia> # 方法1: マスタープロセスでデータを生成し、ワーカープロセスに転送 julia> A = rand(1000,1000) julia> Bref = @spawn A^2 RemoteRef(3,1,24) julia> fetch(Bref) julia> julia> # 方法2: ワーカープロセスでデータを生成して、マスタープロセスから参照 julia> Bref = @spawn rand(1000,1000)^2 julia> fetch(Bref)
このケースの場合は、方法2ではプロセス間のデータ転送がなく、方法1に比べて効率的である。しかし、いつでも方法2のようにワーカープロセスでデータを生成する方法が良いとは限らず、マスタープロセスや他のワーカープロセスでどのようなデータを使用するかについて考慮しながら、その場その場で適したデータ生成・転送方法を選択する必要がある。
並列Map処理とループ
各ワーカープロセス間でデータのやり取りが必要ない場合は、「@parallel for」構文と「pmap関数」の2つの方法によって、並列計算が可能になる。
次の例では、ワーカープロセスでの論理値の乱数を合計200000000回発生させ(Map処理)、その乱数の和を計算(Reduce処理)している。
julia> nheads = @parallel (+) for i=1:200000000 int(randbool()) end 99999830
このようにfor文の前に@parallelや集約処理を記述することによって、各ワーカープロセスでデータのやり取りを必要としない並列計算が可能になる。この「@parallel for」構文は、1つ1つの処理が軽いときに用いる。
次の例では、ワーカープロセスで成分が区間上の
の行列を10個生成し、pmap関数を用いて特異値分解を実行している。
julia> M = {rand(1000,1000) for i=1:10}
julia> pmap(svd, M)
このように1つ1つの処理が重いときは、pmap関数を用いる。
k平均法
並列計算の応用の一つとして、k平均法を実行してみよう。
Juliaでは、Clusteringパッケージの関数を用いてk平均法を実行する。ここでは、各成分が区間[0, 1]の一様乱数からなるの行列
に対して、k平均法により50個(
)のクラスターを生成する処理を初期値を変えながら200回繰り返してみる。
まずは、逐次処理から。
$ julia julia> using Clustering julia> srand(123) julia> x = rand(1000, 1000) julia> k = 50 julia> @time [result_seq = kmeans(x, k; maxiter=50) for i=1:200]; elapsed time: 49.761205104 seconds (1121608128 bytes allocated, 1.58% gc time)
約50秒計算に要していることが分かる。
続いて、4つのワーカープロセスによる並列計算。
$ julia -p 4 julia> @everywhere using Clustering julia> srand(123) julia> @everywhere x = rand(1000, 1000) julia> @everywhere k = 50 julia> @time @parallel [result_par = kmeans(x, k; maxiter=50) for i=1:200] elapsed time: 20.371857333 seconds (11684024 bytes allocated)
約20秒で計算が完了しているから、2.5倍近くスピードアップしていることが分かる。
まとめと残課題
以上で、基本的な並列計算は実行できるようになったと思われる。一方で、今後調査しなければならない課題もいくつかある。
- 並列計算での乱数生成
RではL'ecuyerのアルゴリズムなど、baseパッケージのRNGkind関数で並列計算に向いた乱数生成法が提供されている。このような並列計算用の乱数生成法が存在しないと、乱数を大量に発生させる必要があるケースでは実用的ではなくなってしまう。Juliaでそのような機能は提供されているのだろうか。おそらくこの問題意識に基づいたソースコードがGithub上で公開されていると思われるので、要調査。 - ロードバランシング
Juliaではタスクのスケジューリングを行っているようであるが(公式ドキュメント"Synchronization With Remote References Scheduling"参照)、ロードバランシングに相当する負荷調整も行っているのだろうか。 - タイミングの測定
Rのsnowパッケージには、マスタープロセスからワーカープロセスに処理やデータを送信したタイミング、ワーカープロセスが処理を行った時間、マスタープロセスがワーカープロセスから処理の結果を受信したタイミングを測定する機能が提供されている。Juliaではこの3つのタイミングのうち、最後のタイミングは存在しないと思われる。前者2つのタイミングを測定することは可能だろうか。 - メッセージパッシング以外の並列計算方法
Rでは、ソケット、フォーク、MPI、PVMなど様々な通信方式が提供されている。Juliaでメッセージパッシング以外の並列計算方法はないだろうか。
参考文献
R2DOCXによるレポート作成
RでMicrosoft Wordのレポートを作成するには，いくつかの方法がある．Markdown+Pandocを使った方法，R2wdパッケージを用いた方法などである．ここでは，David Gohel氏によるR2DOCXパッケージを用いたレポート作成について取り上げる．
インストール
R2DOCXパッケージは，github上で公開されている．R2DOCXとともに，R2DOCパッケージもインストールする．Rのバージョンは3.0.0以降でなければならない模様．
> library(devtools) > install_github('R2DOC', 'davidgohel') > install_github('R2DOCX', 'davidgohel')
R2DOCXパッケージを使用するためには，JavaおよびrJavaパッケージがインストールされている環境が必要．
サンプルプログラム
まずは，チュートリアルのサンプルプログラムを実行してみよう．
library(R2DOCX) # 作成するドキュメント名 docx.file = "document.docx" # ドキュメントの生成(Docx形式) doc = new("Docx", title = "My example") # irisデータセットの先頭10行をドキュメントに追加 doc = addTable(doc, iris[1:10,]) # スタイル名のついたテキストをドキュメントに追加 doc = addParagraph(doc, value="Hello World!", stylename="Normal") # 図をドキュメントに追加 doc = addPlot(doc, function() plot(rnorm(10), rnorm(10)), width=10, height=8) # ドキュメントの出力 writeDoc(doc, docx.file)
この結果，以下の内容のファイルが生成される．
提供する機能
R2DOXパッケージが提供する機能は下表のとおりである．
| 機能 | 関数 |
|---|---|
| 目次の作成 | addTOC |
| 図のプロット | addPlot |
| 外部画像データの挿入 | addImages |
| 表の出力 | addTable |
| テキストの出力 | addParagraph |
| 章の見出しの作成 | addHeader |
| 改ページ | pagebreak |
目次の作成
目次を作成するためには，addTOC関数を使用する．
# 作成するドキュメント名 docx.file = "document_toc.docx" # ドキュメントの生成(Docx形式) doc = new("Docx", title="My example") # 目次の作成 doc = addTOC(doc) # 1章の見出しの作成 doc = addHeader(doc, "Title 1", 1) # 1.1章の見出しの作成 doc = addHeader(doc, "Title 1.1", 2) # 改ページ doc = addPageBreak(doc) # 1.2章の見出しの作成 doc = addHeader(doc, "Title 1.2", 2) # 改ページ doc = addPageBreak(doc) # 1.2.1章の見出しの作成 doc = addHeader(doc, "Title 1.2.1", 3) # ドキュメントの出力 writeDoc(doc, docx.file)
以下のドキュメントが作成される．
図の描画
図を描画するには，addPlot関数を使用する．
doc = addPlot(doc, plot, x=rnorm(10), y=rnorm(10), main="main title", legend="graph example", stylename="PlotReference" )
latticeパッケージやggplot2パッケージを用いて描画する場合は，第2引数にprint関数を指定すれば良い．
doc = addPlot(doc, print, x=qplot(Sepal.Length, Petal.Length, data=iris, color=Species, size=Petal.Width), legend="graph example", stylename="PlotReference" )
両方を組み合わせることも可能．
doc = addPlot(doc, function() { print(qplot(Sepal.Length, Petal.Length, data=iris, color=Species, size=Petal.Width)) plot(x=rnorm(10), y=rnorm(10), main="main title") }, legend = "graph example", stylename="PlotReference" )
外部画像データの挿入
外部の画像データを挿入するには，addImage関数を使用する．
doc = addImage(doc, "D:/img/image.png")
表の出力
表を出力するには，addTable関数を使用する．
ヘッダ名は引数header.labelsにlist型で指定する．
library(R2DOCX) # 新規オブジェクトの作成 doc = new("Docx", title = "My example" ) # irisデータセットの1-5行目の出力 doc = addTable(doc , iris[1:5, ], header.labels = list( "Sepal.Length"="Sepal Length", "Sepal.Width"="Sepal Width", "Petal.Length"="Petal Length", "Petal.Width"="Petal Width", "Species" = "Species" ) ) # ファイルの出力 writeDoc(doc, "document.docx")
ヘッダのラベルをマージするには，grouped.cols引数にマージするヘッダを指定する．以下の例では，"Sepal.Length"と"Sepal.Width"をマージして"Sepal"に，"Petal.Length"と"Petal.Width"をマージして"Petal"にしている．
doc = addTable(doc, iris[1:5,], grouped.cols=list("Sepal"=c("Sepal.Length", "Sepal.Width"), "Petal"=c("Petal.Length", "Petal.Width"), "Species" = c("Species") ) )
同一内容のセルを列方向に結合するには，span.columns引数に結合する列名を指定する．
doc = addTable(doc, iris[46:55, c(5,1:4)], span.columns="Species" )
フォントやセルに色を指定するには，col.fontcolors引数やcol.colors引数を指定すれば良い．
doc = addTable(doc, iris[ 1:5, ] col.fontcolors=list(Species=c("red", "blue", "red", "blue", "red")), col.colors=list(Sepal.Length=c("red", "blue", "red", "blue", "red")) )
テキストの出力
単純なテキストの出力には，addParagraph関数を使用する．
doc = addParagraph(doc, "Hello World!", stylename="Normal")
引数stylenameに使用できるオプションについては，チュートリアルの該当ページを参照のこと．
addParagraph関数を用いて，
# テンプレートのテキスト x = "[animal] eats [food]." # フォーマットのプロパティの定義 textProp=textProperties(color="blue") replacement.styles=list(animal=textProp, food=textProp) # 置換するテキストの定義(animalをdonkeyに，foodをgrassに置換) replacements=list(animal="donkey", food="grass") # テキストの追加(出力は"donkey eats grass"，donkeyとgrassが青字) doc = addParagraph(doc, value=x, stylename="Normal", replacements=replacements, replacement.styles=replacement.styles )
ベクトル化も可能．
x = c("text 1", "text 2") doc = addParagraph(doc, value=x, stylename="Normal")
