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

前回までは、Rcpp を使ってファイルの行数をカウントするコードを紹介しました。しかし、紹介したコードでは C/C++ の関数やクラスしか使わなかったため、Rcpp の機能が活用されていない感じになってしまいました。そこで今回は Rcpp の機能を使って、いかに R との連携がスムーズに行われるか見て見たいと思います。その例として、データフレームを疎行列に変換するプログラムを紹介します。

今回紹介するプログラムは Rcpp がもつ R との連携機能を色々活用しています。

DataFrame の要素へのアクセス
DataFrame の行数や列数
NA 値の判定
パッケージの関数の呼び出し
オブジェクトの属性値（列名、行名）へのアクセス
NULL値のセット
標準 C++ のデータ構造を利用

コード例

まずは、さっそくコード例を見てみましょう。ちなみに、ここで紹介しているプログラムでは、カテゴリカル変数のダミー変数への変換は、既に行われているという前提になっておりますので、あしからず。

// [[Rcpp::export]]
S4 asSparseMatrix( DataFrame df, bool na = 0){
  // df : データフレーム、値は数値のみである前提
  // na  df にNAが含まれる場合は 1、 含まれていない場合は 0 にセット
  
  //データの行数と列数  
  int nrow = df.nrows();
  int ncol = df.length();
  
  // R の Vector では要素の追加（push_back）の効率が悪いので
  // 標準 C++ の vector を使用する。
  // R_xlen_t は行番号や列番号の型
  std::vector<R_xlen_t>    rows;
  std::vector<R_xlen_t>    cols;
  std::vector<double>      vals;
  

  //データフレームの全ての要素にアクセスして、
  //0じゃない要素の位置と値を保存する
  if(na){//NAありの場合
    //データがNAを含む可能性がある場合は
    //NAのチェックを行う
    for(R_xlen_t col=0; col<ncol; ++col){
      NumericVector column = df[col];
      for(R_xlen_t row=0; row<nrow; ++row){
        // NA の要素は 0 とみなす
        if(NumericVector::is_na(column[row])) continue;
        else if(column[row]!=0.0){
          rows.push_back(row+1);
          cols.push_back(col+1);
          vals.push_back(column[row]);
        }
      }
    }
  } else {
    //データがNAを含まないと保証された場合は
    //NAのチェックを行わない
    for(R_xlen_t col=0; col<ncol; ++col){
      NumericVector column = df[col];
      for(R_xlen_t row=0; row<nrow; ++row){
        if(column[row]!=0.0){
          rows.push_back(row+1);
          cols.push_back(col+1);
          vals.push_back(column[row]);
        }
      }
    }
  }
  
  // Matrix パッケージから sparseMatrix 関数を読み込み
  // Matrix::sparseMatrix() という形式で呼び出したのと同じ
  Environment env = Environment::namespace_env("Matrix");
  Function sparseMatrix = env["sparseMatrix"];
  // wrap() を使って std::vector を NumericVector に変換している。
  S4 sm = sparseMatrix( Named("i") = wrap(rows),
                        Named("j") = wrap(cols),
                        Named("x") = wrap(vals),
                        Named("dims") = NumericVector::create(nrow,ncol));
  //疎行列に列名をセットする、行名は NULL にする
  List dimnames = List::create(R_NilValue, df.names());
  sm.attr("Dimnames") = dimnames;
  return sm;
}

実行例

x <- sample(c(0,1), size = 5, replace = TRUE)
y <- sample(c(0,1), size = 5, replace = TRUE)

data <- data.frame(x,y)

sm1 <- Matrix::Matrix(as.matrix(data), sparse = TRUE)
sm2 <- asSparseMatrix(data)

identical(sm1, sm2) // TRUE

この関数を作成した背景

Matrix パッケージの疎行列クラス（dgCMatrix）を分析で使用しているのですが、データフレームから疎行列を作成するためには、通常はMatrix(as.matrix(data), sparse = TRUE) という感じで疎行列に変換します。しかし、データが巨大になってくると、行列のサイズの限界により as.matrix() での行列変換ができなくなってしまうのです。それならデータフレームから疎行列を直接変換すればいいじゃん、というのが自然だと思うのですが、残念ながら Matrix パッケージでは対応してくれておりません。

実は、Matrixパッケージには、疎行列を作成するもう一つの関数 sparseMatrix() があり、こちらは全ての非ゼロの要素の行番号・列番号・値を入力することで疎行列を作成することができます。そのためにはデータの全ての要素にアクセスして値が、ゼロかどうか判定し、行番号・列番号と値を記録しなければなりませんが、このように値をいちいちチェックする処理は R の苦手とする処理の１つです。そこで Rcpp の出番というわけです。

標準 C++ データ構造を使用する理由

コード例では、非ゼロの要素の行番号・列番号・値を記録するために、NumericVector の代わりに、標準 C++ のデータ構造 std::vector を用いています。データフレームの中に非ゼロ要素がいくつあるかは、全ての要素をチェックしてからでないとわからないので、あらかじめ決められたサイズのベクトルに値を保持することが難しいです。そこで、ベクトルのサイズを動的に増加させることができると嬉しいのですが、Rcpp の Vector（というか R のベクトル）は、一旦ベクトルを作成した後で、要素のサイズを変更する（v <- c(v, 1)）処理を行うと、ベクトルの要素の値を全てをコピーする処理が発生するので、非常に効率が悪いという問題があります。しかし、標準 C++ の std::vector は、新たな要素の追加の効率に優れていますので、こちらを使用しました。

パフォーマンス

Rcpp のコード例のプログラムのパフォーマンスを評価するため、比較対象として、Matrix(as.matrix(data), sparse = TRUE) の他に、コード例と同じ処理を R で実装したバージョンを用意しました。今時なら dplyr パッケージを利用すればもっと効率がよい実装が可能ですが、Rの基本機能で実装すると以下のようにになると思います。

asSparseMatrix_R <- function(data){
  tmp <- 
    lapply(1:ncol(data), function(col){
      row <- which(data[,col]!=0)
      data.frame(row=row, col=col, val=data[row,col])
    }) 
  res <- do.call("rbind", tmp)
  Matrix::sparseMatrix(i = res[,1],
                       j = res[,2],
                       x = res[,3],
                       dims = c(nrow(data),ncol(data)))
}

それではパフォーマンスを比較してみましょう。

> nrow <- 10000
> ncol <- 100
> m <- matrix(sample(c(0,1), size = nrow*ncol, replace = TRUE)
+             , nrow = nrow)
> data <- as.data.frame(m)
> 
> microbenchmark::microbenchmark(
+   Matrix::Matrix(as.matrix(data), sparse = TRUE),
+   asSparseMatrix(data, na = 0),
+   asSparseMatrix_R(data),
+   times = 10)

Unit: milliseconds
                                           expr       min        lq       mean     median         uq       max neval
 Matrix::Matrix(as.matrix(data), sparse = TRUE)  35.88858   41.0130   97.69263   49.57377   58.56511  306.7089    10
                           asSparseMatrix(data) 113.41807  127.1215  159.82937  134.94027  150.90175  391.1737    10
                         asSparseMatrix_R(data) 863.00960 1104.5381 1248.25324 1184.21548 1485.64696 1579.7606    10

残念ながら、この比較では Matrix() 関数に負けてしまいましたが、まずまずの結果です。もともとMatrix() 関数が使えない状況のためにasSparseMatrix() は作成したので、これで十分でしょう。ここで使用したデータはサイズが7MBくらいのデータですが、数十GB程度のデータに対して、それほど時間もかからず処理できるので、業務でも役に立ってくれています。

というわけで、今回は以上です。いかがでしたでしょうか、ちょっと今回は機能が盛り沢山な例になってしまったので、消化するのが大変かもしれないですね。Rcppの様々な機能の使い方が端的に示されていますので、ゆるゆると学んでいただけば良いかなと思います。

次回は、今回とは逆に疎行列をデータフレームに変換する処理を紹介する予定です。