python - Matrix (scipy sparse) - Matrix (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 毫秒，而最后一个 A*u[:,:10] 是 2.2 秒)。但是，它比我预期的要快得多。相反，矩阵-向量乘法比矩阵-矩阵乘法慢得多。有人可以解释为什么吗？此外，是否有一种有效的方法可以将矩阵-向量乘法的效率提高到与矩阵-矩阵乘法相似的水平？

Answer 1

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

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