RのdoMCによる並列計算（マルチコア化）のメモ

R言語で計算高速化のためにマルチコアCPUで並列計算したい場合がある。並列計算と聞いて抵抗を感じる人も多いが、Rに関しての並列計算プログラムは容易に組むことが可能だ。
ではどうやって並列計算のコードを組むのか、モンティ・ホール問題の試行を計算するプログラムを元にやってみる。

モンティ・ホール問題（Monty Hall problem）とはアメリカのテレビ番組「Let's make a deal」(1984-1986)で行われたゲームに関する確率論の問題である。ゲームの内容は次のとおりである。

  
1. 3つのドア (A, B, C) に（景品、ヤギ、ヤギ）がランダムに入っている。
  
2. プレイヤーはドアを1つ選ぶ。
  
3. モンティは残りのドアのうち1つを必ず開ける。
  
4. モンティの開けるドアは、必ずヤギの入っているドアである。
  
5. モンティはプレーヤーにドアを選びなおしてよいと必ず言う。

(モンティ・ホール問題及びそのゲーム内容に関する引用:モンティ・ホール問題 - Wikipedia)

このルールに則り、ゲームを再現するプログラムを書いた。

montyHall <- function(changeAns = FALSE){
  doorNum <- 1:3
  rightAns <- sample(doorNum, size = 1)
  playerChoice <- sample(doorNum, size = 1)
  montyChoice <- doorNum[doorNum != rightAns]
  montyChoice <- montyChoice[montyChoice != playerChoice]
  montyChoice <- sample(montyChoice, size = 1)
  if(changeAns){
    playerChoice <- doorNum[doorNum != playerChoice]
    playerChoice <- playerChoice[playerChoice != montyChoice]
    playerChoice <- sample(playerChoice, size = 1)
  }
  if(playerChoice == rightAns){result <- 1}else{result <- 0}
  result_data <- c(changeAns, result)
  return(result_data)
}

プログラムの便宜上、ドアの名前を(A, B, C)から(1,2,3)にしている。
また、返り値として、(選びなおす(1)か否か(0),景品か(1)否か(0))としている。
さて、これを通常のforを使った書き方でドアを選びなおす場合と選びなおさない場合をそれぞれ10万回行い、その時間を測る。コードは以下の通り

####### SC code #######
montyResult_Change <- NULL
SC_time_Change <- system.time(
  for(i in 1:repeatNum){
    montyResult_Change <- rbind(montyResult_Change, montyHall(TRUE))
  }
)

montyResult_NoChange <- NULL
SC_time_NoChange <- system.time(
  for(i in 1:repeatNum){
    montyResult_NoChange <- rbind(montyResult_NoChange, montyHall())
  }
)
cat("-----Result-----\nrepeat:",repeatNum,"\nChange    :",sum(montyResult_Change[,2]),"\nNo Change :",sum(montyResult_NoChange[,2]),"\n")

それぞれにかかった時間は、

                 user.self sys.self elapsed user.child sys.child
SC_time_Change      54.284    1.885  56.233      0.000     0.000
SC_time_NoChange    56.459    0.000  56.529      0.000     0.000

のようになった。

次に、マルチコア化したプログラムで同じように10万回試行する。くわしい書き方は後述する。

####### MC code #######
#install.packages("doMC")
library(doMC)
registerDoMC(cores=4)

MC_time_Change <- system.time(
  montyResult_Change <- foreach(i=1:repeatNum,.combine=rbind)%dopar%{
    montyHall(TRUE)
  }
)
MC_time_NoChange <- system.time(
  montyResult_NoChange <- foreach(i=1:repeatNum,.combine=rbind)%dopar%{
    montyHall()
  }
)
cat("-----Result-----\nrepeat:",repeatNum,"\nChange    :",sum(montyResult_Change[,2]),"\nNo Change :",sum(montyResult_NoChange[,2]),"\n")

それぞれにかかった時間は、

                 user.self sys.self elapsed user.child sys.child
MC_time_Change      35.387    0.028  36.517      3.148     0.144
MC_time_NoChange    46.255    0.040  47.395      3.992     0.252

のように、計算速度は約1.5倍ほど速くなった。

では，マルチコア化の手順を書いていく，
まず，packageを読み込む．私の環境はUbuntuまたはMacOS Xなので，packageはdoMCを使用する．Rに読み込ませる場合は以下のコードを実行する，

install.packages("doMC")
library(doMC)

このコードで，依存する（？）関連package（i.e. foreach etc...）がなければ自動的にインストール・読込される．次に，並列計算で使用するcore数を指定する．今回は4core使って計算するので次のように指定する．

registerDoMC(cores=4)

この３つのコードを書いたあとに，foreach文を使って並列計算の内容を書いていく．
並列計算の使用シーンはfor構文で書かれた繰り返し計算を並列化しjobが無いcoreにどんどんjobを割り振りたい時である．
例えば上記の例で言えば，

for(i in 1:repeatNum){
  montyResult_Change <- rbind(montyResult_Change, montyHall(TRUE))
}

という部分を並列にどんどん計算して行きたいわけである．ここで使う関数（構文？）がforeachである．
foreach文はfor構文に少しにているがfor構文と挙動が異なるので注意が必要である．例えば，for文で次のようなコードを実行しようとしよう．

a <- NULL
for(i in 1:10){
  a <- c(a,i)
}

これはベクトルaに1から10までの数字をaの末尾にどんどん追加していくコードである．for文が終わればこのaは

> a
[1] 1 2 3 4 5 6 7 8 9 10

という結果が帰ってくる．これをforeach文（シングルコアver）であらわそう．

a <- foreach(i = 1:10, .combine=c)%do%{
  i
}

シンプルな文ではあるが，色々とforとは違う部分があってわかりにくい部分もあると思う．私がforeachを最初に使った時に疑問に思った部分を中心に解説をしていく．
まず，foreach文は文の{}内で何かの変数に代入してもその変数の情報は{}をその実行している時だけ保持するように設計されている．だから，

a <- NULL
foreach(i = 1:10, .combine=c)%do%{
  a <- i
}

のように書いてもforeach文が終わった後のaにはNULLしか入っていないことに注意が必要である．しかし，foreach文で計算した結果はforeach文の中だけで外部に出力（何かの変数に代入など）することが出来ないかと言うわけではない．外部に出力をしたい場合にはforeachの前にその結果を順次追加していく変数をけば良い．

a <- foreach(i = 1:10, .combine=c)%do%{
  i
}

順次追加していく値は，foreachの{}内で追加したい数をterminalで確認するように打てば良い．例えば，

a <- foreach(i = 1:10, .combine=c)%do%{
  result <- i*10
}

という文をかいたとしても，resultの値はaに追加されない．変数resultの値をaに使いしたい場合は，

a <- foreach(i = 1:10, .combine=c)%do%{
  result <- i*10
  result
}

と組めば良い．ところで，以下のコードではどの変数がaに代入されるだろうか？

a <- foreach(i = 1:10, .combine=c)%do%{
  result <- i*10
  result2 <- i^10
  result
  result2
}

正解はresult2がaに追加される．つまり，一番最後に書いた（確認した？）変数がaに代入される．
次に，".combine" optionに関して書く．
.combineはaに随時代入される形式を指定する．各形式についてRのコードで例えてみよう，

===================================

a <- foreach(i = 1:10,.combine=c)%do%{
  i
}
ならば，
a <- NULL
for(i in 1:10{
  a <- c(a,i)
}

===================================

a <- foreach(i = 1:10,.combine=rbind)%do%{
  i
}
ならば，
a <- NULL
for(i in 1:10{
  a <- rbind(a,i)
}

===================================

a <- foreach(i = 1:10,.combine=cbind)%do%{
  i
}
ならば，
a <- NULL
for(i in 1:10{
  a <- cbind(a,i)
}

===================================

ちなみに，.combineにとくに指定がない場合はlist形式で随時追加される．
あとは次の2点がforとforeachと違う点である．

• for文の"in"はforeachの"="である．
• for(){}はforeach()%do%{}

さて，シングルコア版のforeach文をもう一度だそう．

a <- foreach(i = 1:10, .combine=c)%do%{
  i
}

これのマルチコア化したコードは次のとおりである．

a <- foreach(i = 1:10, .combine=c)%dopar%{
  i
}

つまり，"%do%"を"%dopar%"に変更するだけである．
ちなみに，あまりにも繰り返し回数が少ない，繰り返し構文中の計算量が少ない場合はマルチコア化した構文を書いても遅くなったり，無駄な計算が増えるので注意が必要である．また，例えば次のようなコードを書いたとする，

a <- foreach(i=1:10, .combine=c)%dopar%{
  cat(i,"")
  cat(i,"")
  i
}

これはiの内容を実行中に逐次terminalに表示させる構文であるが，普通のfor構文のような考えなら，

1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10

というように表示されるが，実際には

1 1 5 5 9 9 2 2 6 6 10 10 3 3 7 7 4 4 8 8

のように表示されてしまう（一例）．ちなみにaには

1 2 3 4 5 6 7 8 9 10

というデータが代入される．

なぜこのようにバラバラに表示されてしまうかというと，foreachの{}内のjobは開いたコアから順次割り振られ，コアによってjobの実行速度は変わってくるためである．forで作った時にcat()やprint()を入れているコードはマルチコア化するときにこの点が要注意である．

ちなみに，montyHall()の実行結果だが，ドアを選び直したほうが景品が当たりやすいようである．

詳しい参考：
Package ‘doMC’
Getting Started with doMC and foreach