编辑：参见this question,我在那里学习了如何使用Numba在Python中并行化稀疏矩阵向量乘法,并且能够与Matlab结合使用.
原始问题：
我观察到,Matlab中的稀疏矩阵向量乘法比Python(使用scipy稀疏矩阵)快4到5倍.这是Matlab命令行中的一些详细信息：
>> whos A
Name Size Bytes Class Attributes
A 47166x113954 610732376 double sparse
>> whos ATrans
Name Size Bytes Class Attributes
ATrans 113954x47166 610198072 double sparse
>> nnz(A)/numel(A)
ans =
0.0071
>> whos x
Name Size Bytes Class Attributes
x 113954x1 911632 double
>> myFun = @() A*x; timeit(myFun)
ans =
0.0601
>> myFun = @() ATrans'*x; timeit(myFun)
ans =
0.0120
矩阵ATrans是A的转置.请注意,在Matlab中计算A乘以x大约需要0.06秒,但是如果我使用怪异的“转置技巧”并计算ATrans乘以x(其结果与A乘以x相同) ,计算需要0.012秒. (我不明白为什么这个技巧有效.)
以下是Python命令行的一些计时结果：
In[50]: type(A)
Out[50]:
scipy.sparse.csc.csc_matrix
In[51]: A.shape
Out[51]:
(47166, 113954)
In[52]: type(x)
Out[52]:
numpy.ndarray
In[53]: x.shape
Out[53]:
(113954L, 1L)
In[54]: timeit.timeit('A.dot(x)',setup="from __main__ import A,x",number=100)/100.0
...:
Out[54]:
0.054835451461831324
因此,对于相同的矩阵A,Python A运行时x的运行时间是Matlab运行时的4.5倍左右.如果我以csr格式存储A,那么Python运行时会更糟：
In[63]: A = sp.sparse.csr_matrix(A)
In[64]: timeit.timeit('A.dot(x)',setup="from __main__ import A,x",number=100)/100.0
...:
Out[64]:
0.0722580226496575
这是一些有关我正在使用的python版本和anaconda版本的信息：
In[2]: import sys; print('Python %s on %s' % (sys.version, sys.platform))
Python 2.7.12 |Anaconda 4.2.0 (64-bit)| (default, Jun 29 2016, 11:07:13) [MSC v.1500 64 bit (AMD64)] on win32
问题：为什么在Matlab中这种稀疏的矩阵矢量乘法比在Python中更快？以及如何在Python中使其同样快？
编辑1：这是一个线索.在Python中,如果将线程数设置为1,则密集矩阵-向量乘法的运行时间会受到严重影响,但是稀疏矩阵-向量乘法的运行时间几乎不变.
In[48]: M = np.random.rand(1000,1000)
In[49]: y = np.random.rand(1000,1)
In[50]: import mkl
In[51]: mkl.get_max_threads()
Out[51]:
20
In[52]: timeit.timeit('M.dot(y)', setup = "from __main__ import M,y", number=100) / 100.0
Out[52]:
7.232593519574948e-05
In[53]: mkl.set_num_threads(1)
In[54]: timeit.timeit('M.dot(y)', setup = "from __main__ import M,y", number=100) / 100.0
Out[54]:
0.00044465965093536396
In[56]: type(A)
Out[56]:
scipy.sparse.csc.csc_matrix
In[57]: timeit.timeit('A.dot(x)', setup = "from __main__ import A,x", number=100) / 100.0
Out[57]:
0.055780856886028685
In[58]: mkl.set_num_threads(20)
In[59]: timeit.timeit('A.dot(x)', setup = "from __main__ import A,x", number=100) / 100.0
Out[59]:
0.05550840215802509
因此,对于密集矩阵向量乘积,将线程数设置为1可使运行时间减少大约6倍.但是对于稀疏矩阵向量乘积,将线程数减少为1并不会更改运行时间.
我认为这表明在Python中,稀疏矩阵向量乘法不是并行执行的,而密集矩阵向量乘法则是在利用所有可用的核.您是否同意这个结论？如果是这样,是否有一种方法可以利用所有可用的内核在Python中进行稀疏矩阵矢量乘法？
请注意,默认情况下,Anaconda应该使用英特尔MKL优化：https://www.continuum.io/blog/developer-blog/anaconda-25-release-now-mkl-optimizations
编辑2：
我读了here,在Intel MKL中“对于稀疏矩阵,所有2级[BLAS]操作(除了稀疏三角求解器都是线程化的)”.这向我暗示scipy没有使用Intel MKL来执行稀疏矩阵矢量乘法.似乎@hpaulj(在下面发布的答案中)已经通过检查函数csr_matvec的代码来确认此结论.那么,我可以直接调用Intel MKL稀疏矩阵矢量乘法函数吗？我该怎么做？
编辑3：这是另外一个证据.当我将最大线程数设置为1时,Matlab稀疏矩阵矢量乘法运行时似乎没有变化.
>> maxNumCompThreads
ans =
20
>> myFun = @() ATrans'*x; timeit(myFun)
ans =
0.012545604076342
>> maxNumCompThreads(1) % set number of threads to 1
ans =
20
>> maxNumCompThreads % Check that max number of threads is 1
ans =
1
>> myFun = @() ATrans'*x; timeit(myFun)
ans =
0.012164191957568
这使我质疑我以前的理论,即Matlab的优势是由于多线程.
解决方法:
根据稀疏和密集的混合,我们会得到时间上的变化：
In [40]: A = sparse.random(1000,1000,.1, format='csr')
In [41]: x = np.random.random((1000,1))
In [42]: Ad = A.A
In [43]: xs = sparse.csr_matrix(x)
稀疏与密集产生密集,但稀疏与稀疏产生稀疏：
In [47]: A.dot(xs)
Out[47]:
<1000x1 sparse matrix of type '<class 'numpy.float64'>'
with 1000 stored elements in Compressed Sparse Row format>
In [48]: np.allclose(A.dot(x), Ad.dot(x))
Out[48]: True
In [49]: np.allclose(A.dot(x), A.dot(xs).A)
Out[49]: True
与替代方案相比,稀疏密集的外观看起来还不错：
In [50]: timeit A.dot(x) # sparse with dense
137 µs ± 269 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [51]: timeit Ad.dot(x) # dense with dense
1.03 ms ± 4.32 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [52]: timeit A.dot(xs) # sparse with sparse
1.44 ms ± 644 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
对于更大的矩阵,相对时间可能会有所不同.相同的不同稀疏性.
我没有访问MATLAB的权限,因此无法进行等效测试,尽管我可以在Octave上进行尝试.
这是在具有最新numpy和scipy的基本Linux(ubuntu)计算机上.
我以前的探索Directly use Intel mkl library on Scipy sparse matrix to calculate A dot A.T with less memory用矩阵计算来处理稀疏矩阵.稀疏的必须使用不同的代码.我们必须对其进行跟踪,以查看它是否将任何特殊的东西委派给了基础数学库.
csc格式的计算速度稍慢-可能是因为基本迭代是跨行的.
In [80]: Ac = A.tocsc()
In [81]: timeit Ac.dot(x)
216 µs ± 268 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
A.dot(x)可以追溯到此对编译代码的调用：
In [70]: result = np.zeros((1000,1))
In [71]: sparse._sparsetools.csr_matvec(1000,1000,A.indptr, A.indices, A.data, x, result)
_sparsetools是一个.so文件,可能是从Cython(.pyx)代码编译而成的.
在sparsetools / csr.h中,matvec的代码(模板)为：
void csr_matvec(const I n_row,
const I n_col,
const I Ap[],
const I Aj[],
const T Ax[],
const T Xx[],
T Yx[])
{
for(I i = 0; i < n_row; i++){
T sum = Yx[i];
for(I jj = Ap[i]; jj < Ap[i+1]; jj++){
sum += Ax[jj] * Xx[Aj[jj]];
}
Yx[i] = sum;
}
}
看起来像直接在csr属性(indptr,索引)上进行c迭代,相乘和求和.没有尝试利用优化的数学库或并行核.
https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h