Matrix (scipy sparse) - 矩阵(dense; numpy array) 乘法效率

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​

当我增加密集矩阵 u 乘以稀疏矩阵 A 的大小时，我预计计算时间几乎是线性增加的（例如第二个 A*u[:,:10] 对我来说应该是 220 毫秒，最后一个 @987654325 @ 2.2 秒）。但是，它比我预期的要快得多。相反，矩阵向量乘法比矩阵矩阵乘法慢得多。有人可以解释为什么吗？此外，有没有一种有效的方法来提高矩阵向量乘法与矩阵矩阵乘法相似的效率水平？

Answer 1

如果您查看source code，您可以看到csr_matvec（实现矩阵-向量乘法）在C 代码中实现为简单的求和循环，而csr_matvecs（实现矩阵-矩阵乘法）被实现为对axpy BLAS 例程的调用。根据您的安装链接到的 BLAS 库，这样的调用可能比用于矩阵向量乘法的直接 C 实现更有效。这可能就是您看到矩阵向量乘法如此缓慢的原因。

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