---
title: Rによる繰返し処理の並列化
tags:  #r #並列化  
author: [wada, kazumi](https://convert.docswell.com/user/kazwd2008)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/DEYVWLVYEM.jpg?width=480
description: forやapplyによる繰り返し処理を、パッケージforeachとdoParallelにより並列化する方法の解説
published: December 19, 23
canonical: https://convert.docswell.com/s/kazwd2008/Z4QWP1-2023-12-19-084649
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/DEYVWLVYEM.jpg)

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


# Page. 2

![Page Image](https://bcdn.docswell.com/page/VJNM94MR78.jpg)

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


# Page. 3

![Page Image](https://bcdn.docswell.com/page/YE9N2QNZE3.jpg)

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


# Page. 4

![Page Image](https://bcdn.docswell.com/page/GE8K5GNY7D.jpg)

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


# Page. 5

![Page Image](https://bcdn.docswell.com/page/LELZYGP9JR.jpg)

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


# Page. 6

![Page Image](https://bcdn.docswell.com/page/4JMVNQ4VJW.jpg)

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


# Page. 7

![Page Image](https://bcdn.docswell.com/page/PJR2D8WWE9.jpg)

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


# Page. 8

![Page Image](https://bcdn.docswell.com/page/PEXD18PV7X.jpg)

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


# Page. 9

![Page Image](https://bcdn.docswell.com/page/3EKX2KD57D.jpg)

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


# Page. 10

![Page Image](https://bcdn.docswell.com/page/L7354ZR1E5.jpg)

結果出力の違い
&gt; 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
&gt; is.vector(ot1)
[1] TRUE
.combine=c
&gt; 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
&gt; is.matrix(ot2)
&gt; ot3
[1] TRUE
ベクトルだったり
行列だったりするけど
出力は一つだけ
.combine=&quot;cbind&quot;
.combine=&quot;rbind&quot;
(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
&gt; is.matrix(ot3)
[1] TRUE
10


# Page. 11

![Page Image](https://bcdn.docswell.com/page/87DDQRYK7G.jpg)

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


# Page. 12

![Page Image](https://bcdn.docswell.com/page/VJPXLW6378.jpg)

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


# Page. 13

![Page Image](https://bcdn.docswell.com/page/2EVLQ83NJQ.jpg)

コア毎に違う乱数が割り当てられる
ために結果は変わる
実験1Aの結果(並列化前)
&gt; 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の結果(並列化後)
&gt; 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


# Page. 14

![Page Image](https://bcdn.docswell.com/page/57GXW535EL.jpg)

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


# Page. 15

![Page Image](https://bcdn.docswell.com/page/4EQ83Z5LJP.jpg)

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


# Page. 16

![Page Image](https://bcdn.docswell.com/page/KJ421395E1.jpg)

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


# Page. 17

![Page Image](https://bcdn.docswell.com/page/LE1MG1K27G.jpg)

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


# Page. 18

![Page Image](https://bcdn.docswell.com/page/GEW3K8D2E2.jpg)

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


# Page. 19

![Page Image](https://bcdn.docswell.com/page/47ZVZ8G4E3.jpg)

行列要素の並び順:
Rとfortranは縦並び
&gt; x0 &lt;- 1:100
&gt; (x1 &lt;- 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
&gt; x1[21]
[1] 21
ちなみにCは横並びだそうです。
19


# Page. 20

![Page Image](https://bcdn.docswell.com/page/YJ65ZPYG7V.jpg)

配列要素の場合は?
&gt; x0 &lt;- 1:1000
&gt; (x1 &lt;- 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


# Page. 21

![Page Image](https://bcdn.docswell.com/page/GJ5YWKGXJ4.jpg)

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


# Page. 22

![Page Image](https://bcdn.docswell.com/page/LE354ZR4J5.jpg)

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


# Page. 23

![Page Image](https://bcdn.docswell.com/page/8EDDQRY5EG.jpg)

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


# Page. 24

![Page Image](https://bcdn.docswell.com/page/V7PXLW6DE8.jpg)

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


# Page. 25

![Page Image](https://bcdn.docswell.com/page/2JVLQ83G7Q.jpg)

並列化で嬉しくなるのは?
• 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


# Page. 26

![Page Image](https://bcdn.docswell.com/page/5EGXW5YDJL.jpg)

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


# Page. 27

![Page Image](https://bcdn.docswell.com/page/4JQ83ZMX7P.jpg)

参考資料
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


