为什么MATLAB在矩阵乘法中如此之快？

14 回答

• 1

  在进行矩阵乘法时，使用naive乘法方法，该方法需要 O(n^3) 的时间 .

  存在矩阵乘法算法，需要 O(n^2.4) . 这意味着在 n=2000 ，您的算法需要的计算量是最佳算法的100倍 .
  您应该检查维基百科页面中的矩阵乘法，以获取有关实现它的有效方法的更多信息 .

  回复于 2024-05-01T22:36:34+08:00
• 38

  你需要小心与C进行公平的比较 . 你可以发布显示你用于矩阵乘法的核心内部循环的C代码吗？大多数情况下，我关心你的记忆布局，以及你是否在做浪费的事情 .

  我编写的C矩阵乘法与Matlab一样快，但需要注意 . （编辑：在Matlab使用GPU之前 . ）

  几乎可以肯定Matlab在这些“内置”功能上浪费了很少的周期 . 我的问题是，你在哪里浪费周期？ （没有冒犯的意思）

  回复于 2024-05-01T22:36:34+08:00
• 4

  这是我在使用Tesla C2070的机器上使用MATLAB R2011a Parallel Computing Toolbox的结果：

  >> A = rand(1024); gA = gpuArray(A);
  % warm up by executing the operations a couple of times, and then:
  >> tic, C = A * A; toc
  Elapsed time is 0.075396 seconds.
  >> tic, gC = gA * gA; toc
  Elapsed time is 0.008621 seconds.
  

  MATLAB使用高度优化的库进行矩阵乘法，这就是普通MATLAB矩阵乘法如此之快的原因 . gpuArray 版本使用MAGMA .

  Update using R2014a 在具有Tesla K20c的机器上，以及新的 timeit 和 gputimeit 功能：

  >> A = rand(1024); gA = gpuArray(A);
  >> timeit(@()A*A)
  ans =
      0.0324
  >> gputimeit(@()gA*gA)
  ans =
      0.0022
  

  Update using R2018b 在具有16个物理内核和一个Tesla V100的WIN64机器上：

  >> timeit(@()A*A)
  ans =
      0.0229
  >> gputimeit(@()gA*gA)
  ans =
     4.8019e-04
  

  回复于 2024-05-01T22:36:34+08:00
• 148

  这种问题反复出现，应该比Stackoverflow上的“Matlab使用高度优化的库”或“Matlab使用MKL”更清楚地回答 .

  History:

  矩阵乘法（与矩阵向量，向量 - 向量乘法和许多矩阵分解一起）是线性algrebra中最重要的问题 . 从早期开始，工程师就一直在用计算机解决这些问题 .

  我不是历史专家，但显然那时候，每个人都只用简单的循环重写了他的Fortran版本 . 然后出现了一些标准化，识别出需要解决的大多数线性代数问题的“内核”（基本例程） . 然后，这些基本操作在称为“基本线性代数子程序（BLAS）”的规范中标准化 . 然后，工程师可以在他们的代码中调用这些经过良好测试的标准BLAS例程，使他们的工作变得更加容易 .

  BLAS:

  BLAS从1级（定义标量矢量和矢量矢量运算的第一个版本）演变为2级（矢量矩阵运算）到3级（矩阵矩阵运算），并提供越来越多的"kernels"标准化越来越多基本线性代数运算 . 最初的Fortran 77实现仍然可以在Netlib's website上使用 .

  Towards better performance:

  因此，多年来（特别是在BLAS 1级和2级版本之间：80年代早期），随着向量操作和缓存层次结构的出现，硬件发生了变化 . 这些演进使得有可能大大提高BLAS子程序的性能 . 然后不同的供应商出现了BLAS例程的实现，这些例程越来越高效 .

  我不知道所有的历史实现（当时我还没出生或是个孩子），但是最着名的两个是在21世纪初出现的：英特尔MKL和GotoBLAS . 您的Matlab使用的是英特尔MKL，这是一款非常优秀的优化BLAS，它可以解释您所看到的卓越性能 .

  Technical details on Matrix multiplication:

  那么为什么Matlab（MKL）在 dgemm （双精度通用矩阵 - 矩阵乘法）如此快？简单来说：因为它使用矢量化和良好的数据缓存 . 更复杂的术语：请参阅Jonathan Moore提供的article .

  基本上，当您在提供的C代码中执行乘法运算时，您根本不熟悉缓存 . 由于我怀疑你创建了一个指向行数组的指针数组，因此你在内部循环中对"matice2"： matice2[m][k] 的第k列的访问非常慢 . 实际上，当您访问 matice2[0][k] 时，您必须获得矩阵的数组0的第k个元素 . 然后在下一次迭代中，您必须访问 matice2[1][k] ，这是另一个数组（数组1）的第k个元素 . 然后在下一次迭代中，您访问另一个数组，依此类推......由于整个矩阵 matice2 可以't fit in the highest caches (it' s 8*1024*1024 字节大），程序必须从主内存中获取所需的元素，从而浪费大量时间 .

  如果你只是转换了矩阵，那么访问将在连续的内存地址中，你的代码已经运行得更快，因为现在编译器可以同时加载缓存中的整行 . 试试这个修改过的版本：

  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();
  

  因此，您可以看到缓存局部性如何大大提高代码的性能 . 现在真正的 dgemm 实现将它用于非常广泛的层次：它们对由TLB的大小定义的矩阵块进行乘法（转换后备缓冲区，长话短说：可以有效缓存的内容），以便它们流式传输到处理器确切地处理它可以处理的数据量 . 另一个方面是矢量化，他们使用处理器's vectorized instructions for optimal instruction throughput, which you can' t真正从您的跨平台C代码 .

  最后，人们声称这是因为Strassen或Coppersmith-Winograd算法错误的是，由于上面提到的硬件考虑因素，这两种算法在实践中都不可实现 .

  回复于 2024-05-01T22:36:34+08:00
• 8

  This is why . MATLAB不会像在C代码中那样循环遍历每个元素，从而不执行简单的矩阵乘法 .

  当然我假设您只是使用 C=A*B 而不是自己编写乘法函数 .

  回复于 2024-05-01T22:36:34+08:00
• 1

  Matlab在不久前收录了LAPACK，所以我假设他们的矩阵乘法使用至少那么快的东西 . LAPACK源代码和文档随时可用 .

  你也可以看看Goto和Van De Geijn的论文"Anatomy of High-Performance Matrix Multiplication" at http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

  回复于 2024-05-01T22:36:34+08:00
• 18

  答案是LAPACK和BLAS库使MATLAB在矩阵运算中的速度非常快，而不是MATLAB人员的任何专有代码 .

  使用C代码中的LAPACK和/或BLAS库进行矩阵运算，您应该获得与MATLAB类似的性能 . 这些图书馆应该可以在任何现代系统上免费获得，而且学术界已经开发了几十年的部件 . 请注意，有多个实现，包括一些封闭源，如Intel MKL .

  关于BLAS如何获得高性能的讨论is available here.

  ---

  顺便说一句，直接从c调用LAPACK库是一种严重的痛苦（但值得） . 您需要非常准确地阅读文档 .

  回复于 2024-05-01T22:36:34+08:00
• 1

  根据您的Matlab版本，我相信它可能已经在使用您的GPU了 .

  另一件事; Matlab会跟踪矩阵的许多属性;它的对角线，hermetian等等，并专门研究其基于此的算法 . 也许它的专业化基于您传递的零矩阵，或类似的东西？也许它正在缓存重复的函数调用，这会弄乱你的时间？也许它优化了重复使用的矩阵产品？

  为了防止发生这种情况，请使用随机数字矩阵，并确保通过将结果打印到屏幕或磁盘或其他某些部分来强制执行 .

  回复于 2024-05-01T22:36:34+08:00
• 7

  使用双精度数和一个实数数组而不是三个单独的数组导致我的C＃代码与C / Java几乎相同（使用您的代码：1024 - 更快一点，2048 - 大约140s和4096 - 大约22分钟）

  1024x1024   2048x2048   4096x4096
                  ---------   ---------   ---------
  your C++ (ms)   6137.10     64369.29     551390.93
  my C# (ms)      9730.00     90875.00    1062156.00
  

  这是我的代码：

  const int rozmer = 1024;
      double[][] matice1 = new double[rozmer * 3][];
      Random rnd = new Random();
  
      public Form1()
      {
          InitializeComponent();
  
          System.Threading.Thread thr = new System.Threading.Thread(new System.Threading.ThreadStart(() =>
          {
  
              string res = "";
              Stopwatch timer = new Stopwatch();
              timer.Start();
  
              double temp = 0;
              int r2 = rozmer * 2;
  
              for (int i = 0; i < rozmer*3; i++)
              {
                  if (matice1[i] == null)
                  {
                      matice1[i] = new double[rozmer];
                      {
                          for (int e = 0; e < rozmer; e++)
                          {
                              matice1[i][e] = rnd.NextDouble();
                          }
                      }
                  }
              }
              timer.Stop();
              res += timer.ElapsedMilliseconds.ToString();
  
              int j = 0; int k = 0; int m = 0;
  
              timer.Reset();
              timer.Start();
              for (j = 0; j < rozmer; j++)
              {
                  for (k = 0; k < rozmer; k++)
                  {
                      temp = 0;
                      for (m = 0; m < rozmer; m++)
                      {
                          temp = temp + matice1[j][m] * matice1[m + rozmer][k];
                      }
                      matice1[j + r2][k] = temp;
                  }
              }
              timer.Stop();
              this.Invoke((Action)delegate
              {
                  this.Text = res + " : " + timer.ElapsedMilliseconds.ToString();
              });
          }));
          thr.Start();
      }
  

  回复于 2024-05-01T22:36:34+08:00
• 5

  您是否检查过所有实现都使用了算法的多线程优化？他们使用相同的乘法算法吗？

  我真的很怀疑 .

  Matlab本身并不快，你可能使用了慢速实现 .

  Algorithms for efficient matrix multiplication

  回复于 2024-05-01T22:36:34+08:00
• 7

  “为什么matlab在执行xxx时比其他程序更快”的一般答案是matlab有很多内置的优化函数 .

  使用的其他程序通常没有这些功能，所以人们应用自己的创意解决方案，这比专业优化的代码慢得多 .

  这可以通过两种方式解释：

  1）常见/理论方式：Matlab并没有明显更快，你只是做错了基准测试

  2）现实的方法：对于这个东西，Matlab在实践中更快，因为作为c的语言太容易以无效的方式使用 .

  回复于 2024-05-01T22:36:34+08:00
• 74

  鲜明的对比不仅仅是由于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);
  };
  

  和

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

  应该使用相同的乘法算法，以使翻牌的数量相同 . （对于大小为n的平方矩阵，n ^ 3）

  我要求你计时，以便结果与你之前的（在同一台机器上）相当 . 随着比较一下，您将准确显示内存访问时间的重要性！

  回复于 2024-05-01T22:36:34+08:00
• 1

  MATLAB使用来自英特尔的高度优化的LAPACK实现，称为Intel Math Kernel Library（英特尔MKL） - 特别是dgemm function . 速度该库利用处理器功能，包括SIMD指令和多核处理器 . 他们没有记录他们使用的具体算法 . 如果您从C调用英特尔MKL，您应该会看到类似的性能 .

  我不确定MATLAB用于GPU乘法的库，但可能类似nVidia CUBLAS .

  回复于 2024-05-01T22:36:34+08:00
• 1

  它在C中很慢，因为你没有使用多线程 . 基本上，如果A = BC，它们都是矩阵，A的第一行可以独立于第二行计算，等等 . 如果A，B和C都是n×n矩阵，你可以加速乘法系数n ^ 2，如

  a_ {i，j} = sum_ b_ {i，k} c_ {k，j}

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

  回复于 2024-05-01T22:36:34+08:00