python - Matlab 中的稀疏矩阵向量乘法比 Python 更快吗？

Question

编辑:参见 this question，我在其中学习了如何使用 Numba 在 Python 中并行化稀疏矩阵向量乘法，并且能够与 Matlab 配合使用。

---

原始问题:

我观察到稀疏矩阵向量乘法在 Matlab 中比在 Python 中快 4 或 5 倍(使用 scipy 稀疏矩阵)。以下是来自 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，执行 A 乘以 x 的 Python 运行时间大约是 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 应该默认使用 Intel MKL 优化:https://www.continuum.io/blog/developer-blog/anaconda-25-release-now-mkl-optimizations

---

编辑 2: 我读到 here，在英特尔 MKL 中“对于稀疏矩阵，除稀疏三角求解器之外的所有 2 级 [BLAS] 操作都是线程化的”。这向我表明 scipy 没有使用 Intel MKL 来执行稀疏矩阵向量乘法。似乎@hpaulj(在下面发布的答案中)通过检查函数 csr_matvec 的代码证实了这个结论。那么，我可以直接调用英特尔 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 的优势在于多线程。

Answer 1

根据稀疏和密集的混合，我们得到时间的变化:

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)

sparse 与 dense 产生密集，而 sparse 与 sparse 产生稀疏:

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 ，用矩阵计算处理稀疏矩阵。稀疏与密集必须使用不同的代码。我们必须跟踪它以查看它是否将任何特殊内容委托(delegate)给底层数学库。

---

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、indices)的直接 C++ 迭代、乘法和求和。没有尝试使用优化的数学库或并行内核。

https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h