Laderman的3x3矩阵乘法只有23次乘法,值得吗？

Question

取两个3x3矩阵的乘积A*B=C.天真地,这需要使用标准算法进行 27次乘法.如果一个人很聪明,你可以只使用23次乘法来做到这一点,这是拉德曼于1973年发现的结果.该技术涉及保存中间步骤并以正确的方式组合它们.

现在让我们修改一个语言和一个类型,比如说C++的元素double.如果Laderman算法是硬编码而不是简单的双循环,那么我们是否可以期望现代编译器的性能能够消除算法的差异？

关于这个问题的注释:这是一个编程站点,问题是在时间关键内循环的最佳实践的上下文中提出的; 过早优化这不是.关于实施的提示非常受欢迎.

Answer 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个这些计算.

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

Answer 2

时间测试:

我自己进行了计时测试,结果让我感到惊讶(因此我首先提出问题的原因).缺点是,在标准编译下,laderman速度提高了约225%,但是使用-03优化标志时速度要慢50%!每次在-O3标志期间我都必须在矩阵中添加一个随机元素,或者编译器完全优化掉简单乘法,在时钟精度内采用零时间.由于laderman算法很难检查/仔细检查,我会在下面发布完整的代码供后人使用.

规格:Ubuntu 12.04,Dell 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;
}