須通り_統計_複数の引数を取るRの関数を反復・並列化するmapplyとclusterMap

ホーム | 統計 Top | 複数の引数を取るRの関数を反復・並列化するmapplyとclusterMap

（2015年2月20日初版）なぜかclusterMap関数の日本語情報が極端に少ないので、ちょっと念入りに紹介してみる。

はじめに

数値シミュレーションにおいて、複数パラメータを網羅的に動かしつつ、それらを引数に取るような関数を適用したいことがある。最も原始的なやり方はfor文をネストして

（変数iとjについて各々1から10まで動かした全組み合わせを試す）

for(i in 1:10) { 
    for(j in 1:10) { 
        ここに処理内容 
    } 
}

などとするのだが、Rの場合は残念ながらfor文の実行速度があまり速くないという欠点がある。

では、繰り返し計算の速度をいかにして向上させるのか？幾つかの原則が広く知られている。

ベクトル化、行列化
applyファミリー関数の使用
マルチコア化（マシン内の複数のCPUで並列計算）
C++コードを混ぜる（Rcpp）

他にも方法はあるが（GPGPUとか分散コンピューティングとか）、とりあえずハード面での追加投資が不要な奴はこれくらい。入門のハードルとしては上に挙げたものほど容易である。本ページでは、これらの手法のうち、前三者の導入を手短に解説する（Rcppについては別口で扱う）。

はじめに断っておくと、実行速度がボトルネックにならない事例では、for文でごり押しすればいい。しばしば「Rでfor文を使うべきでない」とまで言う原理主義者も見受けられるが、無理に難しい構文を書き、数ヶ月後にそれを解読する人間側の労力に比べれば、コンピュータの計算リソースはタダ同然に安いものだ。

準備

「複数の引数を取るRの関数」の簡単な例として、台形の面積の公式を以下のように定義しておこう（下をRコンソールにコピペして関数定義）。

# 台形の面積公式。a, b は上底と下底の長さ、h は高さ
trapezoid <- function(a, b, h) {
  return( (a+b)*h/2 )
}

関数をベクトル化する

それほど特別なことをやるわけではない。一般的なRの四則演算子は、複数のベクトルや行列を（大きさが同じもの同士であれば）食べさせたとき、「対応する要素同士の計算」として解釈した結果を返すよう実装されている。たとえば c(1,2,3)+c(5,6,7) は、i番目の要素について足し合わせたベクトルである c(6,8,10)を返す。同じことが「ベクトルや行列を扱える演算子を、直列に組み合わせて作った関数」についても言え、実は上で定義した台形の面積公式も既にベクトル化されている。

> trapezoid(1,3,5) # スカラーで引数セットを与えて計算する
[1] 10
> trapezoid(2,3,2) # こちらもスカラーで
[1] 5
> trapezoid(c(1,2),c(3,3),c(5,2)) # ベクトルで引数セットを与えても計算できる
[1] 10  5

しかしベクトル化が一筋縄で行かないこともある。最たるものが、関数に条件分岐が含まれる場合だ。

# 台形の面積公式を、高さが正であるときだけ計算させるようにした関数（ベクトル化されていない例）
trapezoid.modified <- function(a, b, h) {
  if (h > 0){
    return( (a+b)*h/2 )
  } else {
    return(0)
  }
}

> trapezoid.modified(1,3,5) # スカラーを与えて計算した例
[1] 10
> trapezoid.modified(2,3,-2) # 台形の高さが負数だとゼロを返すようになっている
[1] 0

> trapezoid.modified(c(1,2),c(3,3),c(5,-2)) # ベクトル化できておらずエラーとなる例
[1] 10 -5
 警告メッセージ: 
In if (h > 0) { :  条件が長さが 2 以上なので、最初の 1 つだけが使われます

じゃあどうするのさ。上の例だと、if()の代わりにifelse()関数を使うという回避技が存在する。残念ながらこいつで書き直すと構文が変わってしまうため、常に通用する手口ではないが。

# 台形の面積公式を、高さが正であるときだけ計算させるようにした関数（こちらは正しく動く）
trapezoid.modified2 <- function(a, b, h) {
  h <- ifelse(h > 0, h, 0)
  return((a+b)*h/2)
}

> trapezoid.modified2(c(1,2),c(3,3),c(5,-2))
[1] 10  0

mapplyでループを隠蔽する

上に書いたベクトル化は、しばしば関数の中身を書き直す必要があるし、そもそもベクトル化された関数を見つけねばならず、あまり普遍性が高いとは言えない。そこで次の選択肢が、applyファミリー関数の一種たるmapplyだ。実際にやっていることはリスト作成の隠蔽なので、厳密には「高速化手法」ではないが、Rプログラミングの思想を理解する上でも触れておくべきテクニックだろう。

# 公式ヘルプより mapplyの使い方
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/mapply.html
mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE,
       USE.NAMES = TRUE)

FUN   
計算に用いる関数の名を、FUN=名称文字列 などとして指定する。

...   
関数FUNに渡したい（計算毎に異なる値を取るような）追加の引数を、ベクトルもしくはリスト形式で指定する。引数が複数種ある場合はその都度コンマで区切るか、リストにまとめる。

MoreArgs  
FUNに渡したいが、値そのものは一連の計算を通じて固定で構わない引数は、ここで指定できる（複数種でも可）。リスト形式にまとめる必要がある。たとえば, MoreArgs=list(...), といった具合。

SIMPLIFY  
生の計算結果はリストにまとめられた形で返ってくるが、これを行列やベクトルといったより単純な形式に変換するかどうか。
SIMPLIFY=TRUEとすることは、通常のapplyファミリーにおいて、lapplyの代わりにsapplyを使うのに等しい。この場合、たとえば関数FUNが本来持つ戻り値が1つであれば、mapplyの結果は1本のベクトルにまとめられる。

USE.NAMES   
結果オブジェクトの行名を、引数の組から自動的に取ってくるか否かを指定する論理値。

# 早速使ってみよう
upper.length <- seq(1, 300, by=3)
lower.length <- seq(1, 200, by=2)
result1 <- mapply(FUN=trapezoid, 
                  a=upper.length, 
                  b=lower.length, 
                  MoreArgs=list(h=30), 
                  SIMPLIFY=TRUE, 
                  USE.NAMES=TRUE ) # 律儀に全パラメータの名前を指定したが、必須ではない
summary(result1)

> summary(result1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     30    1886    3742    3742    5599    7455

上記の例だと、上辺の長さはc(1, 4, 7, ..., 300)、下辺はc(1, 3, 5, ..., 200)という、それぞれ長さ100のベクトルで与えている。一方、台形の高さは固定でいいことにして、h=30という定数をMoreArgsに放り込むことで解決。このときmapplyによる計算は都合100回ほど繰り返されており、最初の計算に代入した寸法は c(a, b, h)=c(1, 1, 30)、2回目はc(4, 3, 30) ...（以下略）となっている。

発展的な例として、上辺のパラメータと下辺のパラメータが所与のベクトルとしてそれぞれ与えられており、さらにその間で全組み合わせを網羅的に試したい時は、expand.grid()関数をまず用いてパラメータセットをデータフレーム化し、そいつの各列を指定してmapplyの引数にぶち込む。感度分析においては超強力な手法だが、調子に乗って変数の数を増やしていくと計算時間が乗数で増えることに注意。上の例で言えば、総計算回数は100×100で1万回となる。あるいは2変数なら、outer()を使う手もある。

clusterMapで計算をマルチスレッド化する

最近のパソコンはCPUにプロセッサコアが複数備わっており（マルチコア）、我々の快適な文明生活に寄与するところ大なり。その恩恵に与るべく、parallelパッケージに含まれるclusterMapを使って計算を並列化する。

# 公式ヘルプより clusterMapの使い方（勝手に加筆あり）
# https://stat.ethz.ch/R-manual/R-devel/library/parallel/html/clusterApply.html
clusterMap(cl = NULL, fun, ..., MoreArgs = NULL, RECYCLE = TRUE,
           SIMPLIFY = FALSE, USE.NAMES = TRUE,
           .scheduling = c("static", "dynamic"))

cl  
本パッケージ（library(parallel)のこと）もしくは"snow"パッケージによって作成されるクラスターオブジェクトの名前。

fun
計算に用いる関数の名を、fun="名称文字列"などとして指定する。clusterApplyファミリーの関数には小文字でfunとして与えるものと、大文字のFUNとして与えるものとが混在しているので注意。clusterMapでは小文字である。

...   
関数funに渡したい（計算毎に異なる値を取るような）追加の引数を、ベクトルもしくはリスト形式で指定する。引数が複数種ある場合はその都度コンマで区切るか、リストにまとめる。
なお原文に、"beware of partial matching to earlier arguments" という注意書きがある。つまり引数の名前が部分一致した時点で先走ってfunに渡されてしまうのかな？

MoreArgs  
funに渡したいが、値そのものは固定で構わない引数は、ここで指定できる（複数種でも可）。リスト形式にまとめる必要がある。たとえば, MoreArgs=list(...), といった具合。

RECYCLE   
論理値。ある引数に与えたベクトルが他のものより短かったとき、再帰的にベクトルの冒頭に戻って、足りないデータを適宜埋めていいですか？うむ、よきにはからえ。

SIMPLIFY  
生の計算結果はリストにまとめられた形で返ってくるが、これを行列やベクトルといったより単純な形式に変換するかどうか。これは利便性に応じて使い分けるといい。

USE.NAMES   
結果オブジェクトの行名を、引数の組から自動的に取ってくるか否かを指定する論理値。通常はTRUEのままでいいと思うゾ。

.scheduling   
各コアへの仕事の割り振りを最初に固定するか、動的なload-balancingを行うか（下の記述を参照）。

# 早速使ってみよう
# 並列化の準備
library(parallel) # R 2.14.0以上を要するので注意
cl<-makeCluster(detectCores()) # サブプロセスを cl という名で作っている。

# パラメータセットのネタを作っておく
upper.length <- seq(1, 300, by=1)
lower.length <- seq(1, 200, by=1)

# 折角だからexpand.grid()で、パラメータセットを格納したデータフレームを作って計算ネタにする
om <- expand.grid( upper.length=upper.length, lower.length=lower.length)
print(nrow(om)) # 計算回数＝データフレーム列数を実行前に確認しよう。増やしすぎないよう注意

# 下の行で並列計算を実行。重い（2.5 GHz Core i5搭載のMBPで13秒かかった）ので注意。
# 引数としてデータフレームomの列を指定しているのがこの処理の肝である。
system.time( result.raw <- clusterMap(cl=cl, fun=trapezoid, 
                                      a=om$upper.length, b=om$lower.length, 
                                      MoreArgs=list(h=30), 
                                      SIMPLIFY=F, .scheduling = "dynamic" ) )

# 計算が終わったら結果を取り出して、omの新しい行に書き込む。
om$result <-  sapply( result.raw, FUN=function(x) { x[[1]] }) 

# SIMPLIFY=Fとして計算させたときの、リストからの結果取り出し。
# [[]] 二重カッコで囲むことにより、リストの要素に直接アクセスできる。

# ガベージコレクション（計算後にやっておくと吉）
gc();gc();
invisible(clusterEvalQ(cl, gc())); invisible(clusterEvalQ(cl, gc()));

stopCluster(cl) ## 1日の最後に並列処理プロセスを停止させる

> summary(om)
  upper.length     lower.length        result    
 Min.   :  1.00   Min.   :  1.00   Min.   :  30  
 1st Qu.: 75.75   1st Qu.: 50.75   1st Qu.:2272  
 Median :150.50   Median :100.50   Median :3765  
 Mean   :150.50   Mean   :100.50   Mean   :3515  
 3rd Qu.:225.25   3rd Qu.:150.25   3rd Qu.:4890  
 Max.   :300.00   Max.   :200.00   Max.   :6000

注意点

結果オブジェクトのリストからの取り出しには慣れが必要。最初は小さなデータで練習すること。
計算の途中経過を知ることができないので、やはり最初に組み合わせを絞って、およその所要時間を計ってから本番に入ろう。
計算に用いる関数は、Rコンソール上でコピペして定義する場合には難しいことは考えなくていいが、外部のテキストファイルに保存しているものを source("source.R") などとして読ませた場合は、各クラスターに自動では定義が渡されないことに注意。その場合は、以下のようにclusterCall()関数を使うと、クラスターに命令を渡せる。

Sys.setlocale(locale="C") # ソースコードに日本語コメントが入っている場合はロケール変更。
invisible(clusterCall(cl, source, "source.R")) # クラスターにソースを読ませる。コアの数だけメッセージが出ると面倒なので、invisibleで出力抑制
Sys.setlocale() # ロケールを戻す

並列化における計算効率とload-balancingのはなし

唐突に話は変わるが並列化の身近な例として、小学校の先生に「算数ドリルの1ページ目から4ページ目までを解きなさい」という宿題を与えられた吉岡君、吉川君、吉田君、吉村君の4人が、それぞれ1ページずつ自宅で計算して翌朝、答えを持ち寄るケースを考える。なお本手法の行く先には、計算効率の最適化という昏い淵が横たわっている。たとえば算数ドリルの後ろのページに行くに従って問題が難しくなり、4ページ目担当の吉村君に過度の負担がかかる＝計算により多くの時間を要するかもしれない。かたや1ページ目担当の吉岡君が簡単な計算をあっさり終えて遊んでいたとしたら、将来的に4人の友情に亀裂を入れるのみならず、甚だ非効率的である。

これを解決するのがload-balancingと呼ばれるテクニックで、たとえば4人がLINEで連絡を取り合って、先に計算を終えた吉岡君が、ドリル4ページ目の問題の一部を肩代わりすることも可能だ。clusterMap関数では .scheduling="dynamic" オプションを指定するのだが、実は有効な実装が未だ為されていないという噂もある。

だったら手動で静的なload-balancingをやってしまえばいい。その一つとして、計算のネタとなるデータフレームの行をシャッフルしてからclusterMapに渡す方法を紹介する。パラメータリストの作成時にはしばしば「小さい数字から大きい数字へ」みたいなベクトルを作って組み合わせるため、そのまま渡すと、計算の重いパラメータ領域（例：大きな数が多い）が特定のコアに集中してしまいがちなのだ。

om <- expand.grid( upper.length=upper.length, lower.length=lower.length) # 上記同様に引数セットを格納したデータフレームを作る。
om <- cbind(LID=c(1:nrow(om)), om ) # データフレームの行番号に相当する、一意のIDを割り振る。
randomize.seed <- sample(nrow(om), replace=FALSE) # 1から「omの行数」までの自然数を、ランダムに入れ替えた乱数ベクトルを発生させる。
om <- om[ order(randomize.seed), ] # 上の乱数を使って、データフレームの行をシャッフルする。

#（ここに計算用コードを挿入）
#（計算が終わったら、omに新しい列を作って結果を代入）

om <- om[ order(om$LID), ] # 計算・代入終了後に、格納しておいたIDを使ってデータフレームの並びを復元

安易に並列化できないケース

インタプリタ言語であるが故に「遅い」Rであっても、様々な並列化手法を駆使してそれなりの計算速度を実現できることは、本稿で述べてきた通り。しかし限界もある。ある計算の結果を次の計算に順次代入していく、純然たる繰り返し計算だ。時間ないし生物の世代についてのシミュレーションがしばしば、このケースに当てはまる。

先ほどの宿題の例ならば、算数ドリルの2ページ目が1ページ目の計算結果を使った発展問題で、吉岡君が深夜に解き終えてからでないと吉川君が仕事に掛かれない場合である。

詳細の説明は本稿の範囲を超えているが、対処原則としてはメモリへのアクセス回数をなるべく減らすとよい。

たとえば、case 1 よりも case 2 の方が関数をネストしている分、一般論としては若干高速である。

# case 1
system.time(
for (i in 1:10000) {
  a <- rnorm(100);
  b <- abs(a);
  c <- sort(b);
  d <- c[50];
} )

# case 2
system.time(
for (i in 1:10000) {
  A <- sort(abs(rnorm(100)))[50];
} )

まあ上の例だと速度は1%程度しか変わらないんですけどね。関数をネストするメリット・デメリットは、どちらかといえば計算途中のデータへのアクセス性を求めるか否かであろう。上 case 1 では、計算の途中結果が a, b, c, d というデータオブジェクト（全て数値ベクトル）として、メモリ上のそれぞれ異なる領域にコピーされる。従って、途中データを保存して後から再利用するには便利だが、メモリ上の番地を行ったり来たりして計算する分、僅かに遅い。また仮に計算途中で物理メモリが枯渇しようものなら、HDDにデータをガリガリと書き込み始め、とんでもなく低速になってしまう。

そして計算に掛かる時間が明らかに数日、数週間のオーダーに達するようになったら、さすがにRのみでは限界がある。そんなときはC++やFORTRANといった、より低級な言語の出番だ。現状のR環境だとRcppを用い、C++で書き直したサブルーチンを呼び出すのがベターな選択だろう。筆者はRcppについては勉強中であり、ある程度地平が見えてきたら次回記事で取り上げようと考えている。