Armadillo中的稀疏矩阵与密集矩阵相乘速度意外缓慢

Question

我刚刚发现了一件事情。在Armadillo中，将密集矩阵乘以稀疏矩阵的速度比将稀疏矩阵乘以密集矩阵慢得多（即颠倒顺序）。

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::sp_mat mult_sp_den_to_sp(arma::sp_mat& a, arma::mat& b)
{
    // sparse x dense -> sparse
    arma::sp_mat result(a * b);
    return result;
}

// [[Rcpp::export]]
arma::sp_mat mult_den_sp_to_sp(arma::mat& a, arma::sp_mat& b)
{
    // dense x sparse -> sparse
    arma::sp_mat result(a * b);
    return result;
}

我正在使用RcppArmadillo包将Arma与R进行接口交互；RcppArmadillo.h包含armadillo。以下是在R中对一些相当大的矩阵进行计时的结果：

set.seed(98765)
# 10000 x 10000 sparse matrices, 99% sparse
a <- rsparsematrix(1e4, 1e4, 0.01, rand.x=function(n) rpois(n, 1) + 1)
b <- rsparsematrix(1e4, 1e4, 0.01, rand.x=function(n) rpois(n, 1) + 1)

# dense copies
a_den <- as.matrix(a)
b_den <- as.matrix(b)

system.time(mult_sp_den_to_sp(a, b_den))
#   user  system elapsed 
# 508.66    0.79  509.95 

system.time(mult_den_sp_to_sp(a_den, b))
#   user  system elapsed 
#  13.52    0.74   14.29 

因此，第一个乘法比第二个乘法慢了约35倍（所有时间以秒为单位）。

有趣的是，如果我只是暂时地制作一个稀疏矩阵的副本，性能会大大提高：

// [[Rcpp::export]]
arma::sp_mat mult_sp_den_to_sp2(arma::sp_mat& a, arma::mat& b)
{
    // sparse x dense -> sparse
    // copy dense to sparse, then multiply
    arma::sp_mat temp(b);
    arma::sp_mat result(a * temp);
    return result;
}

system.time(mult_sp_den_to_sp2(a, b_den))
#   user  system elapsed 
#   5.45    0.41    5.86 

这是预期的行为吗？我知道在稀疏矩阵中，你做事情的确切方式对代码的效率有很大的影响，比密集矩阵更加明显。但是速度相差35倍似乎有些太大了。

- Hong Ooi

---

我不会将稀疏矩阵用作结果，但这并不会在您的示例中影响性能。 - F. Privé

Answer 1

这个问题应该会在即将发布的Armadillo 8.500中得到解决，它将被包裹在RcppArmadillo 0.8.5中。具体而言：

• 稀疏矩阵转置速度更快了
• (sparse x dense) 重新实现为 ((dense^T) x (sparse^T))^T，利用了相对较快的(dense x sparse)代码

当我测试时，所需的时间从约500秒降至约18秒，与其他计时相当。

- Hong Ooi

Answer 2

稀疏矩阵和密集矩阵的存储方式非常不同。Armadillo 对于密集矩阵使用CMS（列主存储），对于稀疏矩阵使用CSC（压缩稀疏列）。根据Armadillo文档：

Mat、mat和cx_mat类用于密集矩阵，元素按列主序存储（即逐列存储）。

SpMat、sp_mat和sp_cx_mat类用于稀疏矩阵，元素以压缩稀疏列（CSC）格式存储。

我们首先要了解的是每种格式需要多少存储空间。

根据给定的量element_size（单精度为4字节，双精度为8字节）、index_size（如果使用32位整数，则为4字节；如果使用64位整数，则为8字节）、num_rows（矩阵的行数）、num_cols（矩阵的列数）和num_nnz（矩阵的非零元素数），以下公式给出了每种格式的存储空间：

storage_cms = num_rows * num_cols * element_size
storage_csc = num_nnz * element_size + num_nnz * index_size + num_cols * index_size

关于存储格式的更多细节，请参见wikipedia或netlib。

假设双精度和32位索引，在您的情况下意味着：

storage_cms = 800MB
storage_csc = 12.04MB

当你在进行稀疏矩阵与密集矩阵（或者密集矩阵与稀疏矩阵）的乘法时，你会访问大约812MB的内存，而当你进行稀疏矩阵与稀疏矩阵的乘法时，你只会访问大约24MB的内存。

需要注意的是，这并不包括你写入结果的内存，这可能是一个相当大的部分（在两种情况下都高达800MB），但我不太熟悉Armadillo以及它用于矩阵乘法的算法，因此无法确定它如何存储中间结果。

无论使用什么样的算法，它肯定需要多次访问输入矩阵，这就解释了为什么将密集矩阵转换为稀疏矩阵（只需要一次访问800MB的密集矩阵），然后进行稀疏矩阵乘法（需要多次访问24MB的内存）比进行密集矩阵与稀疏矩阵或者稀疏矩阵与密集矩阵乘法更有效率。

这里还有各种缓存效应，需要知道算法的具体实现和硬件才能充分解释（需要花费大量时间），但以上是一般的概念。

关于为什么 dense x sparse 比 sparse x dense 快，这是因为稀疏矩阵的 CSC 存储格式。正如 Scipy 文档 中所述，CSC 格式对列切片很有效，但对行切片则较慢。 dense x sparse 的乘法算法需要对稀疏矩阵进行列切片，而 sparse x dense 则需要对稀疏矩阵进行行切片。请注意，如果 Armadillo 使用 CSR 而不是 CSC，则 sparse x dense 将会更有效，而 dense x sparse 则不会。

我知道这并不是关于您看到的所有性能影响的完整答案，但应该可以让您大致了解正在发生的事情。要进行适当的分析需要更多的时间和精力，必须包括算法的具体实现以及运行它的硬件信息。

- gflegar

---

然而，在R中，a %*% b_den非常快速。 - F. Privé