行列計算が速くなる？思ってたのと違った

Google DeepMindのAlphaEvolveが、より効率的な行列乗算アルゴリズム見つけたと思ったら、思ってたのと違いました。

CG分野では、4x4行列を多用するので、これの乗算が速くなるのかなと思ったんですが、複素数行列においてとか、条件付きでした。AIとかで使う、とても大きな行列で使う特殊な行列演算が対象のようです。

単純な4x4行列の乗算では、64回の乗算と64回の加算をする。
それが48回になるのかなと思ったけど違う話でした。

一般的に乗算命令は重くて、加算は軽いと言われるので、乗算を減らしたほうがいいのだが、現実的には行列演算用のCPUやGPUの命令があるので、それを使ったほうが速い。専用クラス使ったほうが最適化され速い。

でも、気になったので、色々調べていたら、次の論文を見つけた。
こっちのほうが、47回で、なんだろうと思いながら、

https://www.nature.com/articles/s41586-022-05172-4

論文内の画像でExtended Data Fig. 1: Algorithm for multiplying 4 × 4 matrices in modular arithmetic (Z 2) with 47 multiplications.とある。

https://www.nature.com/articles/s41586-022-05172-4/figures/6

よく意味もわからず、実装したが、結局のところ、mod 2、つまり、0と1の世界の行列乗算が47回の乗算でできるという話のようだった。乗算は減るけど、加算が340回。とんでもねぇ。

この論文の数式が画像埋め込みで、OCRとか駆使しながらなんとか、がんばって実装したので、記念においておく。

さらに、ニュースで比較に出てくるStrassen法のアルゴリズム（49回の乗算）も貼っておく。

```java:Algorithm for multiplying 4 × 4 matrices in modular arithmetic (Z 2) with 47 multiplications
void setup() {
  float[][] m = {{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
  float[][] n = {{1, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};

  m = alphaTensor47Mult(m, n);
  println("47 Mult =");
  println(m[0][0], m[0][1], m[0][2], m[0][3]);
  println(m[1][0], m[1][1], m[1][2], m[1][3]);
  println(m[2][0], m[2][1], m[2][2], m[2][3]);
  println(m[3][0], m[3][1], m[3][2], m[3][3]);

}

float[][] alphaTensor47Mult(float[][] a, float[][] b) {
  float[] h = new float[48]; // h[1]〜h[47]のみ使う
  float[][] c = new float[4][4];

  h[1]  = a[0][0] * b[0][2];
  h[2]  = (a[0][0] + a[2][0] + a[2][2]) * (b[0][0] + b[2][0] + b[2][2]);
  h[3]  = (a[0][0] + a[2][0] + a[2][3]) * (b[0][1] + b[3][1] + b[3][2]);
  h[4]  = (a[0][2] + a[1][0] + a[1][2]) * (b[0][2] + b[0][3] + b[2][3]);
  h[5]  = (a[0][0] + a[2][0]) * (b[0][0] + b[0][1] + b[0][2] + b[2][0] + b[2][2] + b[3][1] + b[3][2]);
  h[6]  = (a[0][2] + a[1][2]) * (b[0][2] + b[0][3] + b[2][1] + b[2][2] + b[2][3] + b[3][1] + b[3][2]);
  h[7]  = (a[0][3] + a[3][2] + a[3][3]) * (b[2][0] + b[2][2] + b[3][0]);
  h[8]  = (a[0][3] + a[3][0] + a[3][3]) * (b[0][2] + b[0][3] + b[3][3]);
  h[9]  = (a[0][2] + a[1][2] + a[1][3]) * (b[2][1] + b[3][1] + b[3][2]);
  h[10] = (a[0][3] + a[3][3]) * (b[0][2] + b[0][3] + b[2][0] + b[2][2] + b[3][0] + b[3][2] + b[3][3]);
  h[11] = a[2][2] * (b[0][0] + b[1][1] + b[1][2] + b[2][0] + b[2][1]);
  h[12] = (a[0][1] + a[2][1] + a[2][2]) * (b[1][1] + b[1][2] + b[2][1]);
  h[13] = a[2][3] * (b[0][1] + b[1][0] + b[1][2] + b[3][0] + b[3][1]);
  h[14] = (a[0][1] + a[2][1]) * (b[1][0] + b[1][1] + b[1][2] + b[2][1] + b[3][0]);
  h[15] = (a[0][1] + a[2][1] + a[2][3]) * (b[1][0] + b[1][2] + b[3][0]);
  h[16] = a[1][0] * (b[0][1] + b[0][3] + b[1][1] + b[1][2] + b[2][3]);
  h[17] = (a[0][1] + a[1][0] + a[1][1]) * (b[0][1] + b[1][1] + b[1][2]);
  h[18] = (a[0][1] + a[1][1]) * (b[0][1] + b[1][1] + b[1][2] + b[1][3] + b[3][3]);
  h[19] = a[1][3] * (b[1][2] + b[1][3] + b[2][1] + b[3][1] + b[3][3]);
  h[20] = (a[0][1] + a[1][2] + a[1][3] + a[2][1] + a[2][2]) * b[2][1];
  h[21] = (a[0][1] + a[1][1] + a[1][3]) * (b[1][2] + b[1][3] + b[3][3]);
  h[22] = a[3][2] * (b[1][2] + b[1][3] + b[2][0] + b[2][3] + b[3][0]);
  h[23] = (a[0][0] + a[0][2] + a[0][3] + a[1][2] + a[1][3] + a[2][0] + a[2][3]) * (b[3][1] + b[3][2]);
  h[24] = (a[0][1] + a[3][1] + a[3][2]) * (b[1][2] + b[1][3] + b[2][3]);
  h[25] = (a[0][1] + a[3][1]) * (b[0][0] + b[1][0] + b[1][2] + b[1][3] + b[2][3]);
  h[26] = (a[0][1] + a[3][0] + a[3][1]) * (b[0][0] + b[1][0] + b[1][2]);
  h[27] = a[0][3] * b[3][2];
  h[28] = (a[0][1] + a[1][0] + a[1][1] + a[2][0] + a[2][3]) * b[0][1];
  h[29] = (a[0][1] + a[1][0] + a[1][2] + a[3][1] + a[3][2]) * b[2][3];
  h[30] = (a[0][1] + a[2][0] + a[2][2] + a[3][0] + a[3][1]) * b[0][0];
  h[31] = a[3][0] * (b[0][0] + b[0][3] + b[1][0] + b[1][2] + b[3][3]);
  h[32] = (a[0][1] + a[2][1] + a[2][3] + a[3][2] + a[3][3]) * b[3][0];
  h[33] = (a[0][1] + a[1][1] + a[1][3] + a[3][0] + a[3][3]) * b[3][3];
  h[34] = (a[1][0] + a[2][0] + a[3][0]) * (b[0][0] + b[0][1] + b[0][3]);
  h[35] = (a[0][1] + a[1][0] + a[1][1] + a[2][1] + a[2][2]) * (b[1][1] + b[1][2]);
  h[36] = (a[0][1] + a[1][3] + a[2][1] + a[3][2]) * (b[1][2] + b[1][3] + b[2][1] + b[3][0]);
  h[37] = (a[0][1] + a[1][0] + a[2][2] + a[3][1]) * (b[0][0] + b[1][1] + b[1][2] + b[2][3]);
  h[38] = (a[1][1] + a[2][1] + a[3][1]) * (b[1][0] + b[1][1] + b[1][3]);
  h[39] = a[0][1] * b[1][2];
  h[40] = a[0][2] * b[2][2];
  h[41] = (a[0][0] + a[0][2] + a[0][3] + a[1][0] + a[1][2] + a[3][0] + a[3][3]) * (b[0][2] + b[0][3]);
  h[42] = (a[0][1] + a[2][1] + a[2][3] + a[3][0] + a[3][1]) * (b[1][0] + b[1][2]);
  h[43] = (a[1][3] + a[2][3] + a[3][3]) * (b[3][0] + b[3][1] + b[3][3]);
  h[44] = (a[1][2] + a[2][2] + a[3][2]) * (b[2][0] + b[2][1] + b[2][3]);
  h[45] = (a[0][0] + a[0][2] + a[0][3] + a[2][0] + a[2][2] + a[3][2] + a[3][3]) * (b[2][0] + b[2][2]);
  h[46] = (a[0][1] + a[1][1] + a[2][3] + a[3][0]) * (b[0][1] + b[1][0] + b[1][2] + b[3][3]);
  h[47] = (a[0][1] + a[1][1] + a[1][3] + a[3][1] + a[3][2]) * (b[1][2] + b[1][3]);

  c[0][0] = h[15] + h[26] + h[2] + h[30] + h[32] + h[39] + h[40] + h[42] + h[45] + h[7];
  c[1][0] = h[11] + h[12] + h[14] + h[20] + h[22] + h[24] + h[25] + h[29] + h[35] + h[36] + h[37] + h[38] + h[44] + h[47];
  c[2][0] = h[11] + h[12] + h[14] + h[15] + h[26] + h[30] + h[39] + h[42];
  c[3][0] = h[15] + h[22] + h[24] + h[25] + h[26] + h[32] + h[39] + h[42];

  c[0][1] = h[12] + h[17] + h[20] + h[23] + h[27] + h[28] + h[35] + h[39] + h[3] + h[9];
  c[1][1] = h[12] + h[17] + h[18] + h[19] + h[20] + h[21] + h[35] + h[39];
  c[2][1] = h[12] + h[13] + h[14] + h[15] + h[17] + h[28] + h[35] + h[39];
  c[3][1] = h[13] + h[14] + h[15] + h[18] + h[19] + h[21] + h[32] + h[33] + h[36] + h[38] + h[42] + h[43] + h[46] + h[47];

  c[0][2] = h[1] + h[27] + h[39] + h[40];
  c[1][2] = h[16] + h[17] + h[18] + h[19] + h[21] + h[39] + h[40] + h[4] + h[6] + h[9];
  c[2][2] = h[11] + h[12] + h[13] + h[14] + h[15] + h[1] + h[2] + h[39] + h[3] + h[5];
  c[3][2] = h[10] + h[22] + h[24] + h[25] + h[26] + h[27] + h[31] + h[39] + h[7] + h[8];

  c[0][3] = h[1] + h[21] + h[24] + h[29] + h[33] + h[39] + h[41] + h[47] + h[4] + h[8];
  c[1][3] = h[16] + h[17] + h[18] + h[21] + h[24] + h[29] + h[39] + h[47];
  c[2][3] = h[16] + h[17] + h[18] + h[25] + h[26] + h[28] + h[30] + h[31] + h[34] + h[35] + h[37] + h[38] + h[42] + h[46];
  c[3][3] = h[21] + h[24] + h[25] + h[26] + h[31] + h[33] + h[39] + h[47];

  return c;
}
```


```java:strassen 4x4
void setup() {
  float[][] m = {
    {1, 2, 3, 4},
    {5, 6, 7, 8},
    {9, 10, 11, 12},
    {13, 14, 15, 16}
  };
  float[][] n = {
    {17, 18, 19, 20},
    {21, 22, 23, 24},
    {25, 26, 27, 28},
    {29, 30, 31, 32}
  };

  float[][] r = strassen4x4(m, n);
  println("Strassen 4x4 Mult =");
  for (int i = 0; i < 4; i++) {
    println(nf(r[i][0], 1, 2), nf(r[i][1], 1, 2), nf(r[i][2], 1, 2), nf(r[i][3], 1, 2));
  }
}

// 2×2行列Strassen積
float[][] strassen2x2(float[][] a, float[][] b) {
  float[][] res = new float[2][2];

  float m1 = (a[0][0] + a[1][1]) * (b[0][0] + b[1][1]);
  float m2 = (a[1][0] + a[1][1]) * b[0][0];
  float m3 = a[0][0] * (b[0][1] - b[1][1]);
  float m4 = a[1][1] * (b[1][0] - b[0][0]);
  float m5 = (a[0][0] + a[0][1]) * b[1][1];
  float m6 = (a[1][0] - a[0][0]) * (b[0][0] + b[0][1]);
  float m7 = (a[0][1] - a[1][1]) * (b[1][0] + b[1][1]);

  res[0][0] = m1 + m4 - m5 + m7;
  res[0][1] = m3 + m5;
  res[1][0] = m2 + m4;
  res[1][1] = m1 - m2 + m3 + m6;

  return res;
}

// 4×4行列 Strassen法2段適用
float[][] strassen4x4(float[][] a, float[][] b) {
  float[][] res = new float[4][4];

  // 2x2ブロックに分割
  float[][] a00 = new float[2][2], a01 = new float[2][2], a10 = new float[2][2], a11 = new float[2][2];
  float[][] b00 = new float[2][2], b01 = new float[2][2], b10 = new float[2][2], b11 = new float[2][2];

  for (int i = 0; i < 2; i++) {
    for (int j = 0; j < 2; j++) {
      a00[i][j] = a[i][j];
      a01[i][j] = a[i][j + 2];
      a10[i][j] = a[i + 2][j];
      a11[i][j] = a[i + 2][j + 2];

      b00[i][j] = b[i][j];
      b01[i][j] = b[i][j + 2];
      b10[i][j] = b[i + 2][j];
      b11[i][j] = b[i + 2][j + 2];
    }
  }

  // Strassen法再帰
  float[][] m1 = strassen2x2(add(a00, a11), add(b00, b11));
  float[][] m2 = strassen2x2(add(a10, a11), b00);
  float[][] m3 = strassen2x2(a00, sub(b01, b11));
  float[][] m4 = strassen2x2(a11, sub(b10, b00));
  float[][] m5 = strassen2x2(add(a00, a01), b11);
  float[][] m6 = strassen2x2(sub(a10, a00), add(b00, b01));
  float[][] m7 = strassen2x2(sub(a01, a11), add(b10, b11));

  // c00
  float[][] c00 = add(sub(add(m1, m4), m5), m7);
  // c01
  float[][] c01 = add(m3, m5);
  // c10
  float[][] c10 = add(m2, m4);
  // c11
  float[][] c11 = add(sub(add(m1, m3), m2), m6);

  // 合成
  for (int i = 0; i < 2; i++) {
    for (int j = 0; j < 2; j++) {
      res[i][j] = c00[i][j];
      res[i][j + 2] = c01[i][j];
      res[i + 2][j] = c10[i][j];
      res[i + 2][j + 2] = c11[i][j];
    }
  }
  return res;
}

// 行列和
float[][] add(float[][] a, float[][] b) {
  float[][] r = new float[a.length][a[0].length];
  for (int i = 0; i < a.length; i++) for (int j = 0; j < a[0].length; j++) r[i][j] = a[i][j] + b[i][j];
  return r;
}
// 行列差
float[][] sub(float[][] a, float[][] b) {
  float[][] r = new float[a.length][a[0].length];
  for (int i = 0; i < a.length; i++) for (int j = 0; j < a[0].length; j++) r[i][j] = a[i][j] - b[i][j];
  return r;
}
```