矩阵（Scipy稀疏矩阵）和矩阵（密集的；Numpy数组）相乘的效率

Question

我是一名从事地球物理反演研究的研究员。这项工作需要解决线性系统：Au = rhs。其中A通常是稀疏矩阵，但rhs和u可以是密集矩阵或向量。为了进行基于梯度的反演，我们需要计算灵敏度，并且这需要进行许多矩阵-矩阵和矩阵-向量乘法。最近我发现了一个矩阵（稀疏）-矩阵（密集）乘法中的奇怪行为，以下是一个例子：

import numpy as np
import scipy.sparse as sp
n = int(1e6)
m = int(100)
e = np.ones(n)
A = sp.spdiags(np.vstack((e, e, e)), np.array([-1, 0, 1]), n, n)
A = A.tocsr()
u = np.random.randn(n,m)

%timeit rhs = A*u[:,0]
#10 loops, best of 3: 22 ms per loop    
%timeit rhs = A*u[:,:10]
#10 loops, best of 3: 98.4 ms per loop
%timeit rhs = A*u
#1 loop, best of 3: 570 ms per loop​

我原本期望当我增加密矩阵乘以稀疏矩阵A的大小时，计算时间会呈几乎线性增长（例如第二个应该是220毫秒，最后一个应该是2.2秒）。然而，它比我预期的要快得多。相反，矩阵向量乘法比矩阵乘法慢得多。有人能解释为什么吗？此外，是否有一种有效的方法来提高矩阵向量乘法的效率，使其与矩阵乘法类似？

- SEOGI KANG

---

你/我需要深入挖掘函数调用栈来处理这些乘法。这并不是一项简单的任务。显然，这不仅仅是迭代u的列并收集值。这是一个“稀疏密集=>密集”的情况，与“稀疏稀疏”或“密集*密集”不同。 - hpaulj

Answer 1

如果你查看源代码，你会发现实现矩阵向量乘法的csr_matvec是在C语言中实现的简单求和循环，而实现矩阵矩阵乘法的csr_matvecs则调用了axpy BLAS例程。根据你的安装链接到的BLAS库不同，这样的调用可能比用于矩阵向量乘法的简单C实现要高效得多。这可能就是你看到矩阵向量乘法如此缓慢的原因。

更改scipy以便在矩阵向量情况下调用BLAS可能对该软件包是一个有用的贡献。

- jakevdp

---

我明白了，那很有道理！谢谢jaekvdp。我想改变那个，但是我没有使用过C语言，所以我不确定我能在这方面做出贡献。嗯...但仍然认识到我们可以改变的地方是一个好的开始！ - SEOGI KANG