7

我有许多 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。

4

1 回答 1

4

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

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

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

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预分配可能不是必需的,但我不知道该怎么做。

于 2013-05-25T20:44:14.610 回答