Rcpp で疎行列をデータフレームに変換する

しばらく体調を崩してしまっていたため、すっかりご無沙してしまいました。ようやく回復してきましたので活動再開しようと思います。

というわけで、だいぶ間があいてしまいましたが、前回の続きということで、Matrixパッケージで作成した疎行列をデータフレームに変換する処理を紹介します。（正直、必要性はあまりないとは思いますが、Rcpp の使い方を示す例としては悪くないと思います。あと、疎行列の実装方法の一つである圧縮列格納方式については、説明する前に力尽きました。）

ちなみに、R の場合には以下のようにします。

# 疎行列の作成
library(Matrix)
sm <- Matrix(sample(c(0,1), size = 25, replace = TRUE)
             , nrow = 5
             , sparse = TRUE)

# データフレームへの変換
df <- as.data.frame(as.matrix(sm))

これを Rcpp で実装することで、以下の Rcpp 機能の使用例を示します。

S4 クラスのスロットへのアクセス
stop() によるエラー処理
要素数や要素番号の型としての R_len_t, R_xlen_t
diff(), seq(), rep() など R と同じ機能の関数の利用
NULL値の判定
String 型文字列の結合（+=）
サイズを指定したデータフレームの作成
リストからデータフレームへの変換

（追記：2016/06/19）64 bit システムでは R のベクトルの要素数の最大値は理論上「2 ^ 64 - 1」なので、int の最大値「2 ^ 31 - 1」よりも大きくなり得ます。なので、基本的には要素数や要素番号の型としては R_len_t（int型と同じ）ではなく R_xlen_t （ptrdiff_t 型と同じ）を使うべきです。しかし、今回のケースで言うと sm.slot("Dim") は IntegerVector （要素の型は int ）になっているので、変数 nrow, ncol などは R_len_t でも良いです。

実装例

// [[Rcpp::export]]
List sparseMatrixToDataFrame(S4 sm){
    
    // Matrixパッケージの疎行列（dgCMatrix）ではないデータが
    // 渡された場合は拒否します
    if(as<std::string>(sm.attr("class"))!="dgCMatrix")
        stop("Input must be dgCMatrix.");
    
    // 行数と列数
    //要素数や要素番号の型には R_xlen_t を使います
    IntegerVector dim = sm.slot("Dim");
    R_xlen_t nrow = dim[0];
    R_xlen_t ncol = dim[1];
    
    // dgCMatrix から値を取り出します
    // dgCMatrix の疎行列の形式は圧縮列格納方式となっています
    NumericVector X = sm.slot("x");//非ゼロ要素の値
    IntegerVector I = sm.slot("i");//非ゼロ要素の行番号（0ベース）
    IntegerVector P = sm.slot("p");//非ゼロ要素の列番号の開始位置（0ベース）
    
    // P から非ゼロ要素の列番号 J (0ベース) を再生します
    IntegerVector J(X.length());
    IntegerVector repeat = diff(P); //各列番号が連続する数
    R_xlen_t n = repeat.length();
    for(R_len_t j=0; j<n; ++j){
        if(repeat[j]==0) continue; //非ゼロ値が存在しない列は飛ばします
        // J に 各非ゼロ値に対応する列番号 j をセットします
        J[seq(R_xlen_t(P[j]),R_xlen_t(P[j+1]-1))] = rep(j, repeat[j]);
    }
    
    // 出力用データを作成します
    // ここでは列数が ncol のデータフレームを作成したいですが
    // DataFrame df(ncol)；とすると、１列目の値が ncol である
    // データフレームができてしまいます、。
    // そこで要素数が ncol であるリストを作成してから
    // 最後にデータフレームに変換します
    List L(ncol);
    
    // リストを初期化するため
    // 要素の値が 0 で長さが nrow のベクトルを作成して
    // そのベクトルをリストの各要素として代入します
    for(R_xlen_t i=0;i<ncol;++i){
        L[i] = NumericVector(nrow);
    }
    
    // 出力用データに値をセットします
    // 行番号 I、列番号 J、の要素に値 X を代入します
    n = I.length();
    NumericVector column;
    for(R_xlen_t k=0; k<n; ++k){
        // column は L の J[k] 番目の要素（列）への参照となります
        column = L[J[k]];
        // column の I[k] 番目の要素に X[k] の値を代入します
        column[I[k]] = X[k];
    }
    
    // 出力用データにセットする行名・列名を作成します
    List Dimnames = sm.slot("Dimnames");
    CharacterVector rownames(nrow);
    CharacterVector colnames(ncol);
    // 行名
    if(Dimnames[0]==R_NilValue){
        //行名がNULLの場合は、1 2 ... という行名にします
        IntegerVector v = seq_len(nrow);
        rownames = as<CharacterVector>(v);
    } else {
        rownames=Dimnames[0];
    }
    // 列名
    if(Dimnames[1]==R_NilValue){
        //列名がNULLの場合は、V1 V2 ... という列名にします
        //R の paste() 関数を呼んだら遅かったのでやめました
        for(int i=0; i<ncol; ++i){
            colnames[i] = (String("V") += String(i+1));
        }
    } else {
        colnames=Dimnames[1];
    }
    
    // これらの属性に値をセットして返すと
    // R ではデータフレームとして扱われます
    L.attr("row.names") = rownames;
    L.attr("names")     = colnames;
    L.attr("class")     = "data.frame";
    
    return L;
}

実行結果

> sm <- Matrix(sample(c(0,1), size = 25, replace = TRUE)
+              , nrow = 5
+              , sparse = TRUE)
> sm
5 x 5 sparse Matrix of class "dgCMatrix"
              
[1,] 1 1 1 1 .
[2,] . 1 . 1 .
[3,] . . 1 1 1
[4,] 1 1 1 1 .
[5,] . 1 1 1 1
> sparseMatrixToDataFrame(sm)
  V1 V2 V3 V4 V5
1  1  1  1  1  0
2  0  1  0  1  0
3  0  0  1  1  1
4  1  1  1  1  0
5  0  1  1  1  1

パフォーマンス比較

それでは R バージョンと速度を比較してみましょう

> sm <- Matrix(sample(c(0,1), size = 1000000, replace = TRUE)
+              , ncol = 10000
+              , sparse = TRUE)
> microbenchmark::microbenchmark(sparseMatrixToDataFrame(sm)
+                                ,as.data.frame(as.matrix(sm))
+                                ,times = 10
+                                )
Unit: milliseconds
                         expr       min       lq      mean    median        uq      max neval
  sparseMatrixToDataFrame(sm)  38.63905  45.4774  77.11882  55.54285  62.58709 288.8416    10
 as.data.frame(as.matrix(sm)) 129.57834 140.7170 258.36984 261.72998 374.73649 391.0923    10

中央値でみるとRcpp版は 55.5 ミリ秒に対して、R版は 261.7 ミリ秒なので、Rcpp のほうが 1/5 の時間で終わっています。実際のところR版も関数はC言語で実装されているので結構速くて抜くのに苦労したのですが、おかげで有用な知見が得られました。