矩陣（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​ 

我在compution時間預計幾乎直線上升，當我越來越密集矩陣的大小u乘以稀疏矩陣A（例如，第二個A*u[:,:10]假設爲220ms，最後一個A*u[:,:10] 2.2s）。但是，它比我預期的要快得多。相反，矩陣向量乘法比矩陣乘法要慢得多。有人能解釋爲什麼嗎？此外，是否有一種有效的方法來提高Matrix-vector乘法的效率與Matrix-Matrix乘法相似？

Answer 1

如果看看source code，你可以看到，csr_matvec（它實現矩陣 - 向量乘法）在C代碼的直接和循環被實施，而csr_matvecs（它實現矩陣的矩陣乘法），爲呼叫被實現到axpy BLAS程序。根據您的安裝所鏈接到的BLAS庫，此類調用可能比用於矩陣向量乘法的直接C實現更有效。這很可能就是爲什麼你看到矩陣向量乘法運算速度太慢。

更改SciPy的，以便它在矩陣矢量情況下調用BLAS可以是將包有益的貢獻。