如何对稀疏矩阵Matrix::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++的eigen矩阵类，详见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在重新排序（转置）方面表现得相当不错。使用特征值进行乘法运算具有一定的优势。

使用Matrix库，我可以将普通矩阵转换为稀疏矩阵。

示例来自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);

如果我想使用eigen将A.csr乘以B.csr，该如何在C++中实现？如果不必要转换类型，我不想这样做。这是一个内存大小的问题。

A.csr %*% B.csr尚未实现。 A.csc %*% B.csc已经可以工作。

我想对不同的选项进行微基准测试，并查看矩阵大小如何最有效。最终，我将拥有一个稀疏度约为1％，具有500万行和列的矩阵...

- mshaffer

---

也许 https://github.com/dselivanov/MatrixCSR 可以提供一些关于如何做这个的提示。 - user20650

Answer 1

“dgRMatrix”交叉乘积功能尚未实现并且不应该实现，因为这会引起不良的实践。在处理稀疏矩阵时需要考虑一些性能问题。”

• 访问主要边缘方向相反的边缘视图非常低效。例如，在dgRMatrix中的列迭代器和在dgCMatrix中的行迭代器需要循环遍历矩阵的几乎所有元素才能找到该列或行中的元素。请参见此Rcpp gallery post获取额外的启示。
• 矩阵交叉积只是所有列组合之间的点积。这意味着在dgRMatrix中使用列迭代器（与在dgCMatrix中使用列迭代器）的罚款将乘以列组合数。
• R中的交叉积函数高度优化，并且（据我的经验）与Eigen，Armadillo，等效的STL变体相比没有显着更快。它们是并行化的，并且Matrix包充分利用了这些优化算法。我使用Rcpp结构编写了C ++并行化的STL交叉积变量，但我没有看到任何性能提升。
• 如果您真的要走这条路，请查看我的Rcpp gallery中有关Rcpp中稀疏矩阵结构的帖子。如果内存是一个问题，则应首选Eigen和Armadillo稀疏矩阵，因为Eigen和Armadillo执行“深拷贝”而不是对已存在于内存中的R对象的引用。
• 在1％的密度下，行迭代器的低效性将大于例如5或10％的密度。我在5％密度下进行大多数测试，并且通常行迭代器的二进制操作比列迭代器慢5-10倍。

可能有一些应用程序适合使用行优先顺序（例如参见Dmitry Selivanov关于CSR矩阵和irlba svd的工作），但这绝对不是其中之一，事实上，你最好进行原地转换以获得CSC矩阵。

简而言之：在行优先矩阵中进行列向量叉积是极其低效的。

- zdebruine

---

非常好的帖子，谢谢——已经相应地点赞了。 - Dirk Eddelbuettel