狠狠撸

狠狠撸Share a Scribd company logo
繰返しの並列処理
-foreach & doParallel-
Windowsで開発しましたが, Mac
OSでもLinuxでも動くと思います
統計環境Rによる
2014/3/10
foreachパッケージ
? 単体で使用しても繰返し処理を高
速化する。
? snowやparallelなどの並列処理
パッケージと組み合わせると繰返
し処理を並列化できる。
? forループとの違い
出力は1つだけ
これが問題!
2
foreach: 出力制御のコツ
制約: 出力は1つだけ
? アウトプットの構造や書き出される
順番を理解し、必要があれば後から自
分で配列化などを施す
? 複数の出力も順番を間違えなければ
扱える
3
例) 回帰パラメータ算出実験
乱数データから回帰係数を算出する
モンテカルロ実験を3回だけ繰返す
‐ xは一様乱数で三変量
‐ yは次の回帰モデルに従い、自由度3
のt分布に従う乱数を誤差項として付与
y=10+5x1+2x2+3x3
自由度の小さいt分布は正規分布よりも裾が長いので、ちょっと
荒れてる感じのデータになります。結果、そのデータから回帰パ
ラメータを推計したとき、真の値と差異が出やすいです。 4
2 4 6 8
40
60
80
100
x1[, 1]
y1
データのイメージ
第一説明変数 x1
目的変数 y1
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
> 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)
[1] TRUE
> ot3
(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
.combine=c
.combine="cbind"
.combine="rbind"
ベクトルだったり
行列だったりするけど
出力は一つだけ
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
コア毎に違う乱数が割り当てられる
ために結果は変わる
> ot1
(Intercept) x11 x12 x13 (Intercept) x11
9.495911 4.987828 2.035046 3.063866 10.212789 5.040721
x12 x13 (Intercept) x11 x12 x13
1.905662 3.046415 8.409125 5.133967 2.103635 3.064619
> 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
実験1Aの結果(並列化前)
実験2の結果(並列化後)
13
で、もちろん速くなるよね?
どんな処理をどのくらいの量だけ並列化するか
によります。この例もそうですが、基本的に遅
くなると考えた方が良いです。 ご利用は慎重に。
> e_time - s_time
Time difference of 0.5630572 secs
実験1Aの結果(並列化前)
> e_time - s_time
Time difference of 10.46505 secs
実験2の結果(並列化後)
並列化前が
圧倒的に速い!
実験2のコードの上
から二行目と下から
二行目の時間測定の
呪文を実験1Aの
コードにも書き加え
て測定しました
14
どうして並列化で遅くなるの?
前の処理の結果に依存する処理は並列化できない
一方、並列化のためには新たな処理(パッケージのロードやコア間通信な
ど)が必要となる
処理の依存関係のない並列化可能な部分の処理が、例えば4コアで75%
時間を短縮できたとしても、全体のなかでその効率化が新たに必要になる
処理時間を上回る分量でなければトータルで速くはならない
リモワ ルフトハンザ
並列化前
I/O 並列化不能部分 並列化可能な部分
?に 並列化のための処理
遅くなる犯人!
並列化後
15
アムダールの法則
並列化可能な部分が全体の50%ならば、どんなにコア数の多い
高価なPCで並列化しても、理論値でさえ2倍の速度が限界です。
P: 並列化可能部分の割合
S: コアの数
1
1 ? ? + ?/?
並列化の効果:
実
測
は
も
っ
と
ず
っ
と
効
率
悪
い
で
す
。
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
なぜ並び順が重要?
? 変数一つの出力の中に
複数の情報を詰め込む
場合、どの情報がどう
いう順番に入っている
かわかっていないと使
えない
? 高速なコードのために、
極力キャッシュミスが
起こりにくいような処
理を心がける
21
キャッシュ
ディスク
メモリ
レジスタ
CPU
低速
安価
大容量
高速
高価
小容量
参考: 「チューニング技法虎の巻」 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
25
Windows VISTA, 32bit, R3.0.0, コア数2の条件でテストした結果
結論: シミュレーションに使える!
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
データの取り出し方
おわり
27
参考資料
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.R-
project.org/package=foreach.
S. Weston and R. Calaway, “Getting started with doParallel and foreach,”
2010, available at http://cran.r-
project.org/web/packages/doParallel/vignettes/gettingstartedParallel.pdf.
S. Weston, “Using the foreach packages,” 2009, available at
http://www.rochester.edu/college/psc/thestarlab/help/UsingTheForeachPac
kage.pdf.

More Related Content

搁による繰り返しの并列処理