r - 如何使用 "sugar"方式对 Rcpp::NumericMatrix 进行逻辑运算？

Question

我必须明智地将矩阵条目与数字进行比较，因此我尝试定义一个 Cxx 函数，例如

src <- '
LogicalMatrix f(NumericMatrix fmat, double t){
  LogicalMatrix result = fmat >= t;
  return result;
}
'
cppFunction(src)

但是会抛出一些异常。是什么原因？那么我怎样才能整齐地做到这一点呢？

Answer 1

@duckmayr 的回答非常准确，并且显示了一个重要的细节:我们最好将实现细节隐藏在一个函数后面，因为毕竟这就是Rcpp Sugar 等人无论如何都为我们做。

但是如果我们首先将矩阵转换为向量，对该向量进行操作，然后恢复矩阵，我们就可以按照@zengchao 的要求依赖 Sugar 操作。这是可行的，因为在内部矩阵只是一个增加了维度的向量(二阶；数组一般化为多于两个)。

但事实证明......该版本比循环(略)昂贵(并且比在列上工作略便宜)。请参阅下面的完整详细信息，但函数 f3() 可能是:

// [[Rcpp::export]]
LogicalMatrix f3(NumericMatrix fmat, double t) {
    IntegerVector dims = fmat.attr("dim");
    NumericVector v(fmat);
    LogicalVector lv = v >= t;
    return LogicalMatrix(dims[0], dims[1], lv.begin()); 
}

但非显而易见的元素 f2() 仍然是最快的:

R> microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 5e4)
Unit: nanoseconds
       expr min  lq    mean median     uq     max neval
  f(mat, 1) 873 992 1322.10   1042 1118.0 1493236 50000
 f2(mat, 1) 823 925 1195.49    975 1049.5 2068214 50000
 f3(mat, 1) 864 977 1288.68   1031 1114.0 1909361 50000
R> 

道德:简单的循环解决方案对临时对象的复制最少，而且速度最快。总的来说，三者之间的速度差异无关紧要。

对于更大的矩阵，不复制临时对象的优势变得更加重要:

R> mat <- matrix(sqrt(1:1000), 1000)

R> microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 1e3)
Unit: microseconds
       expr   min    lq    mean median     uq    max neval
  f(mat, 1) 3.720 3.895 3.99972 3.9555 4.0425 16.758  1000
 f2(mat, 1) 1.999 2.122 2.23261 2.1760 2.2545 17.325  1000
 f3(mat, 1) 3.921 4.156 4.31034 4.2220 4.3270 19.982  1000
R> 

完整代码如下。

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
LogicalMatrix f(NumericMatrix fmat, double t){
    int n = fmat.nrow(), m = fmat.ncol();
    LogicalMatrix result(n, m);
    for ( int j = 0; j < m; ++j ) {
        result(_, j) = fmat(_, j) >= t;
    }
    return result;
}

// [[Rcpp::export]]
LogicalMatrix f2(NumericMatrix fmat, double t){
    int n = fmat.nrow(), m = fmat.ncol();
    LogicalMatrix result(n, m);
    for ( int i = 0; i < n; ++i ) {
        for ( int j = 0; j < m; ++j ) {
            result(i, j) = fmat(i, j) >= t;
        }
    }
    return result;
}

// [[Rcpp::export]]
LogicalMatrix f3(NumericMatrix fmat, double t) {
    int dims[2] = { fmat.nrow(), fmat.ncol() };
    NumericVector v(fmat);
    LogicalVector lv = v >= t;
    return LogicalMatrix(dims[0], dims[1], lv.begin()); 
}

/*** R
mat <- matrix(c(1,2,3,4), 2, 2)
library(microbenchmark)
microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 1e5)

mat <- matrix(sqrt(1:1000), 1000)
microbenchmark(f(mat, 1.0), f2(mat, 1.0), f3(mat, 1.0), times = 1e3)
*/

编辑:我们可以再删除与 f3() 相关的一行，但这对运行时影响不大:

// [[Rcpp::export]]
LogicalMatrix f4(NumericMatrix fmat, double t) {
    IntegerVector dims = fmat.attr("dim");
    LogicalVector lv = NumericVector(fmat) >= t;
    return LogicalMatrix(dims[0], dims[1], lv.begin()); 
}