【发布时间】：2016-04-30 19:45:52
【问题描述】：
我有一个关于从矩阵到向量的子集的问题。用户可以明确地给出索引矩阵（它是一个与 M 大小相同的矩阵，如果不需要条目，则为 0，如果必须提取条目，则为 1）。如果提供了 indexmatrix，那么我们只是将其子集化，如果没有提供 indexmatrix（indexmatrix = NULL），那么我们使用 type1（取 true 或 false）构建它。只有两种类型的索引矩阵是可能的。
我使用了 Subset of a Rcpp Matrix that matches a logical statement
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
// [[Rcpp::export]]
arma::colvec extractElementsRcpp(arma::mat M,
Rcpp::Nullable<Rcpp::NumericMatrix> indexmatrix = R_NilValue,
bool type1 = false) {
unsigned int D = M.n_rows; // dimension of the data
arma::mat indmatrix(D, D); // initialize indexmatrix
if (indexmatrix.isNotNull()) {
// copy indexmatrix to numericmatrix
Rcpp::NumericMatrix indexmatrixt(indexmatrix);
// make indexmatrix into arma matrix indmatrix
indmatrix = Rcpp::as<arma::mat>(indexmatrixt);
} //else {
// get indexmatrix
// Rcpp::NumericMatrix indexmatrixt = getindexmatrix(D, type1)["indexmatrix"];
// // make indexmatrix into arma matrix
// indmatrix = Rcpp::as<arma::mat>(indexmatrixt);
// }
arma::colvec unM = M.elem(find(indmatrix == 1)); // extract wanted elements
return(unM);
}
效果很好，太好了！然而，速度并不是我所希望的。每当提供 indexmatrix 时，C++ 代码都会比普通的 R 代码慢，而我的目标是在速度上有一个很好的改进。我觉得我在复制矩阵太多了。但我是 C++ 新手，还没有找到避免它的方法。
速度对比如下：
test replications elapsed relative
2 extractElementsR(M, indexmatrix = ind) 100 0.084 1.00
1 extractElementsRcpp(M, indexmatrix = ind) 100 0.142 1.69
编辑：R函数定义为
extractElementsR <- function (M, indexmatrix, type1 = FALSE) {
D <- nrow(M)
# # get indexmatrix, if necessary
# if(is.null(indexmatrix)) indexmatrix <- getindexmatrix(D, type1 = type1)$indexmatrix
# extract wanted elements
return (M[which(indexmatrix > 0)])
}
例如可以采用矩阵
M <- matrix(rnorm(1000^2), ncol = 1000)
indexmatrix <- matrix(1, 1000, 1000)
indexmatrix[lower.tri(indexmatrix)] <- 0
作为 M 和索引矩阵。
EDIT2：我在 Rcpp 函数中注释了 else 语句并省略了 R 函数中的默认 NULL 值，因为它对我的问题并不重要。我想在提供 indexmatrix 时提高 Rcpp 函数的速度。但是，我想保留默认的NULL 值（并在必要时创建和索引矩阵）。
【问题讨论】：
-
什么是
getindexmatrix？ -
它是一个构造索引矩阵的函数，根据维度和类型，如果没有索引矩阵传递给函数extractElementsR(cpp)。它的细节在这里并不重要，因为我想专注于在传递索引矩阵时提高速度，同时我希望函数保持索引矩阵的默认值为 NULL。
-
当 type1 = T; 时，您可以将函数 getindexmatrix 视为构造索引矩阵的函数，如示例矩阵中提供的那样。当 type1 = F 时为 1-indexmatrix。当然与 M 具有相同的维度。
-
实际上它很重要，因为没有它你的代码就无法运行。此外，这可能是导致
extractElementsRcpp性能缓慢的另一个因素。请包含此功能，因为您的示例目前无法重现。 -
我明确提到我在提供 indexmatrix 时进行了基准测试，因此这部分代码永远不会运行，它对性能没有贡献。您可以在
extractElementsRcpp中注释else部分并在R 代码中省略NULL默认值而不改变意图。
标签： r performance subset default-value rcpp