Passing an Armadillo sparse matrix by reference with Rcpp

Question

This question is related to this and this. The difference here is that I'm not passing an Rcpp type like NumericVector or NumericMatrix, but an arma::sp_mat.

Is there any way to pass an sp_mat to C++, modify its values, and have the changes show up in the original object in R?

This can be done with a NumericMatrix, for example:

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

The same technique works for an arma::mat dense matrix. But for a sparse matrix, it doesn't work:

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,] .     

Answer 1

Unfortunately there is no auxiliary memory constructor for sparse matrices in Armadillo.

However you can construct sparse matrix like structure in C++ using pointers to R objects. 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());
}

And here is example on how to iterate column by column:

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