使用 Rcpp 通过引用传递犰狳稀疏矩阵

Posted 2023-02-16

【中文标题】使用 Rcpp 通过引用传递犰狳稀疏矩阵

【英文标题】：Passing an Armadillo sparse matrix by reference with Rcpp

【发布时间】：2018-04-04 14:49:23

【问题描述】：

此问题与this 和this 有关。这里的不同之处在于，我传递的不是NumericVector 或NumericMatrix 之类的Rcpp 类型，而是arma::sp_mat。

有没有办法将sp_mat 传递给 C++，修改它的值，并让更改显示在 R 中的原始对象中？

这可以通过NumericMatrix 来完成，例如：

cppFunction("void frob(NumericMatrix& x)

    for(NumericMatrix::iterator it = x.begin(); it != x.end(); ++it)
    
        if(*it != 0) *it = *it + 5;
    
")

M <- Matrix(0, 5, 1, sparse=TRUE)
M[1] <- 1.2345

m <- as.matrix(M)

frob(m)
m
       #[,1]
#[1,] 6.2345
#[2,] 0.0000
#[3,] 0.0000
#[4,] 0.0000
#[5,] 0.0000

同样的技术适用于arma::mat 密集矩阵。但是对于稀疏矩阵，它不起作用：

cppFunction("void frob2(arma::sp_mat& x)

    for(arma::sp_mat::iterator it = x.begin(); it != x.end(); ++it)
    
        *it = *it + 5;
    
", depends="RcppArmadillo")

frob2(M)
M
#5 x 1 sparse Matrix of class "dgCMatrix"

#[1,] 1.2345
#[2,] .     
#[3,] .     
#[4,] .     
#[5,] .     

【问题讨论】：

嗯，这可能是去年夏天对稀疏进口商所做的工作的回归。

@coatless 很有趣，所以这不是预期的行为吗？有没有办法修改矩阵而不进行不必要的复制？

不确定我们曾经有过引用语义。毕竟，稀疏矩阵没有原生类型，所以需要在某个地方复制一些东西......

【参考方案1】：

不幸的是 there is no auxiliary memory constructor 用于犰狳中的稀疏矩阵。

但是，您可以使用指向 R 对象的指针在 C++ 中构造类似稀疏矩阵的结构。 Here is example:

template< typename T>
class MappedCSC 
public:
  MappedCSC();
  MappedCSC(std::uint32_t n_rows,
            std::uint32_t n_cols,
            size_t nnz,
            std::uint32_t * row_indices,
            std::uint32_t * col_ptrs,
            T * values):
    n_rows(n_rows), n_cols(n_cols), nnz(nnz), row_indices(row_indices), col_ptrs(col_ptrs), values(values) ;
  const std::uint32_t n_rows;
  const std::uint32_t n_cols;
  const size_t nnz;
  const std::uint32_t * row_indices;
  const std::uint32_t * col_ptrs;
  T * values;
;

using dMappedCSC = MappedCSC<double>;

Here is how you can extract it:

dMappedCSC extract_mapped_csc(Rcpp::S4 input) 
  Rcpp::IntegerVector dim = input.slot("Dim");
  Rcpp::NumericVector values = input.slot("x");
  uint32_t nrows = dim[0];
  uint32_t ncols = dim[1];
  Rcpp::IntegerVector row_indices = input.slot("i");
  Rcpp::IntegerVector col_ptrs = input.slot("p");
  return dMappedCSC(nrows, ncols, values.length(), (uint32_t *)row_indices.begin(), (uint32_t *)col_ptrs.begin(), values.begin());

以及here is example 关于如何逐列迭代：

Rcpp::NumericMatrix dense_csc_prod(const Rcpp::NumericMatrix &x_r, const Rcpp::S4 &y_csc_r) 
  const arma::dmat x = arma::dmat((double *)&x_r[0], x_r.nrow(), x_r.ncol(), false, false);
  const dMappedCSC y_csc = extract_mapped_csc(y_csc_r);
  Rcpp::NumericMatrix res(x.n_rows, y_csc.n_cols);
  arma::dmat res_arma_map = arma::dmat(res.begin(), res.nrow(), res.ncol(), false, false);
  for (uint32_t i = 0; i < y_csc.n_cols; i++) 
    const uint32_t p1 = y_csc.col_ptrs[i];
    const uint32_t p2 = y_csc.col_ptrs[i + 1];
    // mapped indices are uint32_t, but arma only allows indices be uvec = vec<uword> = vec<size_t>
    // so we need to construct these indices by copying from uint32_t to uword
    const arma::Col<uint32_t> idx_temp = arma::Col<uint32_t>(&y_csc.row_indices[p1], p2 - p1);
    const arma::uvec idx = arma::conv_to<arma::uvec>::from(idx_temp);
    const arma::colvec y_csc_col = arma::colvec(&y_csc.values[p1], p2 - p1, false, false);
    res_arma_map.col(i) = x.cols(idx) * y_csc_col;
  
  return res;

【讨论】：

非常好。把它变成 Rcpp Gallery 的帖子怎么样？