匹配逻辑语句的 Rcpp 矩阵的子集答案

【问题标题】：Subset of a Rcpp Matrix that matches a logical statement匹配逻辑语句的 Rcpp 矩阵的子集
【发布时间】：2018-04-05 23:08:44
【问题描述】：

在 R 中，如果我们有一个数据矩阵，比如一个 100 x 10 的矩阵 X，以及一个 100 元素的向量 t，它的可能值是 (0, 1, 2, 3)，我们可以很容易地找到 X 的子矩阵 y使用简单的语法：

y = X[t == 1, ]

但是，问题是，如何使用 Rcpp 的 NumericMatrix 做到这一点？
（或者，更一般地说，我怎样才能在 C++ 的任何容器中做到这一点？）

多亏德克的提示，看来

NumericMatrix X(dataX);
IntegerVector T(dataT);
mat Xmat(X.begin(), X.nrow(), X.ncol(), false);
vec tIdx(T.begin(), T.size(), false); 
mat y = X.rows(find(tIdx == 1));

可以做我想做的事，但这似乎太冗长了。

【问题讨论】：

标签： c++ r rcpp

【解决方案1】：

我很乐意将其视为糖。不幸的是，我没有资格实施它。这里仍然是我玩过的许多不同的解决方案。

首先，我必须对Gong-Yi Liao 代码进行一些修改才能使其正常工作（colvec 而不是vec for tIdx 和Xmat.rows(... 而不是X.rows(...：

mat Xmat(X.begin(), X.nrow(), X.ncol(), false);
colvec tIdx(T.begin(), T.size(), false); 
mat y = Xmat.rows(find(tIdx == 1));

其次，这里有三个带有基准的函数，所有子集矩阵都基于逻辑语句。这些函数采用 arma 或 rcpp 参数并返回值两个基于 Gong-Yi Liao 的解决方案，一个是简单的基于循环的解决方案。

n(rows)=100, p(T==1)=0.3

                expr   min     lq median     uq    max
1  submat_arma(X, T) 5.009 5.3955 5.8250 6.2250 28.320
2 submat_arma2(X, T) 4.859 5.2995 5.6895 6.1685 45.122
3  submat_rcpp(X, T) 5.831 6.3690 6.7465 7.3825 20.876
4        X[T == 1, ] 3.411 3.9380 4.1475 4.5345 27.981

n(rows)=10000, p(T==1)=0.3

                expr     min       lq   median       uq      max
1  submat_arma(X, T) 107.070 113.4000 125.5455 141.3700 1468.539
2 submat_arma2(X, T)  76.179  80.4295  88.2890 100.7525 1153.810
3  submat_rcpp(X, T) 244.242 247.3120 276.6385 309.2710 1934.126
4        X[T == 1, ] 229.884 236.1445 263.5240 289.2370 1876.980

submat.cpp

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// arma in; arma out
// [[Rcpp::export]]
mat submat_arma(arma::mat X, arma::colvec T) {
    mat y = X.rows(find(T == 1));
    return y;
}

// rcpp in; arma out
// [[Rcpp::export]]
mat submat_arma2(NumericMatrix X, NumericVector T) {
    mat Xmat(X.begin(), X.nrow(), X.ncol(), false);
    colvec tIdx(T.begin(), T.size(), false); 
    mat y = Xmat.rows(find(tIdx == 1));
    return y;
}

// rcpp in; rcpp out
// [[Rcpp::export]]
NumericMatrix submat_rcpp(NumericMatrix X, LogicalVector condition) { 
    int n=X.nrow(), k=X.ncol();
    NumericMatrix out(sum(condition),k);
    for (int i = 0, j = 0; i < n; i++) {
        if(condition[i]) {
            out(j,_) = X(i,_);
            j = j+1;
        }
    }
    return(out);
}


/*** R
library("microbenchmark")

# simulate data
n=100
p=0.3
T=rbinom(n,1,p)
X=as.matrix(cbind(rnorm(n),rnorm(n)))

# compare output
identical(X[T==1,],submat_arma(X,T))
identical(X[T==1,],submat_arma2(X,T))
identical(X[T==1,],submat_rcpp(X,T))

# benchmark
microbenchmark(X[T==1,],submat_arma(X,T),submat_arma2(X,T),submat_rcpp(X,T),times=500)

# increase n
n=10000
p=0.3
T=rbinom(n,1,p)
X=as.matrix(cbind(rnorm(n),rnorm(n)))
# benchmark
microbenchmark(X[T==1,],submat_arma(X,T),submat_arma2(X,T),submat_rcpp(X,T),times=500)

*/

【讨论】：

【解决方案2】：

我所知道的最接近的是find() 函数与Armadillo 中的submat() 函数的组合，可通过RcppArmadillo 访问。

编辑：这当然是我们可以通过补丁添加的东西。如果有人有足够的动力尝试这个，请来到 rcpp-devel 邮件列表。

【讨论】：

是的，添加它需要大量的开发和测试。因此，除非有专门的资金支持，否则它不太可能很快发生