我用CUDA, c++, c#, Java做了一些基准测试,并使用MATLAB进行验证和矩阵生成。当我用MATLAB执行矩阵乘法时,2048x2048甚至更大的矩阵几乎立即被相乘。

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

只有CUDA是有竞争力的,但我认为至少c++会有点接近,而不是慢60倍。我也不知道如何看待c#的结果。算法与c++和Java一样,但从1024年到2048年有了巨大的飞跃。

MATLAB是如何如此快速地执行矩阵乘法的?

c++代码:

float temp = 0;
timer.start();
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

当前回答

它在c++中很慢,因为你没有使用多线程。本质上,如果A = B C,其中它们都是矩阵,则A的第一行可以独立于第二行计算,等等。如果A、B和C都是n × n矩阵,您可以将乘法运算速度提高一个因子n^2,如

A_ {i,j} = sum_{k} b_{i,k} c_{k,j}

如果您使用Eigen [http://eigen.tuxfamily.org/dox/GettingStarted.html],多线程是内置的,线程的数量是可调的。

其他回答

这种问题反复出现,在Stack Overflow上应该比“MATLAB使用高度优化的库”或“MATLAB使用MKL”更清楚地回答。

历史:

矩阵乘法(连同矩阵-向量、向量-向量乘法和许多矩阵分解)是线性代数中最重要的问题。工程师们从早期开始就一直在用计算机解决这些问题。

I'm not an expert on the history, but apparently back then, everybody just rewrote his FORTRAN version with simple loops. Some standardization then came along, with the identification of "kernels" (basic routines) that most linear algebra problems needed in order to be solved. These basic operations were then standardized in a specification called: Basic Linear Algebra Subprograms (BLAS). Engineers could then call these standard, well-tested BLAS routines in their code, making their work much easier.

布拉斯特区:

BLAS从第1级(定义标量-向量和向量-向量运算的第一个版本)发展到第2级(向量-矩阵运算),再到第3级(矩阵-矩阵运算),并提供了越来越多的“核心”,从而标准化了越来越多的基本线性代数运算。最初的FORTRAN 77实现仍然可以在Netlib的网站上找到。

为了更好的表现:

因此,多年来(特别是在BLAS级别1和级别2发布之间:80年代初),随着矢量操作和缓存层次结构的出现,硬件发生了变化。这些演进使得大幅度提高BLAS子例程的性能成为可能。然后,不同的供应商也推出了他们越来越高效的BLAS例程实现。

我不知道所有历史上的实现(那时候我还没出生,也不是孩子),但最著名的两个实现出现在21世纪初:英特尔MKL和GotoBLAS。你的Matlab使用的是英特尔MKL,这是一个非常好的,优化的BLAS,这就解释了你所看到的出色性能。

矩阵乘法的技术细节:

那么为什么Matlab (MKL)在dgemm(双精度一般矩阵-矩阵乘法)上如此之快?简单来说:因为它使用了向量化和良好的数据缓存。用更复杂的术语来说:请参阅乔纳森•摩尔(Jonathan Moore)提供的文章。

Basically, when you perform your multiplication in the C++ code you provided, you are not at all cache-friendly. Since I suspect you created an array of pointers to row arrays, your accesses in your inner loop to the k-th column of "matice2": matice2[m][k] are very slow. Indeed, when you access matice2[0][k], you must get the k-th element of the array 0 of your matrix. Then in the next iteration, you must access matice2[1][k], which is the k-th element of another array (the array 1). Then in the next iteration you access yet another array, and so on... Since the entire matrix matice2 can't fit in the highest caches (it's 8*1024*1024 bytes large), the program must fetch the desired element from main memory, losing a lot of time.

如果您只是调换了矩阵的位置,以便访问相邻的内存地址,那么您的代码将运行得更快,因为现在编译器可以同时在缓存中加载整个行。试试这个修改过的版本:

timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
    for (int q = 0; q < rozmer; q++)
    {
        tempmat[p][q] = matice2[q][p];
    }
}
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * tempmat[k][m];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

So you can see how just cache locality increased your code's performance quite substantially. Now real dgemm implementations exploit that to a very extensive level: They perform the multiplication on blocks of the matrix defined by the size of the TLB (Translation lookaside buffer, long story short: what can effectively be cached), so that they stream to the processor exactly the amount of data it can process. The other aspect is vectorization, they use the processor's vectorized instructions for optimal instruction throughput, which you can't really do from your cross-platform C++ code.

最后,有人声称这是因为Strassen's或Coppersmith-Winograd算法是错误的,这两个算法在实践中都是不可实现的,因为上面提到的硬件方面的考虑。

答案是LAPACK和BLAS库使MATLAB在矩阵运算方面速度惊人,而不是MATLAB的任何专有代码。

在你的c++代码中使用LAPACK和/或BLAS库来进行矩阵运算,你会得到与MATLAB相似的性能。这些库在任何现代系统上都应该是免费的,其中一部分是学术界在几十年里开发出来的。注意,有多种实现,包括一些封闭源代码,如Intel MKL。

这里有关于BLAS如何获得高性能的讨论。


顺便说一句,在我的经验中,直接从c调用LAPACK库是一种严重的痛苦(但值得)。你需要非常精确地阅读文档。

这种鲜明的对比不仅是由于Matlab的惊人优化(正如许多其他答案已经讨论过的那样),而且是由于您将矩阵作为一个对象来表述的方式。

看起来你把矩阵变成了列表的列表?列表的列表包含指向列表的指针,然后包含您的矩阵元素。所包含列表的位置是任意分配的。当循环遍历第一个索引(行号?)时,内存访问时间非常重要。相比之下,为什么不尝试实现矩阵作为一个单一的列表/向量使用下面的方法?

#include <vector>

struct matrix {
    matrix(int x, int y) : n_row(x), n_col(y), M(x * y) {}
    int n_row;
    int n_col;
    std::vector<double> M;
    double &operator()(int i, int j);
};

And

double &matrix::operator()(int i, int j) {
    return M[n_col * i + j];
}

应该使用相同的乘法算法,以使触发器的数量相同。(大小为n的方阵为n^3)

我要求您对它进行计时,以便结果与前面(在同一台机器上)的结果相比较。通过比较,您将准确地显示内存访问时间有多么重要!

这就是原因。MATLAB不像在c++代码中那样,通过遍历每一个元素来执行简单的矩阵乘法。

当然,我假设你只是用C=A*B而不是自己写一个乘法函数。

当做矩阵乘法时,你使用朴素乘法,它需要O(n^3)的时间。

有一个矩阵乘法算法,它需要O(n^2.4)。这意味着当n=2000时,你的算法需要的计算量是最佳算法的100倍。 你真的应该去维基百科上查看矩阵乘法的页面,以获得关于有效实现矩阵乘法的进一步信息。