使用 R 和 Rcpp，如何将两个稀疏 Matrix::csr/csc 格式的矩阵相乘？

如何解决使用 R 和 Rcpp，如何将两个稀疏 Matrix::csr/csc 格式的矩阵相乘？

以下代码按预期工作：

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

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

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

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

解决方法

dgRMatrix 交叉积函数尚未实现是有原因的，事实上，不应实现它们，否则它们会导致不良做法。

处理稀疏矩阵时有一些性能注意事项：

根据主要边缘方向访问边缘视图的效率非常低。例如，dgRMatrix 中的列迭代器和 dgCMatrix 中的行迭代器需要遍历矩阵的几乎所有元素以找到仅在该列或行中的元素。请参阅此 Rcpp gallery post 以获得更多启发。
矩阵叉积只是所有列组合之间的点积。这意味着在 dgRMatrix 中使用列迭代器（与在 dgCMatrix 中使用列迭代器）的代价乘以列组合数。
R 中的交叉乘积函数经过高度优化，并且（根据我的经验）并不比 Eigen、Armadillo 和等效的 STL 变体快得多。它们是并行化的，Matrix 包充分利用了这些优化算法。我已经使用 Rcpp 结构编写了 C++ 并行化 STL 跨产品变体，但我没有看到性能有任何提高。
如果您真的要走这条路，请查看我关于 Rcpp 中稀疏矩阵结构的 Rcpp gallery 帖子。如果需要考虑内存，这比 Eigen 和 Armadillo 稀疏矩阵更可取，因为 Eigen 和 Armadillo 执行深度复制，而不是对内存中已存在的 R 对象的引用。
在 1% 的密度下，行迭代器的低效率将高于 5% 或 10% 的密度。我的大部分测试都以 5% 的密度进行，通常行迭代器的二元运算比列迭代器花费的时间长 5-10 倍。

可能有些应用程序行优先排序（即，请参阅 Dmitry Selivanov 在 CSR 矩阵和 irlba svd 上的工作），但这绝对不是其中之一，事实上，这么多，所以你最好在-place 转换以获得 CSC 矩阵。

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