須通り_統計_C++は未経験だけどRcppで開発したい！(2) 関数自作編

ホーム | 統計 Top | C++は未経験だけどRcppで開発したい！(2) 関数自作編

（2015年5月20日初版）

実践 1：Rcppで簡単な関数を定義してみよう

既存コードのコンパイルが通ることは前回記事で既に確かめているので、満を持してRcppでの開発を始めよう。簡単な計算の例として、「ベクトルを読み込んで要素の最大値を返す」処理を実装してみる。まあRでは既存の max() 関数一発で済む処理だが、練習としては手頃なサイズでそれなりに応用も利くので、是非とも試していただきたい。

ベクトルを読み込んで要素の最大値を返す関数

適当なテキストエディタで新規ファイルを作成。以下の内容を書き込んだら、普段Rのホームディレクトリにしているフォルダに「名前をつけて保存」する。拡張子は.cppとし、とりあえずファイル名を "mysource01.cpp" とでもしておこう。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double returnMaxCpp(const NumericVector& x) {
    const int length = x.length(); // length of given vector x
    double temp_max = x[0];
    for(int i=0; i < length; ++i) {
        if (x[i] > temp_max) {
            temp_max = x[i];
        }
    }

    return(temp_max);
}

なおsourceCppでの評価は以下。コンパイルとR環境へのリンクを自前で行っているので、数秒を要する。

library(Rcpp) # RcppパッッケージをRにロード
sourceCpp("mysource01.cpp") # ローカルのC++ソースファイルを評価

# 使ってみる
> returnMaxCpp(c(-1,0,1,2,3))
[1] 3

Rcppでの関数定義の約束事

#include <Rcpp.h>
using namespace Rcpp;

冒頭2行はおまじないだが、一応解説しておく。最初の行はRcppの標準ヘッダを読み込む命令で、Rcpp.hというファイルの中身はこんな感じ。2行目は「using指令」というもので、このソースファイルでこれから定義する諸関数を"Rcpp"という名前空間に所属させてくださいね、という指令を下している。要するに、もしもR本体で太郎や花子といった名前の関数が既に定義されているとき、同名の関数ができて識別不能になる危険を回避するために、Rcpp家の太郎とか花子とかいった形をとるのだ。

// [[Rcpp::export]]
double returnMaxCpp(const NumericVector& x) {

その次にようやくRcppでの関数（ここでは returnMaxCpp() のこと）の定義が来るが、前にある1行
// [[Rcpp::export]]
は、returnMaxCpp()という関数を素のR環境上からも呼び出し可能にしてくださいね、という指示である。ソースファイル内の、別のRcpp関数からのみ呼び出すサブルーチンなどを書く際には、あえて付けないこともある。 // という2つのスラッシュは、本来C++のソースコードで1行コメントの開始を表す記号なのだが、この場合には意味を持つ1行なので直後の関数定義の行との間に別のコメントを挟んだりしてはいけない（筆者はこれに引っかかりました、ハイ）。

関数定義の中身について、これからRの関数定義と見比べながら説明してゆく。一応アルゴリズムの概要を述べておくと、ベクトルの最初の要素を変数 temp_max （暫定王者の座）に格納し、以降の各要素とタイマンで比較してゆく。より大きな要素があればそいつを新たな暫定王者に据えてゆき、最後に残った値がベクトル全体の最大要素だ。

注意：この関数は必ずしも実用に耐える訳ではない。幾つか問題を挙げると、1) 引数であるベクトルxにNAとなる要素があったら？2) 特にxの最初の要素がNAだったら、比較の結果も常にNAしか返せないのでは？3) 文字列などがあったら予期しない結果になるのでは？といった具合である。手元でコンパイルしたところ、NAを含むベクトルは正しく最大値を返したが、文字列を含むとき「エラー: not compatible with requested type」となった。

RcppとRのソースコードを比べてみよう

違いを挙げ出すときりがないので、先ほどからの例で特に皆さんが気になるであろう箇所を簡潔に。まずRcppでの関数定義を再掲し、続いて同じ処理をR環境で行う関数の例を出す。アルゴリズムはわざと同一にしたので、対応する箇所を見比べるとよい。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double returnMaxCpp(const NumericVector& x) {
    const int length = x.length(); // length of given vector x
    double temp_max = x[0];
    for(int i=0; i < length; ++i) {
        if (x[i] > temp_max) {
            temp_max = x[i];
        }
    }

    return(temp_max);
}

# 上と同じ働きをする関数を pure R で実装した例
returnMaxR <- function(x) {
    length <- length(x) # length of given vector x
    temp_max <- x[1]
    for(i in 1:length) {
        if (x[i] > temp_max) {
            temp_max = x[i]
        }
    }

    return(temp_max)
}

まず double returnMaxCpp(const NumericVector& x) という一行からして、Rと様相が異なる。最大の違いは、Rの変数定義においてユーザが「型」を明示的に与えないことだ。逆に言えば、Rcpp（C++）ではユーザが変数型、ひいてはメモリ領域の割当・管理に責任を負っている。

データ型について

上の例では関数定義そのものが double 型の変数の宣言とも読め、これは関数の戻り値が、実数型の値であることを示している。また引数も const NumericVector& x となっているが、これは x が数値ベクトル（Rcppでは NumericVector というデータ型で扱う）で与えられることを示している。先に少し述べたように const は、returnMaxCpp関数のスコープ内で、x という変数の値を上書きすることを禁じている（この場合はいわゆる「定数」というより「データ保護」の意味合いで使われている）。

なお NumericVector にアンド & 記号が付いているが、これは引数をC++でいうところの「参照渡し」にするものだ。早い話が既存の変数を下敷きにして新たな変数を作る際に、メモリ上の別のアドレスへ値をコピーする（値渡し：Rの変数定義は多くの場合これ）か、コピーはせずに既存のメモリ領域を参照させるかの違いである。参照渡しの方が当然メモリを節約できるし処理も概ね高速だが、自分では新しい変数の値を書き換えたと思っていても、実は参照元の変数の中身を変更してしまっているという副作用もある。意図してそうしたい場合を除き、（小さなデータならば）値渡しにするか、constを付けて上書きを阻止するのが安全だ。

型をユーザーに意識させることの是非は一概には言えない。小規模なデータ処理であれば R の持ついい加減さが、「使い勝手のよさ」や「学びやすさ」となる。反面、大規模なデータを少ない物理メモリ上でやり繰りするとか、反復計算をひたすら高速に繰り返すといった仕事には、C++の厳密なメモリ管理はありがたい。

オブジェクト指向について

これも、RとC++の各言語仕様に起因する大きな差異である。オブジェクト指向そのものの解説は本ページのスコープを超えるし、筆者も十分には理解できていないので今回はパス。こいつが本格的に顔を出すのが、以下の行だ。

    const int length = x.length(); // length of given vector x
    double temp_max = x[0];

冒頭の const int length = x.length(); とは、整数型の変数（ただしconst と頭についているので値の変更は許されず、定数扱い）である length を定義し、それに格納する値を x というベクトルから、x自身の長さ（要素の数）として取得する処理だ。厳密には「Rcppクラスの派生クラスであるRcpp::NumericVectorに属するオブジェクトxに対するメソッドlength」なのだが、細けぇことは気にしないで「NumericVectorという形式のデータに対して、その長さを知る手続きが用意されている」という事実を知っておけばいい。

このようにRcppにおいて . ドットには「ドット演算子」という特別な意味があるので、Rでしばしばやるように変数名の一部（単語の繋ぎ）として使ってはいけない。じゃあどうするべきか。一つの方法は _ アンダーバーで代用することであり、もう一つはキャメルケースの使用である。ハイフン - は、Rにおいてマイナス演算子と区別できないので推奨しない。

次の temp_max についても、変数初期化が、NumericVector のi番目の要素を取り出すという手続きで与えられている。ただし見た目は、R環境でベクトルの成分を x[何番目の要素] として取り出す処理そのまんまだね。一応、C++とRではベクトルの添字を0始まりにするか、1始まりにするかという違いがあるので注意しておこう。

C++のループの書式

以下のループについて、for文の中身 {ここ} 自体はRと殆ど変わらない。せいぜい各命令文の後に、セミコロンが付く程度。

    for(int i=0; i < length; ++i) {
        if (x[i] > temp_max) {
            temp_max = x[i];
        }
    }

問題は for(int i=0; i < length; ++i) という、ループ制御部分の指定方法だ。定型文なので暗記してもらうとして、ループを1～nではなく、0～(n-1)とする必要があるのは要注意。先に述べたようにC++では、長さnのベクトルの最初の要素のアドレスは "0" であり、最後の要素は "n-1" である。よって順にアクセスしたい場合、最初は i=0 として、最後は i < n で終わらせる。

その後ろの ++ はインクリメント演算子。ループ内の計算が1周期分完了するタイミングで、オペランドである i の値を1つずつ増やす。

厳密には ++i は前置インクリメントと呼ばれる。後置インクリメントというのもあって i++ と書ける。両者の間には、i++はインクリメントする直前に i というインスタンスのコピーを作り、++iはコピーしないという違いがあるため、しばしば++iのほうが速いと言われている。どちらを使っても計算結果は変わらない。

計算速度は？

本記事をお読みの諸兄姉方が最も気にする点だろう。（この欲張りさんめ☆）だが一寸した計算では R, Rcpp とも瞬時に完了してしまうので、比較にならない。そこで長さ1億のベクトルを作って、最大値を求めてみる。

> system.time(returnMaxCpp(c(1:10000000)))
   ユーザ   システム       経過
     0.151      0.048      0.198
> system.time(returnMaxR(c(1:10000000)))
   ユーザ   システム       経過
    12.343      0.082     12.325

(参考)
> system.time(max(c(1:10000000)))
   ユーザ   システム       経過
     0.054      0.023      0.077

素のR環境でループをべた書きした関数に比べ、同じアルゴリズムをRcppで実装した関数は2オーダーほど高速である。この「約100倍弱」が、RとRcppの標準的なパフォーマンスの差だと考えて頂ければよい。とはいえ、Rに組み込まれているmax()関数の方がさらに2倍以上速かった。ざんねん！そもそもRの基本関数には、内部的にCで実装されていて既に十分高速なものも多いし、最大値を求めるアルゴリズム自体の巧拙にも起因するのだろう。

実践 2：データ型の違いを意識した計算を行おう

さて前節では Rcpp で関数を書いてみたわけだが、intとかdoubleとかNumericVectorとか様々なデータ型を使い分けており、これまでもっぱらR畑で仕事をしてきた人は戸惑ったかもしれない。とはいえ、最低限以下のデータ型を覚えて何となく使い分ければ、それなりの計算用途には事足りるようだ。

int 整数
double 実数
string 文字列
NumericVector 実数ベクトル
NumericMatrix 実数行列
List リスト

特にRとRcppの違いで注意せねばならないのは、intとdoubleの使い分けである。Rは基本的に数値をnumeric（実数）で扱い、これはC++でいうところの double（倍精度の実数）に他ならない。Rcpp（C++）の整数、たとえば「1」と、実数「1.0」とは明確に別物である。試しに以下をRコンソールで実行してみよ。

library( Rcpp )

test_division <- '
    #include <Rcpp.h>
    using namespace Rcpp;

    // [[Rcpp::export]]
    int divisionAsInt () {
        int FiveAsInteger = 5;
        int ThreeAsInteger = 3;
        int ans = FiveAsInteger/ThreeAsInteger;
        return ans;
    }

    // [[Rcpp::export]]
    double divisionAsDouble () {
        double FiveAsDouble = 5.0;
        double ThreeAsDouble = 3.0;
        double ans = FiveAsDouble/ThreeAsDouble;
        return ans;
    }
'
sourceCpp(code=test_division)

divisionAsInt() # 計算結果は 1 になる
divisionAsDouble() # 計算結果は1.666667 になる

戻り値が複数必要なら属性を使う

これはR環境でも有効な手なのだが、関数から複数の値を返させたい、しかもデータ型の異なる値が混在するのでベクトルにはできないという場合、attributeを使う。こいつはRcppでもほぼ同じ方法で使える。サンプルコードは現在準備中なので暫しお待ちを。

実践 3：外部ライブラリの力を借りて複雑な処理に挑戦しよう

作成中。

実践 4：Rcppの関数をパッケージ化しよう

Rcppコードはパッケージ化せずにsourceCppから使うことも出来るし、開発段階ではその方が手間を省きやすいが、安定動作するようになったサブルーチンは順次パッケージ化するのが何かと安全であろう。サンプルコードは現在準備中なので暫しお待ちを。

補足「なぜRcppなのか」

世のRユーザーの多くは、変数のデータ型とかメモリ管理とか、コンピュータの原始的な部分を意識せずともデータを解析している。それでいて解析の一貫性や再現性を確保することも、コマンド入力やバッチ処理のおかげで、R環境においては容易である。この特徴こそがライトユーザーからコアユーザーまで幅広い層を惹き付け、今日のRを統計解析環境のデファクトスタンダードにまで押し上げた。正にインタプリタ言語の長所が、最大限生かされた例であろう。もちろん柔軟性の代償として、計算速度はコンパイラ言語（C++など）には概ね負ける。ゆえに「軽い作業はRで、重い作業はRcppへ移植して」という用途に応じた使い分けは、シミュレーションモデルの開発において合理的な戦略である。

純粋なC++でライブラリを書く作業に比べ、Rcppの真の強みは、「C++がシームレスにR環境へ統合される」点にある。これこそが、開発フローに一つの必然をもたらす。大規模かつ高速計算が求められるモデルを、最終的にはRcppで組み上げるつもりであっても、構造が確定するまではR言語でこねくり回すのが楽である。言うなればRは、Rcppのプロトタイピングツールとして活用しうるのだ。幸いなことに、現在のpackage Rcppの公式は、この流れを少なからず意識していると思われる。まずR環境でモデル構造を考え、サブルーチン単位で少しずつRcppに書き換えてはsourceCppでインタラクティブに評価し、最終的にはRcppで自作パッケージをリリースするという黄金パターンができつつある。

この開発フローはそれ自体が効率的であると同時に、筆者のようにコンパイラ言語に不慣れである一般研究者に対しても、プログラミング行為の敷居を大きく下げるものだ。まあ一般的にどこまで有効な手法であるかは、これから事例を蓄積せねばなるまい。読者の諸兄姉方におかれては、ぜひともR+Rcppでのモデル開発を試していただき、その所感をドキュメントとして公開してもらえると幸いである。

目次