この文章は、foreach パッケージの Vignette
Using The foreach Package
http://cran.r-project.org/web/packages/foreach/vignettes/foreach.pdf
を適当に翻訳したものです。
1 イントロダクション
R の最も便利な機能の一つとして、インタラクティブなインタプリタが挙げられます。これは、R を学んだり試しに使ってみたりするのを非常に簡単にします。これがあることで、R は電卓のように気軽に使うことができ、算術演算、データセットの表示、プロットの生成、モデルの作成を行うことができます。
新しく R を始めたユーザは、いずれそのうち、いくつかの操作を繰り返し実行する必要が出てきます。ひょっとすると、結果の分布を知るために繰り返しシミュレーションを走らせたいのかもしれません。ひょっとすると、ある関数に様々なバリエーションの引数を渡して実行したいのかもしれません。ひょっとすると、異なるデータセットに対してモデルを作成したいのかもしれません。
反復実行は、手動で行うこともできます。しかし、繰り返しの動作を実行するのは、たとえコマンドライン編集を使ったとしても、きわめて退屈な作業となります。幸いにも、R は電卓よりもはるかに高機能です。R は、計算を繰り返し実行するというような退屈な仕事を、自動的に行うための組み込みの構文を持っています。
R には、ループを構築する様々な方法があります。ループを使えばこの問題は解決します。for ループは最もよく使われるループの方法です。また、repeat 文や while 文も非常に便利です。加えて、apply, lapply, sapply, eapply, mapply, rapply のような、"apply" 関数ファミリーがあります。
foreach パッケージは、R コードを繰り返し実行するための新たなループ作成法を提供します。R におけるループ作成法が、すでにどれを使えばいいか困惑してしまうほど存在するのに、さらに別の方法が本当に必要なのかと疑問を持つかもしれません。foreach パッケージを使用する主な理由は、foreach が並列実行をサポートするためです。これにより、繰り返し処理を複数のプロセッサやコア、またはクラスタ上のノードで実行することができます。もし、繰り返しの各処理が一分以上かかり、それを 100 回以上繰り返したい場合、全体の実行時間は数時間かかることになります。しかし、foreach を使えば、各処理は何百ものプロセッサ上で並列に実行できるため、実行時間を数分にまで削減することができます。
しかし、並列実行だけが foreach パッケージを使用する理由ではありません。後で見るように、時間のかからない処理に対しても、foreach を使用する理由はいくつかあります。
2 はじめよう
foreach パッケージの簡単な使用例を見てみましょう。foreach パッケージはすでにインストールされているとします。まず初めにパッケージをロードする必要があります。
library(foreach)
これにより、foreach が依存するパッケージもロードされることに注意して下さい。
今、foreach を、1 から 3 の値に対して sqrt 関数を繰り返し適用するために使います。結果はリストとして返し、変数 x に代入します1：
x <- foreach(i=1:3) %do% sqrt(i)
x
## [[1]]
## [1] 1
##
## [[2]]
## [1] 1.414214
##
## [[3]]
## [1] 1.732051
これは、for ループに似ていますが、%do% という二項演算子によって書かれているため、少し奇妙に見えます。また、for ループとは異なり、値を返します。この点は非常に重要です。このコードの目的は、計算した結果のリストを返すことにあります。一般的に、%do% を使った foreach は、R の表現式を繰り返し実行するために使われます。結果は、なんらかのデータ構造か、オブジェクトとして返されます。デフォルトではリストを返します。
この例では、sqrt の引数として変数 i を使用したことに注意して下さい。変数 i の値は foreach 関数に対して名前付き引数を使って指定しました。この変数には、例えば a や b など、任意の名前をつけることができます。また、次の例のように、複数の変数を指定して、R の表現式内で使うことができます。
x <- foreach(a=1:3, b=rep(10, 3)) %do% (a + b)
x
## [[1]]
## [1] 11
##
## [[2]]
## [1] 12
##
## [[3]]
## [1] 13
この場合、括弧が必要になることに注意して下さい。中括弧でも OK です：
x <- foreach(a=1:3, b=rep(10, 3)) %do% {
a + b
}
x
## [[1]]
## [1] 11
##
## [[2]]
## [1] 12
##
## [[3]]
## [1] 13
ここでは、a と b が反復変数として使われています。これらの変数は、反復処理の間に変化します。これらの変数が同時に変化することに注意して下さい。
今回の場合、2 つの反復変数の長さは同じですが、必ずしもそうである必要はありません。もし、b に 2 つの値しか与えなかった場合、たとえ a に 1000 個の値を与えたとしても、結果のリストは長さ 2 になります：
x <- foreach(a=1:1000, b=rep(10, 2)) %do% {
a + b
}
x
## [[1]]
## [1] 11
##
## [[2]]
## [1] 12
中括弧の中に複数の行を書いてもよいことに注意して下さい。また、計算の途中で生じた結果を変数に代入することも可能です。しかし、もし代入された変数を、ループのイテレーション間で受け渡ししている場合、そのコードは並列実行するとうまく動きません。これについてはあとで言及します。
3 .combine オプション
ここまでの例では、結果はすべてリストで返されました。リストにはすべての R オブジェクトを含めることができるので、デフォルトとして適切です。しかし、例えば、結果を数値ベクトルとして返したい場合もあります。これは、.combine オプションを使うことで可能です。
x <- foreach(i=1:3, .combine='c') %do% exp(i)
x
## [1] 2.718282 7.389056 20.085537
この結果は数値ベクトルとして返されます。なぜなら、すべての結果を結合するために、R 標準の c 関数が使われているからです。exp 関数は数値を返すため、それを c 関数で結合した結果は長さ 3 の数値ベクトルとなります。
R の表現式がベクトルを返す場合、そのベクトルを行列にするにはどうすればいいでしょうか？それには、cbind 関数が使えます：
x <- foreach(i=1:4, .combine='cbind') %do% rnorm(4)
x
## result.1 result.2 result.3 result.4
## [1,] -0.003413306 0.08932265 -1.01913551 -0.2253865
## [2,] 0.049261237 0.09414726 0.03778331 1.2991611
## [3,] -0.059788349 0.60224948 -0.80578526 0.5049468
## [4,] -1.406062000 -0.86282610 -1.17449742 -0.2191800
この表現式は、4 つの乱数からなるベクトルを返します。それを列として結合することで、4×4 の行列になります。
また、結合関数として "+" や "*" も使えます：
x <- foreach(i=1:4, .combine='+') %do% rnorm(4)
x
## [1] -0.3357704 5.2691451 1.9543764 1.3418343
また、結合関数に、ユーザ定義関数を指定することもできます。この例では、計算した結果を捨てています：
cfun <- function(a, b) NULL
x <- foreach(i=1:4, .combine='cfun') %do% rnorm(4)
x
## NULL
この cfun 関数が、2 つの引数を取ることに注目して下さい。foreach 関数は、c, cbind, rbind 関数が引数をいくつでも取れることを知っており、(デフォルトでは)最大 100 の引数を取ることで、パフォーマンスの改善を行っています。しかし、(例えば "+") のように、他の関数については、2 つの引数しか取れないと仮定して動作します。結合関数が任意の数の引数を取ることができる場合、.multicombine 引数によって指定できます：
cfun <- function(...) NULL
x <- foreach(i=1:4, .combine='cfun', .multicombine=TRUE) %do% rnorm(4)
x
## NULL
結合関数を 10 以下の引数で呼び出したい場合、.maxcombine 引数を使って設定できます：
cfun <- function(...) NULL
x <- foreach(i=1:4, .combine='cfun', .multicombine=TRUE, .maxcombine=10) %do% rnorm(4)
x
## NULL
.inorder オプションは、結果が結合される順序が引数の順序と同じでなければならないかどうかを指定するのに使います。デフォルトは TRUE ですが、もし結合関数が "+" の場合は、.inorder を FALSE にしても大丈夫なことがわかるでしょう。実際には、このオプションが重要になるのは、並列実行した場合に限ります。なぜなら、逐次実行した場合は、結果は常に順序通りに計算されるからです。並列実行の場合、結果は必ずしも順序通りに出るとは限りません。実際、各イテレーションは異なる長さの実行時間を要し、計算結果が出る順序はどのようにでもなりえます。ここでは、違いを示すために、並列実行を行う不自然な例をみてみましょう。この例は、Sys.sleep 関数を使って、早いイテレーションほど長い実行時間を持たせるようにしています。
foreach(i=4:1, .combine='c') %dopar% {
Sys.sleep(3 * i)
i
}
## [1] 4 3 2 1
foreach(i=4:1, .combine='c', .inorder=FALSE) %dopar% {
Sys.sleep(3 * i)
i
}
## [1] 4 3 2 1
2 つのうち最初の結果は、ベクトル c(4, 3, 2, 1) になることが保証されます。2 つ目の結果は、同じ値になっていますが、違う順序になるかもしれません。
4 イテレータ
反復変数には、ベクトルかリストしか指定できないということはありません。反復変数には、「イテレータ」であれば何でも指定できます。イテレータは、iterators パッケージにたくさん用意されています。イテレータとは、抽象化されたデータソースです。ベクトルはイテレータではありませんが、foreach 関数は、ベクトルから自動的にイテレータを作成します。
また、foreach 関数は、リスト、行列、データフレームからもイテレータを自動的に作成します。自然なデータソースである、ファイルやデータベースクエリからもイテレータを作成することができます。iterators パッケージは、呼ばれるたびに指定された数の乱数を返す irnorm 関数を提供します。例えば：
library(iterators)
x <- foreach(a=irnorm(4, count=4), .combine='cbind') %do% a
x
## result.1 result.2 result.3 result.4
## [1,] -1.1290423 -0.09302933 0.8165786 0.8636597
## [2,] 0.3269669 1.61575052 0.6425222 -0.2513625
## [3,] 0.1863954 0.42819766 -1.8843529 0.3478381
## [4,] -0.7434539 1.27246768 -0.6020819 0.6718292
これは、大量のデータを扱う際に便利です。イテレータは、最初に全てのデータを生成するのではなく、計算で必要になったときにオンザフライでデータを生成することを可能にします。
例えば、1000 個のランダムベクトルの合計を求めたいとします：
set.seed(123)
x <- foreach(a=irnorm(4, count=1000), .combine='+') %do% a
x
## [1] 9.097676 -13.106472 14.076261 19.252750
これは非常に少ないメモリ量で済みます。なぜなら、このコードは次の while ループと等価だからです：
set.seed(123)
x <- numeric(4)
i <- 0
while (i < 1000) {
x <- x + rnorm(4)
i <- i + 1
}
x
## [1] 9.097676 -13.106472 14.076261 19.252750
これは、icount 関数を使って、1 から 1000 までの値を生成することによっても可能です。
set.seed(123)
x <- foreach(icount(1000), .combine='+') %do% rnorm(4)
x
## [1] 9.097676 -13.106472 14.076261 19.252750
ただし、このような方法よりも、実際に使用するデータをイテレータで生成する方が望ましい方法です(あとで並列実行のときに見ます)。
最後の例では、iterators パッケージの icount 関数を導入したのに加えて、foreach に引数を名前なしで渡す例にもなっています。これは、値を生成せずに、R の表現式を実行する回数を制御したい場合に便利です。
iterators についてもっと多くのことが言えますが、今はとりあえず並列実行の説明に移りましょう。
5 並列実行
foreach はそれ自体便利なパッケージではありますが、このパッケージの真の利点は並列計算を行うことにあります。これまでに見てきた例を並列実行に変更するには、%do% を %dopar% に置き換えるだけです。ただし、上で行ってきたような、すぐに終わってしまう計算に関しては、並列実行する利点がないでしょう。大量の小さな処理の実行を並列化したとしても、逐次実行した場合に比べてもっと時間がかかるようになるだけです。さらに言えば、それらはすでに十分速いので、速くするモチベーションがありません。しかし、並列に実行しても 1 分以上かかる処理については、モチベーションを持ち始めるでしょう。
5.1 並列ランダムフォレスト
時間のかかる処理の例として、ランダムフォレストを見てみましょう。入力は、行列 x とファクター y とします：
x <- matrix(runif(500), 100)
y <- gl(2, 50)
foreach パッケージのロードはすでに行っているとして、randomForest パッケージをロードする必要があります。
library(randomForest)
## randomForest 4.6-10
## Type rfNews() to see new features/changes/bug fixes.
今、1000 個の木からなるランダムフォレストモデルを作成したいとし、使用するコンピュータには 4 つのコアがあるとします。この問題を 4 つに分割するには、randomForest 関数を ntree=250 にして 4 回実行します。もちろん、各 randomForest の実行結果を結合する必要があります。そのためには、randomForest パッケージの combine 関数がちょうど使えます。
それではやってみたいと思いますが、まずは逐次実行を行います：
rf <- foreach(ntree=rep(250, 4), .combine=combine) %do%
randomForest(x, y, ntree=ntree)
rf
##
## Call:
## randomForest(x = x, y = y, ntree = ntree)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 2
これを並列実行するには、%do% を変更する必要があります。さらに、.package オプションを使って、この式を実行するのに randomForest パッケージが必要であることを伝えます。並列バージョンは次のようになります：
rf <- foreach(ntree=rep(250, 4), .combine=combine, .packages='randomForest') %dopar%
randomForest(x, y, ntree=ntree)
rf
##
## Call:
## randomForest(x = x, y = y, ntree = ntree)
## Type of random forest: classification
## Number of trees: 1000
## No. of variables tried at each split: 2
どこかで並列計算をやった経験のある人は、x と y を扱うために、特別なことをする必要がないことを不思議に思うかもしれません。なぜこれで大丈夫なのかというと、%dopar% 関数は、これらの変数が参照されていることと、現在の環境中で定義されていることを知っているからです。したがって、%dopar% はこれらの変数を並列実行ワーカーに自動的にエクスポートし、foreach で実行される表現式の評価に使います。これは、現在の環境中で定義された関数に対しても同じことが言えます。しかし、今回の場合、使用する関数はパッケージで定義されていたため、.package オプションを使ってワーカーにロードするパッケージを指定する必要がありました。
5.2 並列 apply
次に、R の標準関数である apply の並列バージョンの作り方を見てみましょう。apply 関数は R によって書かれており、たった 100 行程度のコードですが、一目で理解するには少し難しいコードです。しかし、このコードの本質は 2 つの for ループに集約されます。その少し複雑な部分は次のようになります：
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
ans <- vector("list", d2)
for(i in 1:d2) {
tmp <- FUN(array(newX[,i], d.call, dn.call), ...)
if(!is.null(tmp)) ans[[i]] <- tmp
}
ans
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
## [[1]]
## [1] 2.5
##
## [[2]]
## [1] 6.5
##
## [[3]]
## [1] 10.5
##
## [[4]]
## [1] 14.5
これはその部分を関数として書き直したものですが、そうしないと、間違った場所で "..." を使っていると R に文句を言われるでしょう。
これは、foreach を使うと次のように書けます：
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
foreach(i=1:d2) %dopar%
FUN(array(newX[,i], d.call, dn.call), ...)
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
## [[1]]
## [1] 2.5
##
## [[2]]
## [1] 6.5
##
## [[3]]
## [1] 10.5
##
## [[4]]
## [1] 14.5
ただし、この方法では、各並列実行ワーカーすべてに配列 newX を丸ごと渡すことになります。各反復計算に必要となるのは配列の一行だけなので、できれば、この余分なデータ通信は避けたいところです。
この問題を解決する方法の一つに、行列を列ごとに反復させるイテレータを使うという方法があります。
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
foreach(x=iter(newX, by='col')) %dopar%
FUN(array(x, d.call, dn.call), ...)
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
## [[1]]
## [1] 2.5
##
## [[2]]
## [1] 6.5
##
## [[3]]
## [1] 10.5
##
## [[4]]
## [1] 14.5
この方法では、並列実行ワーカーに行列の列を一列ずつ渡しています。しかし、もっと大きなかたまりで渡した方が効率的になります。そのためには、行列の複数の列を返すイテレータ、iblkcol 関数を使います。この場合、R の表現式は、部分行列の各列に対して与えられた関数を実行する必要があります。
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
foreach(x=iblkcol(newX, 3), .combine='c', .packages='foreach') %dopar% {
foreach(i=1:ncol(x)) %do% FUN(array(x[,i], d.call, dn.call), ...)
}
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
## [[1]]
## [1] 2.5
##
## [[2]]
## [1] 6.5
##
## [[3]]
## [1] 10.5
##
## [[4]]
## [1] 14.5
部分行列の各列に関数を適用するために、%dopar% の中で %do% を使っていることに注意して下さい。ここで %do% を使うことは、イテレータが行列 x の添え字になるという点で理に適っています。なぜなら、%do% は %dopar% とは違って x をコピーする必要はないからです。
6 リスト内包表記
プログラミング言語 Python に慣れている人は、foreach パッケージが Python のリスト内包表記とそう違わないものを提供していることに気づくかもしれません。事実、foreach パッケージは when という関数を持っていて、これは起こるべきでない式の評価を防ぐことができます。これは、Python のリスト内包表記の "if" 節と非常によく似ています。例えば、次のように when を使用することで、イテレータから負の値を除外することができます。
x <- foreach(a=irnorm(1, count=10), .combine='c') %:% when(a >= 0) %do% sqrt(a)
x
## [1] 0.6357139 0.7352845 0.9264606 0.5012398 1.0815130 1.1732876
このトピックについて言うべきことはそんなにありませんが、when を使った foreach が、シンプルなクイックソート関数を古典的な Haskell のやり方で書くのに使えるというのを隠してはおけません。
qsort <- function(x) {
n <- length(x)
if (n == 0) {
x
} else {
p <- sample(n, 1)
smaller <- foreach(y=x[-p], .combine=c) %:% when(y <= x[p]) %do% y
larger <- foreach(y=x[-p], .combine=c) %:% when(y > x[p]) %do% y
c(qsort(smaller), x[p], qsort(larger))
}
}
qsort(runif(12))
## [1] 0.04198463 0.05217435 0.07720299 0.09653661 0.34985193 0.36509786
## [7] 0.66330564 0.86476895 0.90657146 0.94659762 0.99242779 0.99457254
これを R 標準の sort 関数よりお勧めするというわけではありません。しかし、これは foreach を使った本当に面白い例です。
7 まとめ
並列計算の多くは、次の 3 つのことをやっています：問題をいくつかの部分に分割する、分割した部分を並列に実行する、結果を一つに結合する。foreach パッケージを使う場合、イテレータは問題の分割を助け、%dopar% 関数は分割された部分を並列に実行し、.combine 関数は結果を一つに結合します。この文書では、foreach パッケージを使うことで、極めて簡単に並列化できるというシンプルな例を示しました。また、もっと複雑な問題の解決方法について、いくつかアイデアを与えました。ただし、このパッケージは非常に新しいものであり、並列計算を行うためのより強力なシステムにするための取り組みを引き続き行う必要があります。
※訳注
1. 並列化バックエンド
この記事では、コードを並列実行するには、%do% を %dopar% に変えるだけと書いていますが、実際はバックエンドの設定を行う必要があります。
詳しくは
が参考になります。
また、pforeach パッケージを使うと、バックエンドの設定を自動で行ってくれるため便利です。
2. iblkcol
関数 iblkcol はユーザ定義関数なので 5 章の最後のコードは実行できません。
iblkcol の定義は次のファイルに載っています。
system.file("doc/foreach.R", package = "foreach")
iblkcol を使う代わりに、iter 関数のchunksize 引数を指定することで、同じことが実現可能です。
applyKernel <- function(newX, FUN, d2, d.call, dn.call=NULL, ...) {
foreach(x=iter(newX, by='col', chunksize=round(nrow(newX)/3)), .combine='c', .packages='foreach') %dopar% {
foreach(i=1:ncol(x)) %do% FUN(array(x[,i], d.call, dn.call), ...)
}
}
applyKernel(matrix(1:16, 4), mean, 4, 4)
3. foreach のオプション
foreach 関数は様々なオプションを取ることができ、これを使うとかなり便利になります。
foreach のオプションについては、下記記事がよくまとまっています。
-
もちろん、
sqrtはベクトル化された関数なので、実際はこのようなことをする必要はありません。しかし、あとで見るように、foreachはベクトル化された関数と一緒に活用できます。 ↩