Scipy 稀疏矩阵 - 密集向量乘法性能 - 块与大型矩阵答案

【问题标题】：Scipy Sparse Matrix - Dense Vector Multiplication Performance - Blocks vs Large MatrixScipy 稀疏矩阵 - 密集向量乘法性能 - 块与大型矩阵
【发布时间】：2013-05-20 18:19:35
【问题描述】：

我有许多 scipy 稀疏矩阵（当前为 CSR 格式），我需要将它们与密集的 numpy 1D 向量相乘。向量被称为G：

print G.shape, G.dtype
(2097152,) complex64

每个稀疏矩阵的形状为(16384,2097152)，并且非常稀疏。密度约为4.0e-6。我有一个包含 100 个稀疏矩阵的列表，称为 spmats。

我可以像这样轻松地将每个矩阵与G 相乘：

res = [spmat.dot(G) for spmat in spmats]

这会按预期生成形状为(16384,) 的密集向量列表。

我的应用程序对性能非常关键，因此我尝试了另一种方法，即首先将所有稀疏矩阵连接成一个大的稀疏矩阵，然后只使用一次 dot() 调用，如下所示：

import scipy.sparse as sp
SPMAT = sp.vstack(spmats, format='csr')
RES = SPMAT.dot(G)

这会产生一个长向量RES，其形状为(1638400,)，并且是上面res 中所有结果向量的串联版本，正如预期的那样。我检查了结果是否相同。

也许我完全错了，但我预计第二种情况应该比第一种更快，因为 numpy 调用、内存分配、python 对象的创建、python 循环等要少得多。我不在乎关于连接稀疏矩阵所需的时间，只有计算结果的时间。然而，根据%timeit：

%timeit res = [spmat.dot(G) for spmat in spmats]
10 loops, best of 3: 91.5 ms per loop
%timeit RES = SPMAT.dot(G)
1 loops, best of 3: 389 ms per loop

我已经检查过我在这两个操作中都没有耗尽内存，并且似乎没有发生任何可疑的事情。我疯了，还是这真的很奇怪？这是否意味着所有稀疏矩阵向量乘积都应该以块的形式完成，一次几行，以使它们更快？据我了解，具有密集向量的稀疏矩阵乘法时间应该与非零元素的数量成线性关系，这在上述两种情况下都没有变化。是什么造成了如此大的不同？

我正在使用 EPD7.3 的具有 4GB 内存的单核 linux 机器上运行

编辑：

这是一个为我重现问题的小例子：

import scipy.sparse as sp
import numpy as n

G = n.random.rand(128**3) + 1.0j*n.random.rand(128**3)

spmats = [sp.rand (128**2, 128**3, density = 4e-6, format = 'csr', dtype=float64) for i in range(100)]
SPMAT = sp.vstack(spmats, format='csr')

%timeit res = [spmat.dot(G) for spmat in spmats]
%timeit RES = SPMAT.dot(G)

我明白了：

1 loops, best of 3: 704 ms per loop
1 loops, best of 3: 1.34 s per loop

这种情况下的性能差异没有我自己的具有某种结构的稀疏矩阵（可能是因为缓存）那么大，但连接矩阵仍然更糟。

我已尝试使用 scipy 10.1 和 12.0。

【问题讨论】：

我无法重现：单点积对我来说快 5 倍。由于我仔细尝试过您所描述的操作，您能否发布一个最小的工作示例以确保我们正在做同样的事情？
@jorgeca 感谢您花时间尝试重现问题。我刚刚用我正在做的工作示例编辑了我的问题。
谢谢。我无法重现您的结果（在 scipy 0.12 中），但对我来说，当 G 具有 dtype=np.complex64 时，列表理解要慢 5 倍（！），正如您所说，当 G 具有 dtype=np.complex128 时，这两种方法都同样快.
我无法使用任何 dtype 重现此问题：在我的机器上（使用 Scipy 0.12.0），对于 float32、float64、complex64、complex128，两个循环的速度大致相同。像这样会导致与机器相关的速度差异的一件事是处理器缓存效率。但我不清楚为什么它会影响这个特殊情况。相关内循环为here； Python 开销可能并不重要。另请注意，除 G.dtype==complex128 之外的其他数据类型需要将 spmat 数据额外强制转换为 complex128。
就我而言，我猜 complex64 的减速源于额外的强制（spmats[0].dot(G) 对于 complex128 需要 4.65 毫秒，对于 complex64 需要 22.7 毫秒）。如果我将 128 更改为较小的值，则差异消失，然后增长（70 慢 20 倍）并最终稳定在慢 5-6 倍左右。

标签： python performance numpy scipy sparse-matrix

【解决方案1】：

我还没有找到问题中提到的奇怪行为的原因，但是我找到了一种显着加快计算速度的方法，这可能对其他人有用。

因为在我的特定情况下，我正在计算 float32 稀疏矩阵和 complex64 密集向量的乘积，所以我可以分别将实部和虚部相乘。这为我提供了 4 倍的加速。

SPMAT.shape == (16384000, 2097152) 需要 2.35 秒：

RES = SPMAT.dot(G)

虽然这只需要 541 毫秒：

RES = n.zeros((SPMAT.shape[0],),dtype=complex64)
RES.real = SPMAT.dot(G.real); RES.imag = SPMAT.dot(G.imag)

结果是一样的。我认为n.zeros 预分配可能不是必需的，但我不知道该怎么做。

【讨论】：