如何基于Eigen实现R语言csr/csc格式稀疏矩阵的C++乘法

以下为可正常运行的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矩阵类实现，使用时可参考Eigen官方文档的相关说明。

在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自带的转置运算性能表现较好
• 基于Eigen实现的矩阵乘法具备一定性能优势

可借助R的Matrix包将普通矩阵转换为稀疏矩阵，示例代码如下：

library(Matrix)
data <- rnorm(1e6)
zero_index <- sample(1e6)[1:9e5]
data[zero_index] <- 0
A <- matrix(data, ncol=1000)

# 转换为CSR格式稀疏矩阵（dgRMatrix）
A.csr <- as(A, "dgRMatrix")
B.csr <- t(A.csr)

# 转换为CSC格式稀疏矩阵（dgCMatrix）
A.csc <- as(A, "dgCMatrix")
B.csc <- t(A.csc)

注：当前R语言中，A.csr %*% B.csr（CSR格式稀疏矩阵乘法）尚未实现，仅A.csc %*% B.csc（CSC格式稀疏矩阵乘法）可正常运行。

• 在C++层基于Eigen实现CSR格式稀疏矩阵（对应R的dgRMatrix）的乘法运算
• 尽可能避免不必要的类型转换，减少额外内存开销
• 后续需通过microbenchmark测试不同实现方案的性能，明确不同矩阵规模下的最高效方案，最终需适配约1%稀疏度、规模达500万行×500万列的超大规模矩阵运算

内容的提问来源于stack exchange，提问作者mshaffer