上げておきます。
Rを高速化するための10の方法
自己紹介
- 吉田康久
- id:syou6162
- syouでおk
- Tsukuba大学の4年生
最近の出来事
- 卒研おわた
- ノンパラメトリック回帰
自分でRのパッケージを作ってみた
Google summer of code
- Rのプロジェクトがあるので応募しようと思っている
ここから本題
- Rを高速化するための10の方法
注意
- 割とアドバンスな内容を話すかも
- だけど、R初心者の人にも役に立つような話も折り混ぜていくので聞いてやってくだしあ><
きっかけ
- 卒論が理論&シミュレーション
- Macbookをがんがん回す
- まともに回すと3日かかる><
3つのステップ
- Rレベルでの高速化
- Cレベルを使った高速化
- 並列化による高速化
まとめ
- Rレベルでの高速化
- 50倍とか1000倍近くの高速化。
- Cレベルを使った高速化
- 50倍とか1000倍近くの高速化。
- 並列化による高速化
- (総時間)/コア数
Rレベルでの高速化
- 「意味のあるときは常に、意味の無いときも常にベクトル化を心がけよ」
matzも言及
さらに、最近データ処理や統計処理に適した言語としてにわかに注目を集めている「R」など、APLのように配列やベクトルの処理を得意とする言語もある。
未来の言語は「APL」? Rubyのまつもと氏が講演 − @IT
Rでのベクトルおさらい
- JavaみたいにInt型の整数とかいう変数は存在しない
- スカラーはない
- ベクトルが基本的なものになっている
> 1:10 [1] 1 2 3 4 5 6 7 8 9 10 > is.vector(1:10) [1] TRUE > 1 [1] 1 > is.vector(1) [1] TRUE
配列
- 配列(Array)とベクトルは別もの
> array(1:12,c(2,3,2)) , , 1 [,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 , , 2 [,1] [,2] [,3] [1,] 7 9 11 [2,] 8 10 12
行列
- "多次元配列"とも"配列の配列"違う
- 行列専用のデータ型(matrix)が用意されている
「意味のあるときは常に、意味の無いときも常にベクトル化を心がけよ」
- RjpWikiより
- 個別に作るより、"まとめて"作る
- 乱数生成もベクトル化された関数の一つ
> n <- 10000 > tmp <- c() > > system.time(for(i in 1:n){ + tmp[i] <- rnorm(1) + }) ユーザ システム 経過 0.422 0.632 1.065 > system.time(rnorm(n)) ユーザ システム 経過 0.001 0.001 0.002
ベクトル化して高速化されるものの典型例
- outer関数
- 2つのベクトルを引数に取り、そのベクトルの要素全ての組み合わせを関数に投げる
例
- 九九の表もすぐにできるよ><
> outer(1:9,1:9,"*") [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1 2 3 4 5 6 7 8 9 [2,] 2 4 6 8 10 12 14 16 18 [3,] 3 6 9 12 15 18 21 24 27 [4,] 4 8 12 16 20 24 28 32 36 [5,] 5 10 15 20 25 30 35 40 45 [6,] 6 12 18 24 30 36 42 48 54 [7,] 7 14 21 28 35 42 49 56 63 [8,] 8 16 24 32 40 48 56 64 72 [9,] 9 18 27 36 45 54 63 72 81
同じことを関数っぽくやってみる
> outer(1:9,1:9, + function(x,y){x*y}) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1 2 3 4 5 6 7 8 9 [2,] 2 4 6 8 10 12 14 16 18 [3,] 3 6 9 12 15 18 21 24 27 [4,] 4 8 12 16 20 24 28 32 36 [5,] 5 10 15 20 25 30 35 40 45 [6,] 6 12 18 24 30 36 42 48 54 [7,] 7 14 21 28 35 42 49 56 63 [8,] 8 16 24 32 40 48 56 64 72 [9,] 9 18 27 36 45 54 63 72 81
どんくらい早くなるの?
- 100倍は軽く違う
> n <- 1000 > > mat <- matrix(n*n,n,n) > > system.time(for(i in 1:n){ + for(j in 1:n){ + mat[i,j] <- i*j + } + }) ユーザ システム 経過 4.668 0.034 4.840 > system.time(outer(1:n,1:n,"*")) ユーザ システム 経過 0.010 0.012 0.023
outer関数の使いどころ
- 3Dのグラフィックスとかでよく使う
x <- seq(-3,3,length.out=30) y <- seq(-3,3,length.out=30) persp(x,y,outer(x,y,function(i,j){ exp(-(i^2+j^2))}), theta=320,phi=20,col=rainbow(50),ticktype="detailed", xlab="x",ylab="y",zlab="density")
その他の例
- 行列の対角要素を取ってくる
- traceの計算とかに使いそう
> n <- 1000 > x <- 1:n > y <- 1:n > mat <- matrix(n*n,n,n) > > system.time(diag(mat)) ユーザ システム 経過 0.007 0.008 0.016 > > d <- c() > system.time(for(i in x){ + for(j in y){ + if(i==j){ + d[i] <- mat[i,i] + } + } + }) ユーザ システム 経過 1.439 0.010 1.490
その他ベクトル化で早くなる関数の例
- 累積和を取るcumsum関数
x <- 50000:1 plot(cumsum(x))
for回したときとの違い
- やたら早い
> x <- 50000:1 > cum.sum.x <- c() > > sum <- 0 > system.time( + for(i in seq(length(x))){ + sum <- sum + x[i] + cum.sum.x[i] <- sum + } + ) ユーザ システム 経過 10.668 9.750 20.455 > system.time(cumsum(x)) ユーザ システム 経過 0.001 0.000 0.001
その他のベクトル単位で計算できる関数
- cumprod
- any
- all
- kronecker
- cummax
- cummin
- diff
- which
- などなど。。。
「Rレベルでの高速化」でのまとめ
- ベクトル化された関数は相当早い
- CやJavaとは違う志向(ベクトル、ベクトル、ベクトル!!)なのでRの流儀に合わせたほうが早い
- ベクトル化された関数を覚えよう
ところで
- ベクトル化すると早いことは分かった
- でも、それ以外にも早いRの関数がある
- ToDo:誰かになんでか聞いてみる
lapplyの中身
- .Internalとか書いてある
> lapply function (X, FUN, ...) { FUN <- match.fun(FUN) if (!is.vector(X) || is.object(X)) X <- as.list(X) .Internal(lapply(X, FUN)) } <environment: namespace:base>
時系列モデル(arima)の中身の抜粋
- .Callとか書いてある
arimaSS <- function(y, mod) { .Call(R_ARIMA_Like, y, mod$phi, mod$theta, mod$Delta, mod$a, mod$P, mod$Pn, as.integer(0), TRUE) }
何が起っているのか?
- RレベルからCの関数を呼び出している!!
- だから早い
- 自分でも書けないのか→書けます→書こう!!
Cレベルを使った高速化
- 拡張できる言語
- C
- C++
- Fortran
- 今回は主にCでの例を紹介
RからCを使う
- 簡単な例から
- ベクトルの内積を計算させる
> sum(1:10 * 10:1) [1] 220
Cのコード
- 引数はポインタで渡す
- 返り値はvoid
- 結果の変数も引数の中で準備
#include <R.h> #include <Rinternals.h> void inner_product(double *x, double *y, int *n, double *result) { int i; double sum = 0; for (i = 0; i < *n; i++){ sum += x[i] * y[i]; } *result = sum; }
コンパイル
- includeとか勝手にしてくれる
/tmp% R CMD SHLIB fuga.c gcc -arch i386 -std=gnu99 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/i386 -I/usr/local/include -fPIC -g -O2 -c fuga.c -o fuga.o gcc -arch i386 -std=gnu99 -dynamiclib -Wl,-headerpad_max_install_names -mmacosx-version-min=10.4 -undefined dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/lib -o fuga.so fuga.o -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation ld: warning, duplicate dylib /usr/local/lib/libgcc_s.1.dylib
Rから呼び出す
- dyn.loadでロード
- それぞれ型を合わせてやる
- resultのところに注意
dyn.load("/tmp/fuga.so") inner_product <- function(x,y){ .C("inner_product", as.double(x), as.double(y), as.integer(length(x)), result=0)$result }
すると…
- 簡単に呼び出すことができる!!
> inner_product(1:10,10:1) [1] 220
デバッグ用の関数
- cのprintfではなくて、Rprintfを使うべし
Rprintf("%f\t%f\t%f\n",x[i],y[i],sum);
さっきまでのは
- ".C"を使った例
- ポインタを渡していた
- 今度はCの中でRのオブジェクトをいじって遊んでみる
lapply関数をCで書いてみる
- ちょっと待て、lapply関数って(ry
- example(lapply)より
- listの要素ごとに関数を適用してくれる
> x <- list(a = 1:10, beta = exp(-3:3), logic = c(TRUE,FALSE,FALSE,TRUE)) > lapply(x,mean) $a [1] 5.5 $beta [1] 4.535125 $logic [1] 0.5 > lapply(x,summary) $a Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 3.25 5.50 5.50 7.75 10.00 $beta Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04979 0.25160 1.00000 4.53500 5.05400 20.09000 $logic Mode FALSE TRUE NA's logical 2 2 0
これの簡易版を作る
- 色々考えることがある
- CでRのlistはどう扱えばいいのか?
- CでRの関数はどう扱えばいいのか?
関数定義
- 引数、戻り値は全て"SEXP"
- "S Expression"だからね><
- Rinternal.hに定義がしてある
- rhoは"環境"(environment)オブジェクト
- eval関数の関係で必要
SEXP mylapply(SEXP list, SEXP fn, SEXP rho)
結果を入れておくオブジェクトを用意
- allocVectorでRのベクトルを用意
- (引数でないSEXPのオブジェクトは)PROTECTという操作をしてやらないといけない
SEXP ans; PROTECT(ans = allocVector(VECSXP, n));
PROTECTマクロ
もし C コード中に R オブジェクトを作ったら、オブジェクトへのポインタに PROTECT マクロを用い R にオブジェクトを使用することを知らせる必要がある。これは R にそのオブジェクトが使用中であることを告げ、そのためそれが破壊されることはなくなる。
Writing R Extensions
使った後は
- UNPROTECTしてあげる
- PROTECTしたオブジェクトの総数が引数
UNPROTECT(2);
listを処理する中身
- Rの先祖はScheme
for(i = 0; i < n; i++) { SETCADR(R_fcall, VECTOR_ELT(list, i)); // S式の肝。(+ 1 2)というような形のリストにする SET_VECTOR_ELT(ans, i, eval(R_fcall, rho)); // ans[i]にeval(R_fcall, rho)の結果を格納 }
ちなみに
- Rinternal.hの中身
- どう見てもlispです、本当に(ry
#define CAR(e) ((e)->u.listsxp.carval) #define CDR(e) ((e)->u.listsxp.cdrval) #define CAAR(e) CAR(CAR(e)) #define CDAR(e) CDR(CAR(e)) #define CADR(e) CAR(CDR(e)) #define CDDR(e) CDR(CDR(e)) #define CADDR(e) CAR(CDR(CDR(e))) #define CADDDR(e) CAR(CDR(CDR(CDR(e)))) #define CAD4R(e) CAR(CDR(CDR(CDR(CDR(e)))))
最後にまたlistの形式に整えていく
- 属性をset
setAttrib(ans, R_NamesSymbol, getAttrib(list, R_NamesSymbol));
全体像
#include <R.h> #include <Rinternals.h> SEXP mylapply(SEXP list, SEXP fn, SEXP rho) { // evalとかの関係で環境も引数に渡してあげないといけない int i, n = length(list); SEXP R_fcall, ans; PROTECT(R_fcall = lang2(fn, R_NilValue)); // 2つのSEXPをconsして、そのTYPEOFをLANGSXPにする PROTECT(ans = allocVector(VECSXP, n)); for(i = 0; i < n; i++) { SETCADR(R_fcall, VECTOR_ELT(list, i)); // S式の肝。(+ 1 2)というような形のリストにする SET_VECTOR_ELT(ans, i, eval(R_fcall, rho)); // ans[i]にeval(R_fcall, rho)の結果を格納 } setAttrib(ans, R_NamesSymbol, getAttrib(list, R_NamesSymbol)); UNPROTECT(2); return(ans); }
実行してみる
- 簡易版lapplyができた!!
> a <- list(a = 1:5, b = c(TRUE,FALSE)) > dyn.load("/tmp/fuga.so") > .Call("mylapply", a, summary, new.env()) $a Min. 1st Qu. Median Mean 3rd Qu. Max. 1 2 3 3 4 5 $b Mode FALSE TRUE NA's logical 1 1 0
SEXPオブジェクトを扱う時のprint
- PrintValueでおk
PrintValue(R_fcall);
CからRを使う
- RのAPI
どういう時に使うか
- 数値解析サブルーチン
- 数学関数
- 乱数生成
例えば
- Writing R Extensionsのp49を参考
分布関数名 | サブルーチン | パラメータ |
---|---|---|
beta | beta | a, b |
non-central beta | nbeta | a, b, lambda |
binomial | binom | n, p |
Cauchy | cauchy | location, scale |
chi-squared | chisq | df |
non-central | chi-squared | nchisq df, lambda |
exponential | exp | scale |
F | f | n1, n2 |
non-central F | {p,q}nf | n1, n2, ncp |
gamma | gamma | shape, scale |
geometric | geom | p |
hypergeometric | hyper | NR, NB, n |
logistic | logis | location, scale |
lognormal | lnorm | logmean, logsd |
negative binomial | nbinom | n, p |
normal | norm | mu, sigma |
Poisson | pois | lambda |
Student's | t | t n |
non-central t | {p,q}nt | df, delta |
Studentized range | {p,q}tukey | rr, cc, df |
uniform | unif | a, b |
Weibull | weibull | shape, scale |
Wilcoxon rank sum | wilcox | m, n |
Wilcoxon signed rank | signrank | n |
より詳細は
- Rmath.hを見てね
- Writing R Extensions
- 20p以降、かな
- An Introduction to the .C Interface to R
Cによる高速化はどのようなものに向いているか?
- 回数をこなさないといけないようなシミュレーション
- モンテカルロ積分
- MCMC
疑問「じゃあ、RじゃなくってCだけでやればよくね?」
- Rでやる意義
- 可視化がしやすいので、プログラムの間違いに気づきやすい
- MCMCとかだと、色々判定が必要
- 系列相関がないか
- 収束しているか
- 信頼性の高い数学的関数
で、結局どのくらい早くなったのさ?
「Cレベルを使った高速化」でのまとめ
- 簡単なものはポインタを使う程度でできる
- C始めて一ヶ月の俺でもできるよ!!
- RのオブジェクトはCだとSEXPオブジェクトとして扱われる
- ドキュメントがしっかしているので、あとは皆がexampleを増やすだけ!!
並列化による高速化
- 拡張パッケージによる並列処理
- snowパッケージ
snowパッケージの特徴
- 今ある関数をそのまま使い回せる
並列化の例
- 中心極限定理
- 簡単に言えば
- 何かしらの乱数があって
- それをたくさん持ってくる
- で、平均を取ると
- 正規分布っぽくなるよ
Tsukuba.R#1でちょっとやった
これの並列化をしたい
- 1回目でn個のサンプルを取ってきて、その平均を出す
- 2回目でn個のサンプルを取ってきて、その平均を出す
- …
- m回目でn個のサンプルを取ってきて、その平均を出す
注意
- なんでも並列化できるわけじゃない
- 互いに依存しあうような処理は並列化できない
- 多少不安定な関数もある
3つの方法
- PVM(Parallel Virtual Machine)クラスタ
- MPIクラスタ
- ソケットクラスタ
- お手軽
ここの焼き直し><
# install.packages("snow") # library(snow) cl <- makeCluster(2, type = "SOCK")
準備
- 必要な関数を準備
- クラスタの内部の計算で使うオブジェクトは、clusterExportで教えてあげないといけない
> rnorm.mean <- function(n){ + mean(rnorm(n)) + } > n <- 1000 > m <- 10000 > clusterExport(cl,"rnorm.mean") > clusterExport(cl,"n")
並列化しない場合
> system.time(sapply(1:m,function(x){rnorm.mean(n)})) ユーザ システム 経過 2.096 0.019 2.232
並列化した場合
- 倍くらい早くなっている!!
- 2コアマシン
> system.time(parSapply(cl,1:m,function(x){rnorm.mean(n)})) ユーザ システム 経過 0.033 0.005 1.310
ちなみに本当に正規分布になっていたのか?
- 密度トレイスを書いてみる
- ヒストグラムのようなもの
means <- parSapply(cl,1:m,function(x){rnorm.mean(n)}) plot(density(means))
分位点プロット
- どうやら正規分布に従っていると言えそうだ!!
qqnorm(means) qqline(means)
欲張ってもだめ><
- コア数以上を指定しても早くならない
- 本当に4コアマシンだったら1/4の実行時間になる
> cl <- makeCluster(2, type = "SOCK") > clusterExport(cl,"rnorm.mean") > clusterExport(cl,"n") > system.time(parSapply(cl,1:m,function(x){rnorm.mean(n)})) ユーザ システム 経過 0.033 0.006 1.393 > > cl <- makeCluster(4, type = "SOCK") > clusterExport(cl,"rnorm.mean") > clusterExport(cl,"n") > system.time(parSapply(cl,1:m,function(x){rnorm.mean(n)})) ユーザ システム 経過 0.031 0.005 1.350
お願い
- もっと多いコア数のマシンとか、クラスタとかスパコンとか使える人
- cl <- makeCluster(30, type = "SOCK")とかで実験して、比較実験してくだしあ><
「並列化による高速化」でのまとめ
- snowパッケージ超便利
- 2コアくらいだとせいぜい倍速
- 計算機環境整っている人の実験求む!!
発表全体のまとめ
- Rレベルでの高速化
- ベクトル、ベクトル、ベクトル!!
- Cレベルを使った高速化
- ポインタとSEXP
- 並列化による高速化
- snowパッケージを使って倍速くらいになった
余った時間で
- 自分のパッケージ作ったりとかどんどんやるといいと思います><
- おわり