【问题标题】：Laderman's 3x3 matrix multiplication with only 23 multiplications, is it worth it?Laderman 的 3x3 矩阵乘法只有 23 次乘法，值得吗？
【发布时间】：2012-06-05 07:59:34
【问题描述】：

取两个 3x3 矩阵 A*B=C 的乘积。天真地，这需要使用 standard algorithm 进行 27 次乘法运算。如果有人聪明，你可以只使用 23 次乘法，a result found in 1973 by Laderman。该技术涉及节省中间步骤并以正确的方式组合它们。

现在让我们修复一种语言和一种类型，比如带有double 元素的 C++。如果 Laderman 算法是硬编码而不是简单的双循环，我们能否期望现代编译器的性能能够消除算法的差异？

关于这个问题的说明：这是一个编程网站，这个问题是在时间关键内循环的最佳实践的上下文中提出的；过早优化这不是。作为 cmets，我们非常欢迎有关实施的提示。

【问题讨论】：

"...我们能否期望现代编译器的性能能够消除算法的差异？"为什么不试试呢？编写两个代码，分别运行 1000 次并比较运行时间。
该问题的一般答案是“否”。世界上仍然需要聪明的算法。
@phs：回答标题中的问题，还是回答注释上方的问题？他们是相反的。
由于 Makarov 还有一种技术只需要 22 次乘法 (original in Russian, PDF;English translation, behind paywall)
后期评论：拉德曼算法的优点只有在递归地用于更大的矩阵时才会显现出来。 Paolo D'Alberto 在他的FastMMW site 中对此进行了实验分析。在乘法复杂度方面，Laderman 算法 O(n^2.854) 优于朴素 O(n^3)，但与 Strassen 算法 O(n^2.807) 相比逊色。具有 21 次乘法 O(n^2.77) 的 3x3 算法将优于 Strassen 的算法。

标签： c++ algorithm linear-algebra matrix-multiplication

【解决方案1】：

关键是掌握平台上的指令集。这取决于您的平台。有几种技术，当您倾向于需要最大可能的性能时，您的编译器将附带分析工具，其中一些内置了优化提示。对于最细粒度的操作，请查看汇编器输出，看看是否有任何改进也在那个水平。

同时指令多个数据命令对多个操作数并行执行相同的操作。这样你就可以拿了

double a,b,c,d;
double w = d + a; 
double x = a + b;
double y = b + c;
double z = c + d;

替换成

double256 dabc = pack256(d, a, b, c);
double256 abcd = pack256(a, b, c, d);
double256 wxyz = dabc + abcd;

因此，当值被加载到寄存器中时，它们会被加载到一个 256 位宽的寄存器中，用于某些具有 256 位宽寄存器的虚构平台。

浮点是一个重要的考虑因素，一些 DSP 可以显着加快对整数的运算速度。 GPU 在浮点方面往往表现出色，尽管有些 GPU 在单精度上的速度要快 2 倍。这个问题的 3x3 情况可以适合单个 CUDA 线程，因此您可以同时流式传输 256 个这些计算。

选择您的平台，阅读文档，实施几种不同的方法并对其进行分析。

【讨论】：

【解决方案2】：

计时测试：

我自己进行了计时测试，结果让我感到惊讶（因此我首先问了这个问题）。简而言之，在标准编译下，laderman 的速度要快 225%，但使用 -03 优化标志时，它会慢 50%！在-O3 标志或编译器完全优化简单乘法期间，我每次都必须在矩阵中添加一个随机元素，在时钟精度内花费零时间。由于laderman 算法很难检查/仔细检查，我将在下面发布完整的代码以供后代使用。

规格：Ubuntu 12.04、戴尔 Prevision T1600、gcc。时间百分比差异：

g++ [2.22, 2.23, 2.27]
g++ -O3 [-0.48, -0.49, -0.48]
g++ -funroll-loops -O3 [-0.48, -0.48, -0.47]

基准代码和 Laderman 实现：

#include <iostream>
#include <ctime>
#include <cstdio>
#include <cstdlib>
using namespace std;

void simple_mul(const double a[3][3], 
        const double b[3][3],
        double c[3][3]) {
  int i,j,m,n;
  for(i=0;i<3;i++) {
    for(j=0;j<3;j++) {
      c[i][j] = 0;
      for(m=0;m<3;m++) 
    c[i][j] += a[i][m]*b[m][j];
    }
  }
}

void laderman_mul(const double a[3][3], 
           const double b[3][3],
           double c[3][3]) {

   double m[24]; // not off by one, just wanted to match the index from the paper

   m[1 ]= (a[0][0]+a[0][1]+a[0][2]-a[1][0]-a[1][1]-a[2][1]-a[2][2])*b[1][1];
   m[2 ]= (a[0][0]-a[1][0])*(-b[0][1]+b[1][1]);
   m[3 ]= a[1][1]*(-b[0][0]+b[0][1]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][2]);
   m[4 ]= (-a[0][0]+a[1][0]+a[1][1])*(b[0][0]-b[0][1]+b[1][1]);
   m[5 ]= (a[1][0]+a[1][1])*(-b[0][0]+b[0][1]);
   m[6 ]= a[0][0]*b[0][0];
   m[7 ]= (-a[0][0]+a[2][0]+a[2][1])*(b[0][0]-b[0][2]+b[1][2]);
   m[8 ]= (-a[0][0]+a[2][0])*(b[0][2]-b[1][2]);
   m[9 ]= (a[2][0]+a[2][1])*(-b[0][0]+b[0][2]);
   m[10]= (a[0][0]+a[0][1]+a[0][2]-a[1][1]-a[1][2]-a[2][0]-a[2][1])*b[1][2];
   m[11]= a[2][1]*(-b[0][0]+b[0][2]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][1]);
   m[12]= (-a[0][2]+a[2][1]+a[2][2])*(b[1][1]+b[2][0]-b[2][1]);
   m[13]= (a[0][2]-a[2][2])*(b[1][1]-b[2][1]);
   m[14]= a[0][2]*b[2][0];
   m[15]= (a[2][1]+a[2][2])*(-b[2][0]+b[2][1]);
   m[16]= (-a[0][2]+a[1][1]+a[1][2])*(b[1][2]+b[2][0]-b[2][2]);
   m[17]= (a[0][2]-a[1][2])*(b[1][2]-b[2][2]);
   m[18]= (a[1][1]+a[1][2])*(-b[2][0]+b[2][2]);
   m[19]= a[0][1]*b[1][0];
   m[20]= a[1][2]*b[2][1];
   m[21]= a[1][0]*b[0][2];
   m[22]= a[2][0]*b[0][1];
   m[23]= a[2][2]*b[2][2];

  c[0][0] = m[6]+m[14]+m[19];
  c[0][1] = m[1]+m[4]+m[5]+m[6]+m[12]+m[14]+m[15];
  c[0][2] = m[6]+m[7]+m[9]+m[10]+m[14]+m[16]+m[18];
  c[1][0] = m[2]+m[3]+m[4]+m[6]+m[14]+m[16]+m[17];
  c[1][1] = m[2]+m[4]+m[5]+m[6]+m[20];
  c[1][2] = m[14]+m[16]+m[17]+m[18]+m[21];
  c[2][0] = m[6]+m[7]+m[8]+m[11]+m[12]+m[13]+m[14];
  c[2][1] = m[12]+m[13]+m[14]+m[15]+m[22];
  c[2][2] = m[6]+m[7]+m[8]+m[9]+m[23];    
}

int main() {
  int N = 1000000000;
  double A[3][3], C[3][3];
  std::clock_t t0,t1;
  timespec tm0, tm1;

  A[0][0] = 3/5.; A[0][1] = 1/5.; A[0][2] = 2/5.;
  A[1][0] = 3/7.; A[1][1] = 1/7.; A[1][2] = 3/7.;
  A[2][0] = 1/3.; A[2][1] = 1/3.; A[2][2] = 1/3.;

  t0 = std::clock();
  for(int i=0;i<N;i++) {
    // A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
    simple_mul(A,A,C);
  }
  t1 = std::clock();
  double tdiff_simple = (t1-t0)/1000.;

  cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl;
  cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl;
  cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl;
  cout << tdiff_simple << endl;
  cout << endl;

  t0 = std::clock();
  for(int i=0;i<N;i++) {
    // A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
    laderman_mul(A,A,C);
  }
  t1 = std::clock();
  double tdiff_laderman = (t1-t0)/1000.;

  cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl;
  cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl;
  cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl;
  cout << tdiff_laderman << endl;
  cout << endl;

  double speedup = (tdiff_simple-tdiff_laderman)/tdiff_laderman;
  cout << "Approximate speedup: " << speedup << endl;

  return 0;
}

【讨论】：

您是否尝试过使用更多优化选项，例如 -funroll-loops 或使用 sse？
不，但我会试一试（我是优化新手），您还有什么建议吗？
玩各种诸如命令行选项之类的东西（只需阅读手册页并尝试一下。gcc 也有一个设置优化标志的属性，所以你可以破解一些可以测试多个这样的标志的东西在一个二进制文件中）和gcc.godbolt.org 上的 __restrict gcc 扩展通常可以很好地了解正在发生的事情。
@PlasmaHH 我已经编辑了答案以反映展开的时间，没有真正的变化。我怀疑我实现简单版本的任何标志最终可能会更快，尽管我会尝试一下。

【解决方案3】：

我预计主要的性能问题将是内存延迟。 double[9] 通常为 72 个字节。这已经是一笔不小的数目，而您正在使用其中的三个。

【讨论】：

72 字节是否太大而无法放入一级缓存？我的印象是你可以将千字节的数据放入缓存中，而且速度非常快。我对CPU中各种操作的时序了解不多；我在这里错过了什么吗？
@Dan：它很容易适应，但会占用多个缓存行。主要关注的是当矩阵不在缓存中时会发生什么。在这种情况下，您可能会从主内存加载 2*64 字节（因为缓存读取整行）。

【解决方案4】：

虽然问题提到了 C++，但我在 C# (.NET 4.5) 中实现了 3x3 矩阵乘法 C=A*B，并在我的 64 位 windows 7 机器上运行了一些基本的时序测试并进行了优化。 10,000,000 次乘法大约花费了

0.556 秒，简单的实现和
0.874 秒，来自另一个答案的拉德曼代码。

有趣的是，laderman 代码比简单的方式慢。我没有使用分析器进行调查，但我猜额外分配比一些额外的乘法更昂贵。

似乎当前的编译器足够聪明，可以为我们进行这些优化，这很好。为了您的兴趣，这是我使用的幼稚代码：

    public static Matrix3D operator *(Matrix3D a, Matrix3D b)
    {
        double c11 = a.M11 * b.M11 + a.M12 * b.M21 + a.M13 * b.M31;
        double c12 = a.M11 * b.M12 + a.M12 * b.M22 + a.M13 * b.M32;
        double c13 = a.M11 * b.M13 + a.M12 * b.M23 + a.M13 * b.M33;
        double c21 = a.M21 * b.M11 + a.M22 * b.M21 + a.M23 * b.M31;
        double c22 = a.M21 * b.M12 + a.M22 * b.M22 + a.M23 * b.M32;
        double c23 = a.M21 * b.M13 + a.M22 * b.M23 + a.M23 * b.M33;
        double c31 = a.M31 * b.M11 + a.M32 * b.M21 + a.M33 * b.M31;
        double c32 = a.M31 * b.M12 + a.M32 * b.M22 + a.M33 * b.M32;
        double c33 = a.M31 * b.M13 + a.M32 * b.M23 + a.M33 * b.M33;
        return new Matrix3D(
            c11, c12, c13,
            c21, c22, c23,
            c31, c32, c33);
    }

其中 Matrix3D 是不可变结构（只读双字段）。

棘手的事情是提出一个有效基准测试，您可以在其中衡量您的代码，而不是编译器对您的代码做了什么（带有大量额外内容的调试器，或者在没有您实际的情况下进行优化）代码，因为从未使用过结果）。我通常会尝试“触摸”结果，这样编译器就无法删除被测代码（例如，检查矩阵元素是否与 89038.8989384 相等，如果相等则抛出）。但是，最后我什至不确定编译器是否会破解这个比较:)

【讨论】：