如何将两个稀疏矩阵相乘:csr/csc格式

Question

以下代码按预期工作：

matrix.cpp

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
// [[Rcpp::export]]
SEXP eigenMatTrans(Eigen::MatrixXd A){
Eigen::MatrixXd C = A.transpose();
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMatMult(Eigen::MatrixXd A, Eigen::MatrixXd B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}
// [[Rcpp::export]]
SEXP eigenMapMatMult(const Eigen::Map<Eigen::MatrixXd> A, Eigen::Map<Eigen::MatrixXd> B){
Eigen::MatrixXd C = A * B;
return Rcpp::wrap(C);
}

这是对矩阵使用C++特征类，请参见https://eigen.tuxfamily.org/dox

在R中，我可以访问这些函数。

library(Rcpp);
Rcpp::sourceCpp('matrix.cpp');  
A <- matrix(rnorm(10000), 100, 100);
B <- matrix(rnorm(10000), 100, 100);
library(microbenchmark);
microbenchmark(eigenMatTrans(A), t(A), A%*%B, eigenMatMult(A, B), eigenMapMatMult(A, B))

这表明R在诉诸(转置(时表现得相当好。乘法与本征有一些优点。

使用矩阵库，我可以将正规矩阵转换为稀疏矩阵。

示例来自https://cmdlinetips.com/2019/05/introduction-to-sparse-matrices-in-r/

library(Matrix);
data<- rnorm(1e6)
zero_index <- sample(1e6)[1:9e5]
data[zero_index] <- 0
A = matrix(data, ncol=1000)
A.csr = as(A, "dgRMatrix");
B.csr = t(A.csr);
A.csc = as(A, "dgCMatrix");
B.csc = t(A.csc);

所以，如果我想用本征将A.csr乘以B.csr，如何在C++中做到这一点？如果不需要的话，我不想转换类型。这是内存大小的问题。

A.csr %*% B.csr尚未实现。A.csc %*% B.csc正在工作。

我想对不同的选项进行微基准测试，看看矩阵大小将如何最有效。最后，我将得到一个大约1%稀疏的矩阵，它有500万行和列。。。

Answer 1

dgRMatrix交叉产品函数尚未实现是有原因的，事实上，它们不应该实现，因为否则它们会导致不良实践。

在处理稀疏矩阵时，有一些性能考虑因素：

• 根据主要的边缘方向访问边缘视图是非常低效的。例如，dgRMatrix中的列迭代器和dgCMatrix中的行迭代器需要循环遍历矩阵的几乎所有元素，才能找到该列或行中的元素。请参阅Rcpp画廊的这篇文章，以获得更多启示
• 矩阵叉积只是列的所有组合之间的点积。这意味着在dgRMatrix中使用列迭代器的代价(与dgCMatrix中的列迭代者相比(乘以列组合的数量
• R中的叉积函数是高度优化的，并且(根据我的经验(没有比Eigen、Armadillo、等效STL变体快得多。它们是并行化的，矩阵包充分利用了这些优化算法。我已经使用Rcpp结构编写了C++并行化的STL跨产品变体，我没有看到性能的任何提高
• 如果你真的要走这条路，请查看我在Rcpp中关于稀疏矩阵结构的Rcpp画廊帖子。如果内存是一个问题，这将优先于特征矩阵和Armadillo稀疏矩阵，因为特征矩阵和Armadillo执行深度复制，而不是引用内存中已经存在的R对象
• 在1%的密度下，行迭代器的低效率将大于在例如5%或10%的密度下。我的大多数测试都是在5%的密度下进行的，通常行迭代器的二进制操作比列迭代器长5-10倍

可能存在行主排序大放异彩的应用程序(即，请参阅Dmitry Selivanov关于CSR矩阵和irlba-svd的工作(，但这绝对不是其中之一，事实上，这太多了，所以您最好进行就地转换以获得CSC矩阵。

tl；dr：行主矩阵中的列叉积是低效率的最后通牒。