来自稀疏矩阵的Rcpp子矩阵

Question

我正在尝试将vec乘以非常大的稀疏矩阵的子集（如脚本所示），但是使用sourceCpp时它无法编译，报告为error: no matching function for call to ‘arma::SpMat<double>::submat(arma::uvec&, arma::uvec&)，如果有人可以帮我一个忙。

#include <RcppArmadillo.h>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
double myf(sp_mat X, vec g, uvec xi){
    double u = g(xi).t() * X.submat(xi, xi) * g(xi);
    return u;
}

Answer 1

所以，正如@RalfStubner提到的，matrix access for sparse matrices is continuous only。也就是说，由于使用了相同的索引，因此对于实际的稀疏矩阵采用的访问方法是对称的。因此，在这种情况下，有必要恢复为(x,y)的标准元素访问器。结果，可以通过单个循环来完成总和减少。

(x,y)

另一种方法，但可能代价更高，将是提取稀疏矩阵的对角线并将其转换为密集向量。从那里开始按元素进行乘法和求和。

#include <RcppArmadillo.h>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
double submat_multiply(const arma::sp_mat& X, 
                       const arma::vec& g, const arma::uvec& xi){

  // Add an assertion
  if(X.n_rows != g.n_elem) {
    Rcpp::stop("Mismatched row and column dimensions of X (%s) and g (%s).",
               X.n_rows, g.n_elem);
  }

  // Reduction
  double summed = 0; 

  for (unsigned int i = 0; i < xi.n_elem; ++i) {
    // Retrieve indexing element
    arma::uword index_at_i = xi(i);
    // Add components together
    summed += g(index_at_i) * X(index_at_i, index_at_i) * g(index_at_i);
  }

  // Return result
  return summed;
}

测试代码：

#include <RcppArmadillo.h>

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
double submat_multiply_v2(const arma::sp_mat& X, 
                          const arma::vec& g, const arma::uvec& xi){

  // Add an assertion
  if(X.n_rows != g.n_elem) {
    Rcpp::stop("Mismatched row and column dimensions of X (%s) and g (%s).",
               X.n_rows, g.n_elem);
  }

  // Copy sparse diagonal to dense vector
  arma::vec x_diag(X.diag());
  // Obtain the subset
  arma::vec g_sub = g.elem(xi); 

  // Perform element-wise multiplication and then sum.
  double summed = arma::sum(g_sub % x_diag.elem(xi) % g_sub);

  // Return result
  return summed;
}