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

【问题标题】：Multithreaded sparse matrix multiplication in MatlabMatlab中的多线程稀疏矩阵乘法
【发布时间】：2014-10-09 09:07:23
【问题描述】：

我正在执行 NxN 稀疏 (~1-2%) 矩阵的几次矩阵乘法，我们称之为 B，使用 NxM 密集矩阵，我们称之为 A（其中 M

现在，矩阵乘法和大多数其他矩阵运算通常在 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 中）的并行化/线程化版本的矩阵运算而不将它们转换为密集形式？我找到了一个旧的suggestion involving .mex files at MathWorks，但似乎链接已经失效并且没有很好的记录/没有反馈？有其他选择吗？

这似乎是对隐式并行功能的相当严格的限制，因为稀疏矩阵存在大量计算量大的问题，并且在这些情况下非常需要超线程功能。

【问题讨论】：

mathworks.com/help/matlab/math/sparse-matrix-operations.html 和 mathworks.com/matlabcentral/answers/…
@Yvon 在链接中，我看到了关于事情如何运作的一般描述，但我无法确定它与问题的相关性。
只是一个愚蠢的事后思考：它有助于使整个矩阵稀疏吗？
@DennisJaheruddin 它在加速方面有所帮助，但在内存方面并不实用。这就是问题的原因。
@Daniyar M 大的信息被格式问题隐藏，已编辑问题以解决此问题。 -- 尽管如此，从完整到稀疏应该只使矩阵存储两倍大，所以除非你接近内存限制，否则它可能是一个有趣的方法。

标签： multithreading matlab matrix-multiplication sparse-matrix

【解决方案1】：

MATLAB 已经将 Tim Davis 的 SuiteSparse 用于其对稀疏矩阵的许多操作（例如 see here），但我认为它们都不是多线程的。

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

毕竟design of sparse matrices 的目标与常规密集矩阵不同，后者通常更重要的是高效的内存存储。

我做了一个快速的search online，并在那里找到了一些实现：

sparse BLAS、spBLAS、PSBLAS。例如，Intel MKL 和 AMD ACML 确实支持稀疏矩阵
在 GPU 上运行的cuSPARSE、CUSP、VexCL、ViennaCL 等。

【讨论】：

支持稀疏矩阵的库不一定支持稀疏矩阵的多线程乘法。
@Daniyar：可能不是参考实现，但英特尔 MKL 中的优化实现对于密集和稀疏例程都是多线程的：software.intel.com/en-us/articles/…。请注意提及本质上是并行的 GPU 库。
@Daniyar：顺便说一下，MATLAB 使用并行计算工具箱mathworks.com/help/distcomp/sparse.html 支持distributed-memory 场景（与共享内存多线程相反）中的并行稀疏矩阵。这是基于ScaLAPACK 库。
如果在单台计算机上使用并行计算，会不会允许多线程？（正如我在回答中建议的那样）
@DennisJaheruddin：你说的是哪个实现？英特尔 MKL，然后是的。但是 MATLAB 不使用 MKL 中的稀疏例程，而是使用 SuiteSparse 库。至于 PCT 工具箱distributed 数组，它们与多进程而不是多线程并行化（想想MPI）

【解决方案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);
}

【讨论】：

+ 1 用于提供代码，即使这个naive algorithm 效率低下，运行时间为立方O(m*p*n)。看看这与英特尔 MKL 中的优化（和多线程）实现相比如何会很有趣，即例程：mkl_dcsrmm（矩阵矩阵乘积C=A*B 和A CSR 格式的通用稀疏矩阵和B和C密集矩阵）
@Amro 这太天真了，因为它使用 matlab 的数据结构。非朴素版本可能表现更好，但这取决于矩阵的稀疏性。
@Daniyar 你能提供一个基准吗？最好是一个显示速度更快的典型案例，以及一个基本功能更快的案例？
@DennisJaheruddin 这应该至少和默认的 matlab 乘法一样快，因为它是多线程的。
@Daniyar 也许它“应该是”，但很高兴看到它确实是（以及快了多少）。

【解决方案3】：

在matlab central 上提出了同样的问题，并给出了这个答案：

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

这基本上意味着你运气不好。

但是我可以想到一些技巧来实现更快的计算：

如果您需要进行多次乘法运算：一次进行多次乘法运算并并行处理？
如果您只想做一次乘法：将矩阵切成小块（例如上半部分和下半部分），并行计算各部分，然后合并结果

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

【讨论】：