非隣接サブマトリクス選択のRcppArmadilloとRランニングスピードの比較

Question

私の主な目標は、行と列に2組のバイナリベクトルを使用して非連続サブマトリクスを選択することです。これは、Rcpp、RcppArmadillo、RcppEigenを使ってC++で実装しているMCMCループのために必要な多くのステップの1つです。

これを行うには、（1）RcppArmadilloを使用する、（2）RcppからR関数を呼び出し、（3）Rを直接使用して結果をC++に渡す方法が考えられます。最後の選択肢は私にはまったく便利ではありませんが。

次に、これら3つのシナリオのパフォーマンススピードを比較しました。興味深いことに、ダイレクトRコードは他の2つのコードよりはるかに高速です！ Rcppから正確なR関数を呼び出すと、Rから直接呼び出すよりもはるかに遅いです。私は、これらの関数が、このolder postの例で示唆されているように比較的同じ実行速度を持つことを期待していました。

とにかく、タイミングの結果はちょっと奇妙なようです。その理由に関するコメント？私は、エルキャピタンのOS、2.5GHz Intel Core i7でMacBook Proを使用しています。それは私のシステム、Mac OSX、またはRcppが私のマシンにインストールされている方法に関係していますか？

ありがとうございます！ここで

は、コードは次のとおりです。

CPP章：

#include <RcppArmadillo.h> 

// [[Rcpp::depends(RcppArmadillo)]] 
using namespace Rcpp; 
using namespace arma; 

// (1) Using RcppArmadillo functions: 
// [[Rcpp::export]] 
mat subselect(NumericMatrix X, uvec rows, uvec cols){ 

    mat XX(X.begin(), X.nrow(),X.ncol(), false); 
    mat y = XX.submat(find(rows>0),find(cols>0)); 
    return (y); 
} 

// (2) Calling the function from R: 
// [[Rcpp::export]] 
NumericalMatrix subselect2(NumericMatrix X, NumericVector rows, NumericVector cols){ 

    Environment stats; 
    Function submat = stats["submat"]; 
    NumericMatrix outmat=submat(X,rows,cols); 
    return(wrap(outmat)); 
} 

Rセクション：

library(microbenchmark) 

# (3) My R function: 
submat <- function(mat,rvec,cvec){ 
return(mat[as.logical(rvec),as.logical(cvec)]) 
} 

# Comparing the performances: 

// Generating data: 
set.seed(432) 
rows <- rbinom(1000,1,0.1) 
cols <- rbinom(1000,1,0.1) 
amat <- matrix(1:1e06,1000,1000) 

//benchmarking: 
microbenchmark(subselect(amat,rows,cols), 
      subselect2(amat,rows,cols), 
      submat(amat,rows,cols)) 

結果：

      expr  min  lq  mean median  uq  max neval 

    subselect(amat, rows, cols) 893.670 1566.882 2297.991 1675.282 2184.783 8462.142  100 
subselect2(amat, rows, cols) 928.418 1581.553 3554.805 1657.454 2060.837 138801.050  100 
    submat(amat, rows, cols) 36.313 55.748 66.782 62.709 73.975 136.970  100 

Answer 1

ここで扱う価値のあるものがいくつかあります。まず、アルマジロ機能のパフォーマンスに重大な影響を与えたベンチマークの設計には微妙な間違いがありました。subselect観察：

set.seed(432) 
rows <- rbinom(1000, 1, 0.1) 
cols <- rbinom(1000, 1, 0.1) 

imat <- matrix(1:1e6, 1000, 1000) 
nmat <- imat + 0.0 

storage.mode(imat) 
# [1] "integer" 

storage.mode(nmat) 
# [1] "double" 

microbenchmark(
    "imat" = subselect(imat, rows, cols), 
    "nmat" = subselect(nmat, rows, cols) 
) 
# Unit: microseconds 
# expr  min  lq  mean median  uq  max neval 
# imat 3088.140 3218.013 4355.2956 3404.4685 4585.1095 21662.540 100 
# nmat 139.298 167.116 223.2271 209.4585 238.6875 533.035 100 

---

Rは、しばしば、浮動小数点値として整数リテラル（例えば、1、2、3、...）扱いがが、配列オペレータ:はこれにいくつかの例外の一つである、

storage.mode(c(1, 2, 3, 4, 5)) 
# [1] "double" 

storage.mode(1:5) 
# [1] "integer" 

です。そのため、式matrix(1:1e6, 1000, 1000)は、numericではなく、integerの行列を返します。 subselectはIntegerMatrixではなくNumericMatrixを予期しており、後者のタイプを渡すとディープコピーがトリガーされるため、上記のベンチマークでは桁違いの差が生じるため、問題があります。

---

第二に、Rの関数submatとにより、基礎となるアルゴリズムの違いと考えられるバイナリインデキシングベクトルの分布を超えるC++関数subselect、の相対的な性能との間の顕著な違いがあります。より疎な索引付け（1より大きい0の割合）の場合、R関数が優先されます。より高密度の索引付けの場合、その逆が真です。これは、以下のプロットに示すように、行列のサイズ（または場合によっては次元）の関数であるようにも見えます.012,06,015,、成功パラメータ0.0,0.05,0.10、...、0.95,1.0を使用して行と列のインデックスベクトルが生成されますまず、1e3 x 1e3行列、次に1e3 x 1e4行列を使用します。このコードは最後に含まれています。

---

---

ベンチマークコード：

library(data.table) 
library(microbenchmark) 
library(ggplot2) 

test_data <- function(nr, nc, p, seed = 123) { 
    set.seed(seed) 
    list(
     x = matrix(rnorm(nr * nc), nr, nc), 
     rv = rbinom(nr, 1, p), 
     cv = rbinom(nc, 1, p) 
    ) 
} 

tests <- lapply(seq(0, 1, 0.05), function(p) { 
    lst <- test_data(1e3, 1e3, p) 
    list(
     p = p, 
     benchmark = microbenchmark::microbenchmark(
      R = submat(lst[[1]], lst[[2]], lst[[3]]), 
      Arma = subselect(lst[[1]], lst[[2]], lst[[3]]) 
     ) 
    ) 
}) 

gt <- rbindlist(
    Map(function(g) { 
     data.table(g[[2]])[ 
      ,.(Median.us = median(time/1000)), 
      by = .(Expr = expr) 
     ][order(Median.us)][ 
      ,Relative := Median.us/min(Median.us) 
     ][,pSuccess := sprintf("%3.2f", g[[1]])] 
    }, tests) 
) 

ggplot(gt) + 
    geom_point(
     aes(
      x = pSuccess, 
      y = Relative, 
      color = Expr 
     ), 
     size = 2, 
     alpha = 0.75 
    ) + 
    theme_bw() + 
    ggtitle("1e3 x 1e3 Matrix") 

## change `test_data(1e3, 1e3, p)` to 
## `test_data(1e3, 1e4, p)` inside of 
## `tests <- lapply(...) ...` to generate 
## the second plot