Matlab中的多线程稀疏矩阵乘法

Question

我正在执行几个NxM稀疏（约1-2％）矩阵（称为B）和一个NxM密集矩阵（称为A，其中M < N）的矩阵乘法。N和M都很大，大约有数千个。我正在运行Matlab 2013a。

通常，在Matlab中，矩阵乘法和大多数其他矩阵操作隐式地并行化，即它们自动使用多个线程。如果两个矩阵中有任何一个是稀疏的，则似乎不是这种情况（请参见例如this StackOverflow discussion - 没有针对预期问题的答案 - 和 this largely unanswered MathWorks thread）。这对我来说是一个非常不愉快的惊喜。

我们可以通过以下代码验证多线程对稀疏矩阵操作没有影响：

clc; clear all; 

N = 5000;         % set matrix sizes
M = 3000;       
A = randn(N,M);   % create dense random matrices
B = sprand(N,N,0.015); % create sparse random matrix
Bf = full(B);     %create a dense form of the otherwise sparse matrix B

for i=1:3 % test for 1, 2, and 4 threads
  m(i) = 2^(i-1);
  maxNumCompThreads(m(i)); % set the thread count available to Matlab
  tic                      % starts timer
    y = B*A; 
  walltime(i) = toc;       % wall clock time
  speedup(i) = walltime(1)/walltime(i);
end

% display number of threads vs. speed up relative to just a single thread
[m',speedup']

这将产生以下输出，说明在稀疏操作中使用1、2和4个线程没有区别：

threads   speedup
1.0000    1.0000
2.0000    0.9950
4.0000    1.0155

如果我将B替换为其密集形式（在上面称为Bf），则会显著提高速度：

threads   speedup
1.0000    1.0000
2.0000    1.8894
4.0000    3.4841

（说明Matlab中稠密矩阵的矩阵运算确实是隐式并行化的）

所以，我的问题是：是否有任何方法可以在不将稀疏矩阵转换为密集形式的情况下访问并行化/线程化版本的矩阵操作（在Matlab中进行稀疏矩阵操作）？ 我发现一个旧的建议涉及MathWorks的.mex文件, 但似乎链接已经失效，而且文档不太完整/没有反馈？有其他选择吗？

这似乎是隐式并行功能的一个相当严重的限制，因为稀疏矩阵在计算量大的问题中很常见，并且这些情况下高度渴望超线程功能。

- Thomas Christensen

---

@Yvon 在这些链接中，我看到了对事物运作方式的一般描述，然而我无法确定它们与问题的相关性。 - Dennis Jaheruddin

只是一个愚蠢的想法：将完整矩阵变为稀疏矩阵是否有帮助？ - Dennis Jaheruddin

@DennisJaheruddin 这有助于加快速度，但在内存方面并不实用。这就是提出这个问题的原因。 - Daniyar

@Daniyar 由于格式问题，M 很大的信息被隐藏了，我已经编辑了问题以解决这个问题。-- 不过，从完整到稀疏只会使矩阵存储变为两倍大，所以除非你接近内存限制，否则这可能是一个有趣的方法。 - Dennis Jaheruddin

@DennisJaheruddin 哦。不过问题一样。 - Daniyar

Answer 1

MATLAB已经使用SuiteSparse by Tim Davis进行许多处理稀疏矩阵的操作（例如此处），但我认为这两个都不是多线程的。

通常，稀疏矩阵上的计算受内存限制而不是CPU限制。因此，即使您使用多线程库，我怀疑您不会在性能方面看到巨大的好处，至少不可与专门用于密集矩阵的库相比...

毕竟，稀疏矩阵的设计目标与常规的密集矩阵有所不同，其中高效的内存存储通常更加重要。

---

我快速 在网上搜索，发现有几个实现：

• 稀疏BLAS, spBLAS, PSBLAS。例如，Intel MKL和AMD ACML对稀疏矩阵有一定的支持。
• cuSPARSE, CUSP, VexCL, ViennaCL等可在GPU上运行。

- Amro

---

支持稀疏矩阵的库不一定支持稀疏矩阵的多线程乘法。 - Daniyar

@Daniyar：也许不是参考实现，但像Intel MKL中那样优化的程序可以同时处理密集和稀疏例程的多线程：https://software.intel.com/en-us/articles/parallelism-in-the-intel-math-kernel-library。还要注意GPU库本质上是并行的。 - Amro

@Daniyar：顺便说一下，MATLAB确实支持在分布式内存场景下的并行稀疏矩阵（与共享内存多线程相对应），使用Parallel Computing Toolbox http://www.mathworks.com/help/distcomp/sparse.html。这是基于[ScaLAPACK](http://www.netlib.org/scalapack/index.html)库实现的。 - Amro

如果在单台计算机上使用并行计算，那么它不会允许多线程吗？（正如我在我的答案中建议的那样） - Dennis Jaheruddin

@DennisJaheruddin：你在谈论哪种实现？如果是Intel MKL，那么是的。但MATLAB不使用MKL的稀疏例程，而是使用SuiteSparse库。至于PCT工具箱中的“分布式”数组，则使用多进程而不是多线程并行化（类似于MPI）。 - Amro

Answer 2

我最终使用OpenMP编写了自己的mex文件来实现多线程。代码如下。编译时不要忘记使用-largeArrayDims和/openmp（或-fopenmp）标志。

#include <omp.h>
#include "mex.h"
#include "matrix.h"

#define ll long long

void omp_smm(double* A, double*B, double* C, ll m, ll p, ll n, ll* irs, ll* jcs)
{
    for (ll j=0; j<p; ++j)
    {
        ll istart = jcs[j];
        ll iend = jcs[j+1];
        #pragma omp parallel for
        for (ll ii=istart; ii<iend; ++ii)
        {
            ll i = irs[ii];
            double aa = A[ii];
            for (ll k=0; k<n; ++k)
            {
                C[i+k*m] += B[j+k*p]*aa;
            }
        }
    }
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *A, *B, *C; /* pointers to input & output matrices*/
    size_t m,n,p;      /* matrix dimensions */

    A = mxGetPr(prhs[0]); /* first sparse matrix */
    B = mxGetPr(prhs[1]); /* second full matrix */

    mwIndex * irs = mxGetIr(prhs[0]);
    mwIndex * jcs = mxGetJc(prhs[0]);

    m = mxGetM(prhs[0]);  
    p = mxGetN(prhs[0]);
    n = mxGetN(prhs[1]);

    /* create output matrix C */
    plhs[0] = mxCreateDoubleMatrix(m, n, mxREAL);
    C = mxGetPr(plhs[0]);

    omp_smm(A,B,C, m, p, n, (ll*)irs, (ll*)jcs);
}

- Daniyar

---

即使这个朴素算法的运行时间为立方级 O(m*p*n)，如果提供代码，也会得到加分。有趣的是，将其与Intel MKL中优化的（多线程）实现进行比较，即mkl_dcsrmm（使用CSR格式的一般稀疏矩阵A和密集矩阵B和C进行矩阵乘积C=A*B的例程）。 - Amro

@Amro 这种方法有些天真，因为它使用了Matlab的数据结构。非天真版本可能会表现更好，但这取决于矩阵的稀疏程度。 - Daniyar

@Daniyar，您能否包含一个基准测试？最好是显示出在典型情况下这个方案更快的一个，以及显示出基本功能更快的一个。 - Dennis Jaheruddin

@DennisJaheruddin 这应该至少与默认的Matlab乘法一样快，因为它采用了多线程技术。 - Daniyar

@Daniyar 或许它“应该是”，但看到它真的是什么样子（以及速度有多快）会很好。 - Dennis Jaheruddin

Answer 3

在Matlab中心上，有人问了同样的问题，并给出了这个答案：

I believe the sparse matrix code is implemented by a few specialized TMW engineers rather than an external library like BLAS/LAPACK/LINPACK/etc... 

这基本上意味着，你没有运气了。

---

然而，我可以想到一些技巧来实现更快的计算：

1. 如果你需要进行多个乘法操作：同时进行多个乘法操作并以并行方式处理？
2. 如果你只想进行一个乘法操作：将矩阵分成几份（例如上半部分和下半部分），并行计算各部分，并在之后合并结果。

可能这些解决方案不会像正确实现的多线程那样快，但希望你仍然能够获得加速。

- Dennis Jaheruddin