週刊 spaghetti_source

Totally Monotone Matrix Searching (SMAWK algorithm)

00:25

前回と引き続き，普通にやると O(n^2) だが，ある性質を満たすと O(n) パターンのアルゴリズムを紹介します．

Monotone / Totally Monotone

• f(i,j) がmonotoneであるとは，i < i' のとき argmin f(i,*) < argmin f(i',*) を満たすことをいう．
• f(i,j) がtotally monotoneであるとは，任意の部分行列がmonotoneであることをいう．

要するに，monotoneというのは「行ごとの最小値の配置が（行列として見たとき）右下に単調に下がっていく」もので，totally monotoneはそれの強いバージョン（部分行列に限ったとき，その中の行最小値配置が単調）です．以下totally monotoneを簡単にTMと略します．当然ですがtotally monotoneならばmonotoneです．

前回のKnuth-Yao speedupで用いたMonge性との関係ですが，「fがMongeならfはTM」という性質が成り立ちます（逆は不成立）．TM性よりもMonge性のほうが判定容易なこともあり，応用上はMonge性を示してTM性とすることがよくあります．

また，最近Bein-Golin-Larmore-Zhang(2009)が，KY-speedupが適用できる問題にはTM-speedup（SMAWK）が適用できる，という結果を示しているため，今回の内容と前回の内容は，本質的には同じものといえるかもしれません．

Aggarwal-Klawe-Moran-Shor-Wilber(1987)は，これらの概念を提案し，TM行列の行最小値配置を求めるO(n+m)アルゴリズムを構築しました．以下これを説明します．元論文では最大値配置でやっていますが，同じ事なので，自分の好みで最小値配置でやります．

1. monotone minima

まず，SMAWKとは直接関係ありませんが，（totallyでない）monotone fの行最小値配置を求めるアルゴリズムから考えます．monotoneはTMよりも弱い条件なので，O(n+m)よりも時間がかかります．

monotone minimaの手続きは，以下の分割統治で与えられます．

1. n/2 (:= im) 行目の最小値を線形探索し，座標を jm とする．
2. (im,jm) の左上と，(im,jm) の右下部分について再帰的に最小値配置を求める．

正当性は「monotone = 最小値配置が右下に下がっていく」のだから，ある行で最小値座標が判明すれば，その右上と左下は調べる必要がない，というだけです．

計算量は f(n,m) = f(n/2,m1) + f(n/2,m2) + O(m) (m1+m2=m-1) を解いて O(m log n) ．実はこれは理論下界でもあります．

int minima[MAXSIZE];
void rowMinima(int ib, int ie, int jb, int je) {
  if (ib  > ie) return;
  if (ib == ie) {
    int jm = jb;
    int fm = f(ib,jm), fj;
    for (int j = jb+1; j <= je; ++j) 
      if (fm > (fj = f(ib, j))) { fm = fj; jm = j; }
    minima[ib] = jm;
    return;
  }
  int im = (ib+ie)/2;
  rowMinima(im, im, jb, je);
  rowMinima(ib, im-1, jb, minima[im]);
  rowMinima(im+1, ie, minima[im], je);
}

2. TM minima (SMAWK algorithm)

次にTM性を仮定してより高速なアルゴリズムを設計します．以下の手続きが，提案者たちの頭文字をとって，現在SMAWKアルゴリズムと呼ばれているものです．Larmore先生 の http://www.egr.unlv.edu/~larmore/Courses/CSC477/monge.pdf が非常にわかりやすい解説（具体例）なので，自分の記事よりもこのpdfを読むのが良い説があります．

SMAWKアルゴリズムも基本的には分割統治ですが，monotone minimaよりも巧妙です．手続きは初期化を含めて以下の４ステップから構成されます．

Initialize

row = {1,2,...,n}, col = {1,2,...,m}

Reduce

col のうち最小値を含む可能性のあるものを取り出して col' とする．

row を 1 つ飛ばし（偶数番目のみ抽出）で作った集合を row' とする．

Recursion

row', col' からなる部分行列について再帰的に行最小値配置を求める（TM行列の部分行列もTM）．

Interpolate

最小値配置の求まっていない部分（row-row'）について最小値配置を決定する．

以下これらの内容を説明します．InitializeとRecursionは説明不要でしょう．

Reduce

SMAWKアルゴリズムの一番の要所がReduceです．原理は以下の２つの性質に基づきます．

• (a) もし f(i,j) < f(i,k) (j < k) ならば f(i',j) (i' > i) は捨てて良い．
• (b) もし f(i,j) > f(i,k) (j < k) ならば f(i',k) (i' < i) は捨てて良い．

捨てた部分に行最小値があるかもしれませんが，行最小値がある列が全て捨てられることは無いので問題ありません．この２つの性質を使って要素をどんどん捨てていき，全部の要素が捨てられた列を捨てることで，列の間引きを行います．

実際のアルゴリズムは以下のようになります．

1. 列集合を管理するスタック S を用意し，１列目から順番に列を追加しようと試みる
2. k 列目を追加する際，j = S.top について，f(i,j) ≧ f(i,k) ならば j 列目を捨てる（i = |S|）
3. |S| < |row| ならば S に k を追加する

S に k を追加する ＝ f(i,j) < f(i,k) が成立しているということなので，S に追加された列は (b) の条件で列の上側が捨てられています．さらに，あとでその列が比較されて 2 の条件が満たされたときは (a) の条件で列の下側も捨てられ，結局すべての要素が捨てられたことになります（よって，Sからpopしてよい）．

Reduce自体は O(n+m) で実行でき，Reduceが終わったときの |S| は高々行数と同じとなることが計算量評価で重要となります．

Interpolate

偶数番目の最小値配置がわかっているときに，奇数番目の最小値配置を求めるのがInterpolateです．これは前回のMonge性を使う加速と同じく，探索の上下界がわかっていることを利用します．

奇数番目の行 i の列最小値の存在範囲は，TM性より i-1行目 の列最小値の位置より右で，i+1行目の列最小値の位置より左です．よってすべての奇数番目の行について，各範囲で線形探索しても，合計で O(n+m) しかかかりません．

計算量評価

ReduceとInterpolateでO(n+m)，再帰のパラメタは(n/2,n) です．すなわち f(n,m) = f(n/2,n) + O(n+m) です．monotone minimaの再帰式との違いは，再帰のパラメタがnになっている点で，このおかげで合計計算量が O(n+m) になります．正確な証明は帰納法によります．

void SMAWK(vector<int> js, int ib, int ie, int id) {
  if (ib > ie) return;

  vector<int> js2;
  for (int q = 0, i = ib; q < js.size(); ++q) {
    while (!js2.empty() && f(i, js2.back()) >= f(i, js[q])) {
      js2.pop_back(); i -= id; }
    if ( js2.size() != (ie-ib)/id ) {
      js2.push_back(js[q]); i += id; }
  }

  SMAWK(js2, ib+id, ie, id*2);

  for (int i = ib, q = 0; i <= ie; i += id*2) {
    int jt = (i + id <= ie ? minima[i+id] : js.back());
    int fm = 99999999, fq;
    for (; q < js.size(); ++q) {
      if ( fm > (fq = f(i, js[q])) ) {
        fm = fq;
        minima[i] = js[q];
      }
      if (js[q] == jt) break;
    }
  }
}
void rowMinimaTM(int ib, int ie, int jb, int je) {
  vector<int> js;
  for (int j = jb; j <= je; ++j) js.push_back(j);
  SMAWK(js, ib, ie, 1);
}

実測

この手の物だと「複雑だから定数項がでかくて O(n) より O(n log n) のほうが実用範囲では早い」なんてことがよくあるので，実測が必須です．

例題として，次の「最近点問題」の特殊版を考えます．

・入力：平面上の点集合 P, Q，ただし P は x 軸上に乗っている

・出力：各 p ∈ P に対する最近点 q(p) ∈ Q

f(p,q) := d(p,q) とおくと，P, Q が x 座標でソートされているとき，f は Totally Monotone となり，したがって SMAWK が適用できます．この問題で P, Q をランダム生成したときの計算時間を以下に示します．この結果より，計算量でまさる SMAWK のほうが実際に高速であることがわかります．

n	Mono	SMAWK
500,000	0.109	0.078
1,000,000	0.234	0.172
5,000,000	1.232	0.827
10,000,000	2.699	1.747

とはいえ，有利なのは1.5倍程度，競技プログラミングでこの差が問題になることはまず無いので，理解が容易でシンプル，かつ適用範囲が（少しだけ）広いMonotone minimaだけメモしておいても十分事足りるかと思います．