読者です 読者をやめる 読者になる 読者になる

Tsukuba.R#4での発表資料

Tsukuba.R R

上げておきます。

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

自己紹介

  • 吉田康久
  • Tsukuba大学の4年生

最近の出来事

  • 卒研おわた
  • ノンパラメトリック回帰


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

どういう時に使うか

  • 数値解析サブルーチン
  • 数学関数
    • 乱数生成

例えば

分布関数名 サブルーチン パラメータ
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

より詳細は

Cによる高速化はどのようなものに向いているか?

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

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

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

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

並列化による高速化

snowパッケージの特徴

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

並列化の例

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

これの並列化をしたい

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

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 

お願い

  • もっと多いコア数のマシンとか、クラスタとかスパコンとか使える人
  • cl <- makeCluster(30, type = "SOCK")とかで実験して、比較実験してくだしあ><

「並列化による高速化」でのまとめ

  • snowパッケージ超便利
  • 2コアくらいだとせいぜい倍速
  • 計算機環境整っている人の実験求む!!

発表全体のまとめ

  • Rレベルでの高速化
    • ベクトル、ベクトル、ベクトル!!
  • Cレベルを使った高速化
    • ポインタとSEXP
  • 並列化による高速化
    • snowパッケージを使って倍速くらいになった

余った時間で

  • 自分のパッケージ作ったりとかどんどんやるといいと思います><
  • おわり