c++ - 使用 R 和 Rcpp,如何将两个稀疏 Matrix::csr/csc 格式的矩阵相乘?
问题描述
以下代码按预期工作:
矩阵.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 万行和列...
解决方案
dgRMatrix 叉积函数尚未实现是有原因的,事实上,它们不应该被实现,否则它们会导致不好的做法。
使用稀疏矩阵时有一些性能注意事项:
- 针对主要边缘方向访问边缘观点是非常低效的。例如,dgRMatrix 中的列迭代器和 dgCMatrix 中的行迭代器需要遍历矩阵的几乎所有元素,才能仅在该列或行中找到元素。请参阅此Rcpp 画廊帖子以获得更多启发。
- 矩阵叉积只是所有列组合之间的点积。这意味着在 dgRMatrix 中使用列迭代器(相对于 dgCMatrix 中的列迭代器)的损失乘以列组合的数量。
- R 中的叉积函数经过高度优化,并且(根据我的经验)并不比 Eigen、Armadillo 和等效的 STL 变体快得多。它们是并行化的,Matrix 包充分利用了这些优化算法。我已经使用 Rcpp 结构编写了 C++ 并行化 STL 交叉产品变体,但我看不到性能有任何提高。
- 如果您真的要走这条路,请查看我在 Rcpp 中的稀疏矩阵结构上的Rcpp 画廊帖子。如果内存是一个问题,这比 Eigen 和 Armadillo 稀疏矩阵更受欢迎,因为 Eigen 和 Armadillo 执行深度复制而不是对内存中已经存在的 R 对象的引用。
- 在 1% 的密度下,行迭代器的低效率将大于 5% 或 10% 的密度。我以 5% 的密度进行大部分测试,通常二进制操作对于行迭代器比列迭代器花费的时间长 5-10 倍。
可能存在行优先排序的应用程序(即参见 Dmitry Selivanov 在 CSR 矩阵和 irlba svd 上的工作),但事实上,这绝对不是其中之一,所以你最好进行就地转换得到一个 CSC 矩阵。
tl; dr:行主要矩阵中的按列交叉积是低效率的最后通牒。
推荐阅读
- amazon-web-services - 找不到 AWS 身份池的规则设置
- javascript - 如果宽度溢出,如何使reactjs将图例重新绘制到下一行?
- windows - 禁用 Windows Server 2016 的自动 Windows 更新,使用 powershell 更改注册表项
- c - gst_init (&argc, &argv) 与 gst_value_table 崩溃为 nullptr
- mysql - 查询 mysql 以获取 2 个表,其中包含每月已售商品的摘要,包括未售出的商品
- python - 无法在python中导入文件
- python - 在 Python 中添加缺失的句点
- neo4j - Neo4j:使用 Cypher 创建具有唯一节点的图
- python - 我在从 Python 中的另一个文件导入列表时遇到困难
- vue.js - 无法读取 vuex 中 Store._callee 处未定义的属性“getProduct”(操作名称)