RのforループをC++で高速化する（Rcppパッケージ）

"時間が掛かるからRでforループを使うな"、"applyファミリーを使え（ほかにsapply, lapplyなど）"とよく言われる。かといって無理に使うと難解なコードになる場合もあるし、せっかく実現しても計測してみたら、むしろforループの方が早かったなんてこともあった（今回の動機）。並列化（複数のCPUコアを用いる）するという手も考えたが、大きなー繋がりのタスクを１回やるだけならよいが、入れ子になっていると呼び出す時間の方が律速になったりもする。

applyファミリーやforeachなどで置換えにくい計算例として、値を毎回更新するようなプログラムがあるだろう。例えば以下のような計算例（"0.9* "さえ無ければapply(x, 2, cumsum)で一発なのだが、sapplyとapplyを組み合わせて書いたらforより遅くなった…）。

NN <- 10^6
x <- matrix(0, NN, 10) 
# NN回分の結果を格納するための入れ物（1回あたりは長さ10のベクトル）
system.time( for(r in 2:(nrow(x))) x[r,] <-  0.9*x[r-1,] + rnorm(10) )
# ひとつ前の結果に0.9を掛け、正規乱数を足していくマルコフ連鎖
#  ユーザ   システム       経過
#    5.457      0.363      5.778

手元の環境では6秒弱掛かった。

いろいろ調べて行き着いたのが Rcpp というパッケージ（要インストール）。RからC++（"しーぷらぷら"と呼ぶみたい、なのでcppなのだろう）で作成した関数を呼び出せるというもの。Cで作られているというパッケージは最近よく見かけるようになったが、自分でカスタマイズできるのはすごい。もっとも、C++を書けるならばですが。それでも部分的にだけでも記述できればだいぶ速くすることができるはず。目的のforループ作成を達成するのに丸一日掛かりましたが、最終目的の計算に掛かる予想時間が数ヶ月だったことを考えれば大幅な時間短縮になりました。

手順として、まず開発環境のセットアップが必要です。Macの場合はXcode（AppStoreからフリーでダウンロード）とX11（ユーティリティ内にあるX11のアイコンをダブルクリックして立ち上げてライセンスにOKすれば大丈夫なはず）。Windows版はすみませんよく分かりません…コマンドプロンプトでC++をコンパイルする必要があると思います。
次に肝心の、C++で記述したコードが必要です。ただし随所に純粋のC++と異なる書き方を要する箇所があります（私はここでつまづきました）。C++で記述したファイルをRの作業ディレクトリに置いて読み込むと、Rの関数として使用できるようになります。





まずは最小限のコード例：（以下、C++のコードは紫にしておきます）
require(Rcpp) # Rcppパッケージを呼び出す（この行はRに直打ち）

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

double testF(double x) {
    return(x*x);
}

紫で書いた6行分をコードエディタなどで書いて test.cpp というファイル名で作業ディレクトリに保存します（ .cpp ってC++ファイルの拡張子なのだと思います）。
C++用コードの初めの3行は決まり文句のようなものです。3行目は用途によって変える場合があるようですが、とりあえず忘れてよさそう。
まずRとの大きな違いはすべての変数型を逐一定義すること。この性質は嫌いじゃないです。むしろ勝手に変数型が変わるのはRが嫌われる点でもありますね。

double は実数の意味、整数なら int 、のように、使う変数の前に半角スペースを置いて書いて定義する必要があります。

testF は、いま定義しようとしている自作関数の名前です。括弧の中で x という変数が与えられた時に { } の処理を行う変数を作成することを宣言しています。

return(x*x); これは x*x の結果を返すことを示しています。もうひとつRと大きな違いとして各行の最後に ;（セミコロン）が入ることです（中括弧 { の後ろ以外）。


次にR上で以下のようにして、この自作関数を呼び出します（初回の呼び出しには数秒掛かるかもしれない）。
sourceCpp("test.cpp")

すると、testF 関数が使用できるようになります。例えば…
testF(1.41421356) # 一夜一夜に人見頃
結果はちゃんと 2 が返ってくるはず。


ちなみにRの上でC++コードを走らせることも可能です。""で挟んで、code=code.testで指定してやるだけです。
sourceCpp(code="
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]

double testF(double x) {
    return(x*x);
} ")

ただし、Rcpp.h以外のファイルを読み込む必要がある場合のやり方が今ひとつ分からず。inlineパッケージのcxxfunctionを使ってできそうなのですが。






次の例では段階を上げて行きます。
rnorm 関数をC++で書いてみましょう。なお、乱数の発生アルゴリズムとしてRの場合はMersenne Twister法が使用されていますが、C++でやろうとすると一工夫必要でした。これについては、詳しい解説とコードも配布されているサイトがありましたのでリンクを貼っておきます。MT.hファイルと内包されるinit_genrand関数、genrand_real3関数はこちらから取得しました。（追記：他にも、匿名さんから頂いたコメントで乱数発生アルゴリズムのまとめサイトがあるとのことです。）

r_normal という名前でC++版rnormを作成します。
下準備としてMt.hと名付けたファイルをRの作業ディレクトリに置きます。下記がRcppから読み込む用のC++コードになります（もはや純粋なC++コードではないので、こう呼びます）。

#include <Rcpp.h>
#include "MT.h"

using namespace Rcpp;
// [[Rcpp::export]]

NumericVector r_normal( int NN, double mu, double sigma ) {  
  NumericVector x(NN);
init_genrand((unsigned)time(NULL));

for(int i=0; i<NN; i++) {
x[i] =  mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );
}
return(x);
}

新たに加わった include "MT.h" で上記のファイルを詠みこみます

#include <stdio.h>
#include "MT.h"

次に、この行ですが、これがクセモノです。
NumericVector r_normal( int NN, double mu, double sigma ) {

NumericVector は返り値がベクトルになるように定義しますという意味です（Rcpp特有の表現：cf. Rcpp Quick reference）。そして r_normal という名前の自作関数を作ろうとします。

int NN, double mu, double sigma は、ベクトルの要素数 NN（整数）、平均が mu（実数）、標準偏差が sigma（実数）の値を用いることを宣言しています。

NumericVector x(NN);
次にこれですが、xという、要素数が NN のベクトル（初期値は全部0になっている）を用意することを意味しています（Rcpp特有の定義のようだ）。

init_genrand((unsigned)time(NULL));
これは乱数のタネの指定、現在時刻を秒精度で使用します。

for(int i=0; i<NN; i++) {

C++版のforループ。整数 i を 0からNN未満まで適用しますという意味。0から始めるのが基本である模様（誤解があるかもしれないけれど、とりあえずこのコードではOK）。

x[i] =  mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );

xのi番目の要素に正規乱数を格納するコード。右辺はごちゃごちゃしていますが、genrand_real3が一様乱数を返す関数です。詳しい説明が参照元にあります。
return(x); で x を返り値で返すようにすれば、r_normal関数が使えるようになります。

rnormと速度比較をしてみます。

NN <- 10^7
system.time(rnorm(n, 0, 1))
#   ユーザ   システム       経過
#    0.677      0.006      0.679

system.time(r_normal(n, 0, 1))
#   ユーザ   システム       経過
#    0.363      0.005      0.367

速度は2倍弱といったところ。もう少し期待したかったですが、それでもベクトル化されているRコードよりもさらに高速です。


一方、Rの関数をC++側に渡すことも可能です。Rの関数をC++で計算させることができます。やり方としては、cppFunction関数を用いてrnormのC++版、rnormCppを定義します。
cppFunction( "
NumericVector rnormCpp(int N, double me, double sd) { 
return(rnorm(N, me, sd)); 
}" )

定義の仕方などはだいたい同じですが、rnormが使われていることに注意です。これはC++の関数ではなく、たしかにRのrnorm関数です。
これを通すと、rnormCppという関数が使えるようになります。気になる速度ですが…

system.time(rnormCpp(10^7, 0, 1))
#   ユーザ   システム       経過
#     0.577      0.019      0.594

Rのrnormより少し速くはなりましたが、C++で１から作るより遅いです。それでも手軽に計算速度を向上させることができるので、何かと試したくなる技です。

ではいよいよ、本題のforループを r_markov 関数として記述してみます（同様にファイル保存します）。

#include <Rcpp.h>
#include "MT.h"

using namespace Rcpp;
// [[Rcpp::export]]

NumericMatrix r_markov( int repl, int length, double mu, double sigma ) {
NumericMatrix x(repl, length);
int i,j;

init_genrand((unsigned)time(NULL));

for(i=1; i<repl; i++) {
for(j=0; j<length; j++) {
x(i,j) =  0.9*x(i-1,j) + mu + sigma*sqrt( -2.0*log(genrand_real3()) ) * sin( 2.0*M_PI*genrand_real3() );
}
}
return(x);
}

繰り返しを各行に割り当てるため、int repl、各回のベクトル要素数は int length です。
ひとつ前のr_normal関数との違いは、まず毎回の結果を行列に格納している点です。
NumericMatrixが行列型を示しています。
一行目でもxをrepl行、length列の行列として定義します。NumericMatrix x(repl, length);
毎回の繰り返しでは、ひとつ前の結果に0.9を掛けて正規乱数を足していくので、iのforループは1から始まるようにしています（1行目は0のままで放置、Rコードで除去します）。
もう一つ、分かりにくい違いが添字の書き方です。さっきはx[i] = ...と肩のある括弧を使いましたが、行列の場合はx(i,j)と普通の括弧を使用するようです。C++の解説を見るとx[i][j]とありますが、こうするとエラーになりました。x[i,j]もダメです。



それでは冒頭のR版forループと速度比較してみます。
NN <- 10^6
system.time(r_markov(NN+1, 10, 0, 1)[-1,])
 #  ユーザ   システム       経過
 #   0.439      0.024      0.460

R版では5.778秒だったので10倍以上早いです。C++の速さを思い知らされました。

Rcpp、少しずつ部分的にC++化できるので、練習にも持って来いだと思いました。たぶん純粋C++に移行することはなくハイブリッドでやっていくことになりそうな気がします。