稀疏矩阵乘法的执行时间

Question

我运行ning 在产品扩散研究中的一个关键步骤涉及稀疏矩阵与向量的乘法。下面的代码重现了 Rcpp 和 R 版本的代码。

#define ARMA_64BIT_WORD
#include <RcppArmadillo.h>
//[[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]
arma::vec multiply(arma::sp_mat A, arma::vec nodes_status) {

  return A * nodes_status;
}


/***R

library(rbenchmark)
library(igraph)

g <- sample_smallworld(dim = 1, size = 1e5, nei = 12, p = 0.65)

nodes <- V(g)
A <- get.adjacency(g)

nodes.status <- sample(x = c(0, 1), 
                       size = length(nodes), 
                       replace = TRUE, 
                       prob = c(0.9, 0.1))

multiplyR <- function(A, n) {return(as.vector(A %*% n))}

sum(multiply(A, nodes.status) == multiplyR(A, nodes.status)) == 1e5 # check if the results are the same

benchmark(multiply(A, nodes.status), multiplyR(A, nodes.status), order = "relative")[, 1:4]

*/

当我运行这段代码时，两个函数的答案相符，但执行时间如下：

                    test     replications elapsed relative
2 multiplyR(A, nodes.status)          100    1.30    1.000
1  multiply(A, nodes.status)          100    3.66    2.815

这里的代码有什么问题？对于我的情况，Rcpp 中是否有更有效的乘法习语？

Answer 1

没有。此处记录的问题有两个：

将 R 中的稀疏矩阵转换为 RcppArmadillo 对象所需的转换时间大于 Matrix 对 C 例程的使用。
复制与参考创建的成本。

关于 2.，arma::sp_mat 的构造使用副本，因为它没有引用（例如 &）作为后缀。特别注意：

#define ARMA_64BIT_WORD
#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::vec sp_multiply_copy(arma::sp_mat A, arma::vec nodes_status) {
  return A * nodes_status;
}

// [[Rcpp::export]]
arma::vec sp_multiply_ref(const arma::sp_mat& A, const arma::vec& nodes_status) {
  return A * nodes_status;
}

从这里可以看出两者之间的细微性能差异：

benchmark(sp_multiply_copy(A, nodes.status),
          sp_multiply_ref(A, nodes.status),
          multiplyR(A, nodes.status), order = "relative")[, 1:4]

#                                test replications elapsed relative
# 3        multiplyR(A, nodes.status)          100   1.240    1.000
# 2  sp_multiply_ref(A, nodes.status)          100   2.766    2.231
# 1 sp_multiply_copy(A, nodes.status)          100   3.141    2.533

话虽这么说，我们 return 第一点：sparse 矩阵的 Matrix 函数是高度优化的，直接使用 C. 可以在 R/products.R, src/Csparse.c, src/dtrMatrix.c 查看上述例程的示例。因此，Matrix 操作将更加高效。

现在，不是说C++ 无法获得速度。特别是，如果在 C++ 中重复使用矩阵对象与传递引用的实例化（例如 &）相乘，那么它应该比调用更快Matrix 的乘法例程。

稀疏矩阵乘法的执行时间

Execution time of sparse matrix multiplication

r

armadillo

rcpp