Rcpp 中许多对象的交错结果

Question

我需要将出现在列表中的矩阵和稀疏矩阵逐行写入文件，我正在做这样的事情：

#include <RcppArmadillo.h>

// [[Rcpp::export]]
bool write_rows (Rcpp::List data, Rcpp::CharacterVector clss, int n) {

  int len = data.length();
  for(int i = 0; i<n; i++) {
    for(int j=0; j<len; j++) {
      if (clss[j] == "matrix") {
        Rcpp::NumericMatrix x = data[j];
        auto row = x.row(i);
        // do something with row i
      } else if (clss[j] == "dgCMatrix") {
        arma::sp_mat x = data[j];
        auto row = x.row(i);
        // do something different with row i
      }
    }
  }

  return true;
}

可以在 R 中调用此函数：

data <- list(
  x = Matrix::rsparsematrix(nrow = 1000, ncol = 1000, density = 0.3), 
  y = matrix(1:10000, nrow = 1000, ncol = 10)
  )

clss <- c("dgCMatrix", "matrix")

write_rows(data, clss, 1000)

该函数接收具有相同行数的矩阵或稀疏矩阵列表，并逐行写入这些矩阵，即。首先写入 data 中所有元素的第一行，然后写入所有元素的第二行，依此类推

我的问题是，这条线 arma::sp_mat x = data[i]; 似乎对性能产生了巨大影响，因为我似乎隐式地将列表元素 data[j] 转换为犰狳稀疏矩阵 n次。

我的问题是： 无论如何我可以避免这种情况吗？有没有更有效的解决方案？我试图通过查看 readr 的源代码找到解决方案，因为他们也逐行编写列表元素，但他们也对每一行进行转换（in this line for example，但也许这不是不会影响性能，因为它们处理 SEXPS？

Answer 1

根据说明，结果似乎应该交错每个矩阵的行。您仍然可以这样做，同时避免多次转化。

这是原始代码，经过修改以生成一些实际输出：

// [[Rcpp::export]]
arma::mat write_rows(Rcpp::List data, Rcpp::CharacterVector clss, int nrows, int ncols) {

    int len = data.length();
    arma::mat result(nrows*len, ncols);

    for (int i = 0, k = 0; i < nrows; i++) {
        for (int j = 0; j < len; j++) {
            arma::rowvec r;

            if (clss[j] == "matrix") {
                Rcpp::NumericMatrix x = data[j];
                r = x.row(i);
            }
            else {
                arma::sp_mat x = data[j];
                r = x.row(i);
            }

            result.row(k++) = r;
        }
    }

    return result;
}

以下代码创建一个转换对象向量，然后根据需要从每个对象中提取行。每个矩阵只进行一次转换。我使用包含密集和稀疏垫的结构，因为它比处理联合要简单得多；而且我不想拖入 boost::variant 或需要 C++17。由于我们只想处理 2 类，因此开销很小。

struct Matrix_types {
    arma::mat m;
    arma::sp_mat M;
};

// [[Rcpp::export]]
arma::mat write_rows2(Rcpp::List data, Rcpp::CharacterVector clss, int nrows, int ncols) {

    const int len = data.length();
    std::vector<Matrix_types> matr(len);
    std::vector<bool> is_dense(len);
    arma::mat result(nrows*len, ncols);

    // populate the structs
    for (int j = 0; j < len; j++) {
        is_dense[j] = (clss[j] == "matrix");
        if (is_dense[j]) {
            matr[j].m = Rcpp::as<arma::mat>(data[j]);
        }
        else {
            matr[j].M = Rcpp::as<arma::sp_mat>(data[j]);
        }
    }

    // populate the result
    for (int i = 0, k = 0; i < nrows; i++) {
        for (int j = 0; j < len; j++, k++) {
            if (is_dense[j]) {
                result.row(k) = matr[j].m.row(i);
            }
            else {
                arma::rowvec r(matr[j].M.row(i));
                result.row(k) = r;
            }
        }
    }

    return result;
}

运行关于一些测试数据：

data <- list(
    a=Matrix(1.0, 1000, 1000, sparse=TRUE),
    b=matrix(2.0, 1000, 1000),
    c=Matrix(3.0, 1000, 1000, sparse=TRUE),
    d=matrix(4.0, 1000, 1000)
)

system.time(z <- write_rows(data, sapply(data, class), 1000, 1000))
#   user  system elapsed 
# 185.75   35.04  221.38 

system.time(z2 <- write_rows2(data, sapply(data, class), 1000, 1000))
#   user  system elapsed 
#   4.21    0.05    4.25 

identical(z, z2)
# [1] TRUE

Rcpp 中许多对象的交错结果

Interleaving results from many objects in Rcpp

r

rcpp