Tsukuba.R#4での発表資料

上げておきます。

Rを高速化するための10の方法

id:syou6162

自己紹介

吉田康久
- 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")

Quartz 2 [*]

その他の例

行列の対角要素を取ってくる
- 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))

Quartz 2 [*]

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
  - 50倍くらい早くなったり

疑問「じゃあ、RじゃなくってCだけでやればよくね?」

Rでやる意義
- 可視化がしやすいので、プログラムの間違いに気づきやすい
- MCMCとかだと、色々判定が必要
  - 系列相関がないか
  - 収束しているか
- 信頼性の高い数学的関数

で、結局どのくらい早くなったのさ?

「Cレベルを使った高速化」でのまとめ

簡単なものはポインタを使う程度でできる
- C始めて一ヶ月の俺でもできるよ!!
RのオブジェクトはCだとSEXPオブジェクトとして扱われる
ドキュメントがしっかしているので、あとは皆がexampleを増やすだけ!!

並列化による高速化

拡張パッケージによる並列処理
snowパッケージ

snowパッケージの特徴

今ある関数をそのまま使い回せる

並列化の例

中心極限定理
簡単に言えば
- 何かしらの乱数があって
- それをたくさん持ってくる
- で、平均を取ると
- 正規分布っぽくなるよ

Tsukuba.R#1でちょっとやった

Rを使った統計勉強法【中心極限定理編】 - Seeking for my unique color.

これの並列化をしたい

1回目でn個のサンプルを取ってきて、その平均を出す
2回目でn個のサンプルを取ってきて、その平均を出す
…
m回目でn個のサンプルを取ってきて、その平均を出す

注意

なんでも並列化できるわけじゃない
- 互いに依存しあうような処理は並列化できない
多少不安定な関数もある

3つの方法

PVM(Parallel Virtual Machine)クラスタ
MPIクラスタ
ソケットクラスタ
- お手軽

ここの焼き直し＞＜

snowパッケージを使って、並列処理をやってみた - Seeking for my unique color.

# 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))

Quartz 2 [*]

分位点プロット

どうやら正規分布に従っていると言えそうだ!!

qqnorm(means)
qqline(means)

Quartz 2 [*]

欲張ってもだめ＞＜

コア数以上を指定しても早くならない
本当に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