稀疏 x 密集矩阵乘以 Armadillo 出乎意料地慢

Posted 2023-02-16

【中文标题】稀疏 x 密集矩阵乘以 Armadillo 出乎意料地慢

【英文标题】：Sparse x dense matrix multiply unexpectedly slow with Armadillo

【发布时间】：2018-04-11 09:13:37

【问题描述】：

这是我刚刚遇到的。出于某种原因，在犰狳中将密集矩阵与稀疏矩阵相乘比将稀疏矩阵和密集矩阵相乘（即颠倒顺序）要慢得多。

// [[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 倍的速度差异似乎相当大。

【问题讨论】：

我不会使用稀疏矩阵作为结果，但这并不会真正影响您示例中的性能。

【参考方案1】：

稀疏矩阵和稠密矩阵的存储方式截然不同。 Armadillo 使用 CMS（主要列存储）存储密集矩阵，使用 CSC（压缩稀疏列）存储稀疏矩阵。来自犰狳的文档：

垫子垫子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 位 indeces，在您的情况下，这意味着：

storage_cms = 800MB
storage_csc = 12.04MB

因此，当您乘以稀疏 x 密集（或密集 x 稀疏）矩阵时，您将访问约 812MB 的内存，而在乘以稀疏 x 稀疏矩阵时仅访问约 24MB 的内存。

请注意，这不包括您写入结果的内存，这可能是很大一部分（在两种情况下都高达 ~800MB），但我不太熟悉犰狳以及它用于矩阵的算法乘法，所以不能确切地说它是如何存储中间结果的。

不管是什么算法，它肯定需要多次访问两个输入矩阵，这就解释了为什么将密集矩阵转换为稀疏矩阵（只需要一次访问 800MB 的密集矩阵），然后做一个稀疏 x 稀疏乘积 (需要多次访问 24MB 内存）比密集 x 稀疏和稀疏 x 密集产品更有效。

这里还有各种各样的缓存效果，这需要了解算法的确切实现和硬件（以及大量时间）才能正确解释，但以上是大致思路。

至于为什么dense x sparse比sparse x dense快，是因为稀​​疏矩阵的CSC存储格式。如scipy's documentation 中所述，CSC 格式对于列切片是有效的，而对于行切片则很慢。 dense x sparse 乘法算法需要对稀疏矩阵进行列切片，sparse x dense 需要对稀疏矩阵进行行切片。请注意，如果犰狳使用 CSR 而不是 CSC，sparse x dense 会很有效，而dense x sparse 则不会。

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

【讨论】：

然而，a %*% b_den 在 R 中速度很快。

【参考方案2】：

这应该会在即将发布的Armadillo 8.500 中得到解决，该问题将很快包含在 RcppArmadillo 0.8.5 Real 中。具体来说：

稀疏矩阵转置要快得多 (sparse x dense) 重新实现为 ((dense^T) x (sparse^T))^T，利用相对快速的 (dense x sparse) 代码

当我测试它时，所用时间从大约 500 秒下降到大约 18 秒，这与其他时间相当。