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

Question

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

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

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

Answer 1

我假设 "tidy way" 你的意思是避免循环以支持使用 Rcpp 中提供的语法糖。由于 sugar 为向量提供了一个比较器，而不是为矩阵提供了一个值（请参阅 here and here），我认为您现在可以做的最多 "tidy way" 是（仅）在列（或行）上循环，即，无需遍历列和行：

// [[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;
}

> f(fmat, 1.0)
      [,1]  [,2]
[1,]  TRUE FALSE
[2,] FALSE  TRUE

> f(fmat, -1.0)
     [,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE

> f(fmat, 2.0)
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE

但是，我建议避免额外的循环并不能真正为您带来任何可读性（实际上可能会损害代码的某些读者的可读性）；考虑循环遍历行和列的函数：

// [[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;
}

我真的不明白这有多难打字，它似乎在本质上是等效的（平均执行时间略低，但中值稍高——参见下面的基准测试），至少对于某些读者，我敢打赌，您在做什么会更容易清楚。

就是说，如果跳过循环对您有帮助，我认为这是您目前能做的最好的事情。

library(microbenchmark)

> microbenchmark(loop = f(fmat, 1.0), nonloop = f2(fmat, 1.0), times = 1e4)
Unit: microseconds
    expr   min    lq     mean median    uq      max neval cld
    loop 6.564 7.402  9.77116  7.612 8.031 9173.952 10000   a
 nonloop 6.425 7.123 10.01659  7.333 7.682 4377.448 10000   a

> microbenchmark(nonloop = f2(fmat, 1.0), loop = f(fmat, 1.0), times = 1e4)
Unit: microseconds
    expr   min    lq      mean median    uq      max neval cld
 nonloop 6.356 7.124 10.179950  7.333 7.544 4822.066 10000   a
    loop 6.775 7.404  9.588326  7.613 7.892 4278.971 10000   a

Answer 2

@duckmayr 的回答确实是 spot-on，并且显示了一个重要的细节：我们不妨将实现细节 隐藏在函数后面 因为毕竟这就是全部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()); 
}

但是 non-obvious element-wise 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() 相关的一行，但对 run-time:

影响不大

// [[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()); 
}

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

How can I do logical operations on Rcpp::NumericMatrix using a "sugar" manner?

r

rcpp