如何使用“糖”方式对 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()); 
}