週刊 spaghetti_source

Matrix Chain Multiplication / Optimal triangulation (approximation)

17:13

いつか理解して簡単化したいと思っていながら全然進まないMCMの話題です．

連鎖行列積問題（matric chain multiplication）とは，適当なサイズの長方行列 M1, M2, ... が与えられたとき，M1 x M2 x ... x Mn を計算するための最良の方法（カッコ付け）を求める問題です．ここで p x q 行列と q x r 行列の積は pqr で計算できるとします．たとえば M1 = 1x100, M2 = 100x100, M3 = 100x10 だった場合，(M1 M2) M3 だと 10000 + 1000 の計算コスト，M1 (M2 M3) だと 100000 + 1000 の計算コストになります．この場合は前者が最適解で，次善の解と比べて約10倍の差になります．

この問題は動的計画法の説明で用いられ，補足として「実は O(n log n) で解ける[Hu-Shing 1982 & 1984」とコメントされるのが様式美です．O(n log n) は比較モデルのもとで最善です [Bradford-Choppella-Rawlins 1993]．

本当はここで Hu-Shing の O(n log n) を紹介したいのですが，原論文がややこしくて，現状積ん読状態です．（誰か一緒に読みませんか？）……私見ですが，現在 O(n log n) を厳密に理解している人はほとんど居ないのではないか，と思っています（結果に触れる論文は多いが，内容について触れる論文が全く無い）．

そこで今回は近似アルゴリズムを紹介します．

連鎖行列積問題のままでは扱いにくいので，これを最適三角形分割（optimal triangulation）とよばれる次の問題に帰着します：正n角形の各頂点に値 w(1), w(2), ..., w(n) が対応しているものが与えられたとき，これの三角形分割であってコストが最小のものを求めよ．ここで三角形分割のコストは現れる三角形 (i,j,k) すべてについて f(i,j,k) := w(i) w(j) w(k) を足しあわせたものとして定義する．

連鎖行列積問題と最適三角形分割問題が同値であることを冒頭の例を変形することで確認しておきます．1x100, 100x100, 100x10 行列に対して重みが 1, 100, 100, 10 なる四角形を対応させます．四角形の三角形分割は２種類あり，それらが連鎖行列積の (M1 M2) M3 と M1 (M2 M3) に対応します．

さて，連鎖行列積問題の近似アルゴリズムは次のアイデアに基づきます：

• (1) とても横長の行列があれば，先に計算する．
  • 例： 3x4, 4x6, 6x10000, 10000x3, 3x2 なら最初に (6x10000) x (10000x3) を計算する．
• (2) そのようなものがなければ小さい方から貪欲に計算する．
  • 例： 2x4, 4x3, 3x5, 5x4 なら，2 が一番多く評価に現れるように (((2x4) x (4x3)) x (3x5)) x (5x4) とする．

これを最適三角形分割の言葉で言い直すと次のようになります：

• (1) とても重い点があれば，その両側を結ぶ．
• (2) そのようなものがなければ一番軽い点から扇状に分割する．

どの程度重ければ左右を直結してよいかが問題となりますが，次の定理を使います．

Theorem. 任意の添字 * について 1/w(*) + 1/w(i) ≦ 1/w(i-1) + 1/w(i+1) が成り立つとき (i-1) と (i+1) を直結する最適解が存在する．

証明：三角形分割を１つ固定し，その中の四角形 (i-1), i, (i+1), p を考える．三角形分割なので辺は (i-1)-(i+1) もしくは i-p である．前者であるための条件は，前者の分割コスト ≦ 後者の分割コストであり，整理すると 1/w(p) + 1/w(i) < 1/w(i-1) + 1/w(i+1) を得る．最適三角形分割ではどの頂点 p が i の対角点になるかわからないが，すべての頂点でこの条件が成り立っていれば十分．

上の不等式は，* として w(s) 最小の点を取るのが一番厳しいことは容易にわかります．これをあわせて，次のアルゴリズムを得ます[Chin 78, Hu-Shing 81]

s := argmin w(*)
while 不等式を満たす頂点 i が存在
  i-1, i+1 を結んで三角形 (i-1,i,i+1) を切り落とす
s を軸とする扇を出力

不等式を満たす i を探す部分の計算量がかかるように思えますが，実のところ，前から順番に見ていくだけで十分です．簡単のため多角形を回転して s = 1 にとります．S を「最終的に扇の羽になる頂点」の集合をあらわすスタックとすると次の手続きが得られます．

S = []
for i = 1, 2, ...
  while 1, Sの先頭2つ, i が不等式を満たす
    S の先頭要素を削除
  S に i を追加
1 を軸とする扇を出力

この手続きの誤差率は1/(2√3+3)≒0.1547...となります．これは，この手続きに対する最悪インスタンスを構成して最悪の値を取ったものです．計算量は同じ点が 2 回スタックに入らないので O(n)です．

※歴史的には，この手続きは Chin が（近似率の悪い）バージョンを提案し，それを見たHu and Shingがすぐに改善できることに気づいた，という流れになっています．

ランダム入力などでやってみると，かなりうまく合うことがわかります．以下の実装を用いた自分の実験では，n = 100, w[i] \in [100,150] uniform のランダム列に対して平均で近似率 1.005 程度の解が出ています．

これだけ簡単な手続きでこれだけ良い近似が出るのだから，厳密アルゴリズムでこれを使って枝刈りするとうまくいかないかなあ，と一時期思っていましたが，あまりうまく行かないモノですねえ．

実装例

typedef long long Int;
bool cond(Int a, Int b, Int c, Int d) {
  // 1/a + 1/c <= 1/b + 1/d  iff b-d separation
  return b*c*d + a*b*d <= a*c*d + a*b*c;
}
Int approxTriangulation(vector<Int> w) {
  int n = w.size(), k = 0, S[n];
  rotate(w.begin(), min_element(w.begin(), w.end()), w.end());
  Int cost = 0;
  for (int i = 1; i < n; ++i) {
    while (k >= 2 && cond(w[0], w[S[k-2]], w[S[k-1]], w[i])) {
      cost += w[S[k-2]] * w[S[k-1]] * w[i];
      --k;
    }
    S[k++] = i;
  }
  for (int i = 0; i+1 < k; ++i) 
    cost += w[0] * w[S[i]] * w[S[i+1]];
  return cost;
}

参考文献

• L. E. Deimel, Jr and T. A. Lampe (1979): "An invariance theorem concerning optimal computation of matrix chain products", North Carolina State University Technical Report #TR79-14.
• F. Y. Chin (1978): "An O(n) algorithm for determining a near optimal computation order of matrix chain products", Communications of the ACM, vol.21, no.7, vol.544-549.
• T. C. Hu and M. T. Shing (1981): "An O(n) algorithm to find a near optimal partition of a conve polygon", Journal of Algorithms, vol.2, pp.122-138.
• T. C. Hu and M. T. Shing (1982): "Computation of matrix chain products, Part I" SIAM Journal on Computing, vol.11, no.2, pp.362-373.
• T. C. Hu and M. T. Shing (1984): "Computation of matrix chain products, Part II" SIAM Journal on Computing, vol.13, no.2, pp.241-261.
• P. G. Bradford, V. Choppell, and G. J. E. Rawlins: "On Lower Bounds for the Matrix Chain Ordering Problem", Indiana University, Technical Report TR391.