稀疏 x 密集矩阵与犰狳相乘出乎意料地慢

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

Answer 1

稀疏矩阵和密集矩阵的存储方式截然不同。 Armadillo 对密集矩阵使用 CMS（column-major 存储），对稀疏矩阵使用 CSC（压缩稀疏列）。来自 Armadillo 的文档：

Mat
mat
cx_mat
Classes for dense matrices, with elements stored in column-major ordering (ie. column by column)

SpMat
sp_mat
sp_cx_mat
Classes for sparse matrices, with elements stored in compressed sparse column (CSC) format

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

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

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

有关存储格式的更多详细信息，请参阅 wikipedia, or netlib。

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

storage_cms = 800MB
storage_csc = 12.04MB

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

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

无论什么算法，它肯定需要多次访问两个输入矩阵，这就解释了为什么将密集矩阵转换为稀疏矩阵（只需要一次访问 800MB 的密集矩阵），然后进行稀疏 x 稀疏product（需要多次访问24MB内存）比dense x sparse和sparse x dense product更高效。

这里还有各种各样的缓存效果，需要了解具体的算法实现和硬件知识（还需要很多时间）才能解释清楚，不过以上是大概的思路。

至于为什么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 则不会。

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

Answer 2

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

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

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

稀疏 x 密集矩阵与犰狳相乘出乎意料地慢

Sparse x dense matrix multiply unexpectedly slow with Armadillo

c++

r

sparse-matrix

armadillo

rcpp