Matlab中的稀疏矩阵向量乘法是否比Python中更快？

2019-11-11 04:59:06 阅读：285 来源： 互联网

标签：scipy sparse-matrix anaconda matlab python

编辑：参见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

标签：scipy,sparse-matrix,anaconda,matlab,python
来源： https://codeday.me/bug/20191111/2017347.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

Matlab中的稀疏矩阵向量乘法是否比Python中更快？