Rによる繰返し処理の並列化

統計環境Rによる 繰返しの並列処理 －foreach & doParallel－ 主にWindows機で動作確認した、1台のPC上でソケッ トを利用してメモリ共有型でコアを並列化する方法で すが、他のOSでも動くと思います 2014/3/10

foreachパッケージ • 単体で使用しても繰返し処理を高 速化する。 • snowやparallelなどの並列処理 パッケージと組み合わせると繰返 これが問題! し処理を並列化できる。 • forループとの違い 出力は1つだけ 2

foreach: 出力制御のコツ 制約: 出力は1つだけ • アウトプットの構造や書き出される 順番を理解し、必要があれば後から自 分で配列化などを施す • 複数の出力も順番を間違えなければ 扱える 3

例) 回帰パラメータ算出実験 乱数データから回帰係数を算出する モンテカルロ実験を3回だけ繰返す ‐ xは一様乱数で三変量 ‐ yは次の回帰モデルに従い、自由度3 のt分布に従う乱数を誤差項として付与 y=10+5x1+2x2+3x3 自由度の小さいt分布は正規分布よりも裾が長いので、ちょっと荒 れてる感じのデータになります。結果、そのデータから回帰パラ メータを推計したとき、真の値と差異が出やすいです。 4

40 60 y1 80 100 目的変数 y1 データのイメージ 2 4 6 x1[, 1] 8 第一説明変数 x1 5

並列化の準備として まずforeachの結果の 戻し方を確認します 実験: 1A, 1B, 1C 6

例) 回帰パラメータ算出実験(1A)
注) まだパラレル化はしていません

rm(list=ls(all=TRUE))# ワークスペースのクリア
require(foreach); require(MASS)
# for rt
set.seed(15)
theta1 <- c(10, 5, 2, 3) # 回帰パラメータ "y=10+5x1+2x2+3x3"
ot1 <- foreach(i=1:3, .combine=c) %do% {
x1 <- matrix(runif(300, min=1, max=10), nc=3, nr=100) # xの値
x2 <- cbind(rep(1, dim(x1)[1]), x1) # 1列目に要素が全て1の列を付与
y1 <- x2 %*% theta1 + rt(100, df=3) # yは推計値＋残差
lm(y1 ~ x1)$coeff
}
ot1
is.vector(ot1)
7

例) 回帰パラメータ算出実験(1B)
rm(list=ls(all=TRUE))# ワークスペースのクリア
require(foreach); require(MASS)
# for rt
set.seed(15)
theta1 <- c(10, 5, 2, 3) # 回帰パラメータ "y=10+5x1+2x2+3x3"
ot2 <- foreach(i=1:3, .combine="cbind") %do% {
x1 <- matrix(runif(300, min=1, max=10), nc=3, nr=100) # xの値
x2 <- cbind(rep(1, dim(x1)[1]), x1) # 1列目に要素が全て1の列を
付与
y1 <- x2 %*% theta1 + rt(100, df=3) # yは推計値＋残差
lm(y1 ~ x1)$coeff
}
ot2
is.matrix(ot2)
8

例) 回帰パラメータ算出実験(1C)
rm(list=ls(all=TRUE))# ワークスペースのクリア
require(foreach); require(MASS)
# for rt
set.seed(15)
theta1 <- c(10, 5, 2, 3) # 回帰パラメータ "y=10+5x1+2x2+3x3"
ot3 <- foreach(i=1:3, .combine="rbind") %do% {
x1 <- matrix(runif(300, min=1, max=10), nc=3, nr=100) # xの値
x2 <- cbind(rep(1, dim(x1)[1]), x1) # 1列目に要素が全て1の列を
付与
y1 <- x2 %*% theta1 + rt(100, df=3) # yは推計値＋残差
lm(y1 ~ x1)$coeff
}
ot3
is.list(ot3)
9

結果出力の違い > ot1 (Intercept) x11 x12 x13 (Intercept) x11 10.899300 4.914961 1.978661 2.960806 10.562035 5.005228 x12 x13 (Intercept) x11 x12 x13 1.962990 2.939047 9.786081 5.035381 1.947302 3.102226 > is.vector(ot1) [1] TRUE .combine=c > ot2 result.1 result.2 result.3 (Intercept) 10.899300 10.562035 9.786081 x11 4.914961 5.005228 5.035381 x12 1.978661 1.962990 1.947302 x13 2.960806 2.939047 3.102226 > is.matrix(ot2) > ot3 [1] TRUE ベクトルだったり 行列だったりするけど 出力は一つだけ .combine="cbind" .combine="rbind" (Intercept) x11 x12 x13 result.1 10.899300 4.914961 1.978661 2.960806 result.2 10.562035 5.005228 1.962990 2.939047 result.3 9.786081 5.035381 1.947302 3.102226 > is.matrix(ot3) [1] TRUE 10

doParallelパッケージを使って 同じ処理を並列化してみます 実験: 1A → 2 11

例) 回帰パラメータ算出実験(2)
rm(list=ls(all=TRUE))
# ワークスペースのクリア
s_time <- Sys.time()
# 開始時間記録のための呪文
require(foreach); require(MASS)
# for rt
require(doParallel)
type <- if (exists("mcfork", mode="function")) "FORK" else "PSOCK"
cores <- getOption("mc.cores", detectCores())
cl <- makeCluster(cores, type=type) # コア数を自動検出
registerDoParallel(cl)
# コア数で並列化
set.seed(15)
theta1 <- c(10, 5, 2, 3) # 回帰パラメータ "y=10+5x1+2x2+3x3"
ot1 <- foreach(i=1:3, .combine=c) %dopar% {
x1 <- matrix(runif(300, min=1, max=10), nc=3, nr=100) # xの値
x2 <- cbind(rep(1, dim(x1)[1]), x1) # 1列目に要素が全て1の列を付与
y1 <- x2 %*% theta1 + rt(100, df=3) # yは推計値＋残差
lm(y1 ~ x1)$coeff
}
stopCluster(cl)
# 並列化終了のための呪文
e_time <- Sys.time()
# 終了時間記録のための呪文
ot1
12

コア毎に違う乱数が割り当てられる ために結果は変わる 実験1Aの結果(並列化前) > ot1 (Intercept) 10.899300 x12 1.962990 x11 x12 x13 (Intercept) x11 4.914961 1.978661 2.960806 10.562035 5.005228 x13 (Intercept) x11 x12 x13 2.939047 9.786081 5.035381 1.947302 3.102226 実験2の結果(並列化後) > ot1 (Intercept) 9.495911 x12 1.905662 x11 x12 x13 (Intercept) x11 4.987828 2.035046 3.063866 10.212789 x13 (Intercept) x11 x12 x13 3.046415 8.409125 5.133967 2.103635 5.040721 3.064619 13

で、もちろん速くなるよね？ 実験1Aの結果(並列化前) > e_time - s_time Time difference of 0.5630572 secs 実験2の結果(並列化後) > e_time - s_time Time difference of 10.46505 secs 並列化前が 圧倒的に速い! 実験2のコードの上 から二行目と下から 二行目の時間測定の 呪文を実験1Aの コードにも書き加え て測定しました どんな処理をどのくらいの量だけ並列化するか によります。この例もそうですが、基本的に遅 くなると考えた方が良いです。 ご利用は慎重に。 14

どうして並列化で遅くなるの? 並列化前 I/O 並列化不能部分 並列化可能な部分 ¼に 並列化後 並列化のための処理 リモワ ルフトハンザ 遅くなる犯人! 前の処理の結果に依存する処理は並列化できない 一方、並列化のためには新たな処理(パッケージのロードやコア間通信 など)が必要となる 処理の依存関係のない並列化可能な部分の処理が、例えば4コアで75% 時間を短縮できたとしても、全体のなかでその効率化が新たに必要になる 処理時間を上回る分量でなければトータルで速くはならない 15

アムダールの法則 1 並列化の効果: 1 − 𝑃 + 𝑃/𝑆 P: 並列化可能部分の割合 S: コアの数 並列化可能な部分が全体の50%ならば、どんなにコア数の多い 高価なPCで並列化しても、理論値でさえ2倍の速度が限界です。 実 測 は も っ と ず っ と 効 率 悪 い で す 。 16

TO DO or NOT TO DO 自分一人で作業した方が早いか、そのへ んの人を何人か連れて来て仕事を説明して 作業を分けてでも手伝ってもらった方が早 く終わるような作業内容なのかを見極める。 • 一つずつ順を追って処理しなければいけない仕 事ならば、淡々と一人で片づけるしかない。 • 分業可能な作業であれば、作業の説明や仕分け やまとめなどの新たな手間が増えるのと、分業 の効率を天秤にかけて判断する。 並列処理は、分業させるためのコード書 きも大変なので、それに見合う効率化が 見込めるときにだけ役立つ技なのです。 17

具体的な判断手順 1. コードを精査し、並列処理できる部分とそうでない部分を 切り分けて、並列処理部分は極力少ない塊にまとめる 2. 並列処理できる部分とそうでない部分を別々に時間計測す る 3. 並列処理できる部分が、全体の何パーセントになるのかを 計算し、アムダールの法則のグラフで理論効率を確認 4. 実測効率は理論効率よりもはるかに悪いことを覚悟する それでもやりたければ並列化コードを書く 5. 大量データを取扱う場合、この並列化はメモリ共有型なの で、作業を切り分ける大きさに留意する 大きく切れば早いがメモリ不足の危険があり、小さく切り すぎれば遅くなる まずシングルコア版を作成し、そこから並列化版を作ると良い です。両者の結果が変わらないことが確認できれば安心です。 18

行列要素の並び順: Rとfortranは縦並び > x0 <- 1:100 > (x1 <- matrix(x0, nc=10)) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 11 21 31 41 51 61 71 81 91 [2,] 2 12 22 32 42 52 62 72 82 92 [3,] 3 13 23 33 43 53 63 73 83 93 [4,] 4 14 24 34 44 54 64 74 84 94 [5,] 5 15 25 35 45 55 65 75 85 95 [6,] 6 16 26 36 46 56 66 76 86 96 [7,] 7 17 27 37 47 57 67 77 87 97 [8,] 8 18 28 38 48 58 68 78 88 98 [9,] 9 19 29 39 49 59 69 79 89 99 [10,] 10 20 30 40 50 60 70 80 90 100 > x1[21] [1] 21 ちなみにCは横並びだそうです。 19

配列要素の場合は? > x0 <- 1:1000 > (x1 <- array(x0, c(10,10,10))) ,,1 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 1 11 21 31 41 51 61 71 81 91 [2,] 2 12 22 32 42 52 62 72 82 92 [3,] 3 13 23 33 43 53 63 73 83 93 [4,] 4 14 24 34 44 54 64 74 84 94 [5,] 5 15 25 35 45 55 65 75 85 95 [6,] 6 16 26 36 46 56 66 76 86 96 [7,] 7 17 27 37 47 57 67 77 87 97 [8,] 8 18 28 38 48 58 68 78 88 98 [9,] 9 19 29 39 49 59 69 79 89 99 [10,] 10 20 30 40 50 60 70 80 90 100 ,,2 20

なぜ並び順が重要? • 変数一つの出力の中に 複数の情報を詰め込む 場合、どの情報がどう いう順番に入っている かわかっていないと使 えない • 高速なコードのために、 極力キャッシュミスが 起こりにくいような処 理を心がける CPU レジスタ 高速 高価 小容量 キャッシュ メモリ ディスク 低速 安価 大容量 参考: 「チューニング技法虎の巻」 https://www.hpci-office.jp/pages /407 21

• https://www.hpci-office.jp/pages/407

並列化によるシミュレーション 実験: 2 → 3 22

例) 回帰パラメータ算出実験(3)

rm(list=ls(all=TRUE))
# ワークスペースのクリア
s_time <- Sys.time()
# 開始時間記録のための呪文
require(foreach); require(MASS)
# for rt
require(doParallel)
type <- if (exists("mcfork", mode="function")) "FORK" else "PSOCK"
cores <- getOption("mc.cores", detectCores())
cl <- makeCluster(cores, type=type)
# コア数を自動検出
registerDoParallel(cl)
# コア数で並列化
set.seed(15)
theta1 <- c(10, 5, 2, 3) # 回帰パラメータ "y=10+5x1+2x2+3x3"
ot1 <- foreach(i=1:4, .combine=c) %dopar% {
# コア数の倍数が効率良いです
cf2 <- matrix(NA, nc=2500, nr=4) # 縦に並べないとわけがわからなくなる!!!
for (j in 1:2500) {
x1 <- matrix(runif(300, min=1, max=10), nc=3, nr=100) # xの値
x2 <- cbind(rep(1, dim(x1)[1]), x1) # 1列目に要素が全て1の列を付与
y1 <- x2 %*% theta1 + rt(100, df=3) # yは推計値＋残差
cf2[,j] <- lm(y1 ~ x1)$coeff
}
cf2
}
stopCluster(cl)
# 並列化終了のための呪文
e_time <- Sys.time()
# 終了時間記録のための呪文
ot1
23

比較用のシングルコア版コード(1A’)
rm(list=ls(all=TRUE))
# ワークスペースのクリア
s_time <- Sys.time()
# 開始時間記録のための呪文
require(foreach); require(MASS)
# for rt
set.seed(15)
theta1 <- c(10, 5, 2, 3) # 回帰パラメータ "y=10+5x1+2x2+3x3"
ot1 <- foreach(i=1:4, .combine=c) %do% {
cf2 <- matrix(NA, nc=2500, nr=4) # 縦に並べないとわけがわからなくなる!!!
for (j in 1:2500) {
x1 <- matrix(runif(300, min=1, max=10), nc=3, nr=100) # xの値
x2 <- cbind(rep(1, dim(x1)[1]), x1) # 1列目に要素が全て1の列を付与
y1 <- x2 %*% theta1 + rt(100, df=3) # yは推計値＋残差
cf2[,j] <- lm(y1 ~ x1)$coeff
}
cf2
}
e_time <- Sys.time()
# 終了時間記録のための呪文
e_time - s_time
24

並列化で嬉しくなるのは? • 400回 シングル Time difference of 2.374 secs 並列化 Time difference of 8.911 secs • 4000回 シングル Time difference of 7.158 secs 並列化 Time difference of 11.222 secs • 10000回 シングル Time difference of 15.31853 secs 並列化 Time difference of 15.63356 secs • 40000回 シングル Time difference of 55.244 secs 並列化 Time difference of 36.901 secs • 100000回 シングル Time difference of 2.264043 mins 並列化 Time difference of 1.329983 mins Windows VISTA, 32bit, R3.0.0, コア数2の条件でテストした結果 25

結論: シミュレーションに使える!
データの取り出し方
length(ot1)
[1] 40000
> ot1x <- matrix(ot1, nr=4) # 整形
> dim(ot1x)
[1]
4 10000
> apply(ot1x, 1, mean)
[1] 10.008551 4.999126 1.999420 3.000582
> apply(ot1x, 1, median)
[1] 10.010047 4.998599 1.999008 3.000921
# モデルの理論値は10, 5, 2, 3
26

参考資料 Revolution Analytics. “doParallel: Foreach parallel adaptor for the parallel package,” R package version 1.0.1. 2012, available at http://CRAN.R-project.org/package=doParallel. Revolution Analytics, “foreach: Foreach looping construct for R,” R package version 1.4.0. 2012, available at http://CRAN.Rproject.org/package=foreach. S. Weston and R. Calaway, “Getting started with doParallel and foreach,” 2010, available at http://cran.rproject.org/web/packages/doParallel/vignettes/gettingstartedParall el.pdf. S. Weston, “Using the foreach packages,” 2009, available at http://www.rochester.edu/college/psc/thestarlab /help/UsingTheFo reachPackage.pdf. おわり 27

• http://cran.r-project.org/package=doParallel
• http://cran.r-project.org/package=foreach
• http://cran.r-project.org/web/packages/doParallel/vignettes/gettingstartedParallel.pdf
• http://www.rochester.edu/college/psc/thestarlab/help/UsingTheForeachPackage.pdf