Scipy稀疏矩阵与密集向量的乘法性能 - 块状矩阵 vs 大型矩阵

Question

我有一些scipy稀疏矩阵（目前在CSR格式中），需要与密集的numpy 1D向量相乘。 这个向量叫做G：

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

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

我可以像这样轻松地将每个矩阵与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

我已经检查过，无论是在任何操作中，我都没有内存不足的情况，并且似乎没有任何可疑的事情发生。这是真的奇怪吗？这是否意味着所有的稀疏矩阵向量积都应该分块进行，每次只处理几行，以使它们更快？

就我所知，稀疏矩阵与密集向量相乘的时间应该是非零元素个数的线性函数，在上述两种情况下未发生改变。是什么导致了这样的差异？

我正在单核Linux机器上运行，拥有4GB的RAM，使用EPD7.3。

编辑：

以下是一个可以重现问题的小例子。

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。

- ali_p

---

我无法复现那个问题：使用单点乘积的速度比你所描述的方法快了5倍。既然我已经仔细尝试了你所述的方法，能否请您提供最小的可工作示例以确保我们执行相同的操作？ - jorgeca

@jorgeca 感谢您抽出时间尝试重现问题。我刚刚编辑了我的问题，并提供了一个可工作的示例。 - ali_p

谢谢。我无法在scipy 0.12中重现您的结果，但是对于我来说，当您所说的G具有dtype = np.complex64时，列表推导式要慢5倍（！），而当G具有dtype = np.complex128时，两种方法的速度都相同。 - jorgeca

我无法在我的机器上重现这个问题（使用Scipy 0.12.0），对于float32、float64、complex64和complex128，两个循环的速度大致相同。导致这种机器相关的速度差异的原因之一是处理器缓存效率。但我不清楚为什么它会影响到这个特定的情况。相关的内部循环在这里；Python开销可能不重要。还要注意，除了G.dtype==complex128之外的其他数据类型需要将spmat数据额外强制转换为complex128。 - pv.

在我的情况下，我猜测complex64的减速源于额外的强制转换（spmats [0] .dot（G）对于complex128需要4.65毫秒，而对于complex64需要22.7毫秒）。如果我将128更改为较小的值，则差异消失，然后增长（对于70最多慢20倍），最终稳定在5-6倍较慢左右。 - jorgeca

Answer 1

我还没有找到与问题有关的奇怪行为的原因，但我已经找到了一种可以显著加快计算速度的方法，这可能对其他人有用。

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

使用SPMAT.shape == (16384000, 2097152)，此操作需要2.35秒：

RES = SPMAT.dot(G)

虽然这只需 541ms：

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

结果是相同的。 我认为n.zeros的预分配可能并不必要，但我不知道还有什么其他方法可以做到这一点。

- ali_p