週刊 spaghetti_source

Dinkelbach algorithm

09:36

関数 f(x) と g(x) の比 f(x)/g(x) の最小化問題のことを分数計画（fractional programming），比の最小化などとよびます．発散を防ぐため g(x) > 0 を仮定します．今回紹介する Dinkelbach algorithm は，連続・離散どちらについても，広い範囲で適用できる比の最小化に対するアルゴリズムです．

比の最小化に対する基本的な方針は，[Jagannathan 1966] と [Dinkelbach 1967] が確立したパラメトリックアプローチです．これは t* = min f(x)/g(x) と思って分母を払う気持ちで h(t) = min { f(x) - t g(x) | x \in X } というパラメタ t をもつ関数を考えると

• h(t) > 0 iff t は最適値より大
• h(t) = 0 iff t は最適値
• h(t) < 0 iff t は最適値より小

という関係がなりたつ，というものです．この関係から，比の最小化は h(t) の零点計算に帰着するのですが，h(t) は狭義単調減少する凹関数なので（∵ h は「x でパラメタ付けられた傾きが負の線型関数」の min），零点計算がわりと簡単に行えます．h(t) の関数値評価では x に関する最小化問題を解く必要があり，一般にはもとの最小化とどちらが簡単かは言えませんが，よく扱われる問題では h(t) の関数値評価のほうが簡単です（eg., f, g が両方共線型なら f - t g も線型，f 凸, g 凹, t ≧ 0 なら f - t g 凸など）．

さて，h(t) = 0 を満たす t の計算法で主流のものが 3 つあります．

• 二分探索 [see: Radzik 1992]
• Dinkelbach algrithm [Dinkelbach 1967]
• Parametric search [Meggido 1979]

二分探索は単純に h(t) = 0 を二分探索するもので普通ですが，特に早くはありません．Parametric Search は二分探索を t を不定元だと考えて実行して t の分岐木を作るタイプのアルゴリズムで，理論計算量では最良のものができることが多いのですが，実装が複雑化して定数が大きい [Cole 1987] と言われています．最近は少し実用的になったみたい [van Oostrum-Veltkamp 2004] ですが，自分の印象では，あまり実装したくないタイプです（※ この論文はタイトルが「Parametric search made practical」．つまり，実用的になったよ！で論文が書ける程度には複雑）．

Dinkelbach algorithm は，実装容易かつ実用上早い手法です．手続きを述べると以下のようになります．

初期値 x[0] を与える
for k = 1, 2, ...
  t[k] := f(x[k]) / g(x[k])
  x[k+1] := argmin { f(x) - t[k] g(x) | x \in X }

収束性として次の結果が知られています．

• コンパクト集合上の連続関数に対して，このアルゴリズムは収束
• 線型などの多くのケースで超一次収束（cf. 二分探索は一次収束）

• 有限集合上の離散関数に対し，このアルゴリズムは収束
• {0,1}線型や劣モ/優モなどの多くの例で反復回数が強多項式（cf. 二分探索は最大関数値がかかる）

歴史

連続：特殊な場合については，Dinkelbach よりも前から同種の手続きが知られていました．例えば線型の場合は [Saaty-Gass 1954] や [Isbell-Marlow 1956]，二次の場合は [Ritter 1962] があります．[Jagannathan 1966] と [Dinkelbach 1967] は一般の連続関数について議論を行い，凸/凹 の場合について Dinkelbach が上のアルゴリズムを定式化しました．一般の収束証明は [Jeflea 2003] にあります（簡単なので，もっと前から知られていそうですが）．

離散：[Radzik 1992] が {0,1} 整数線型計画に対してこのアルゴリズムの計算量が強多項式で抑えられることを証明しました（それまでは parametric search が主流）．その後は具体的な問題に対するケーススタディが多いように見受けられます．また，ごく最近劣モジュラ / 優モジュラに対する同種の結果 [Kawahara-Nagano-Okamoto 2011] が出ていて，個人的には興味深いです．

他にも周辺はいろいろありますが，ありすぎて，あんまり簡潔にまとまっていないのが現状だと思います．「Recent Development と Fractional Programming」で検索して出てくるサーベイや，Handbook of Global Optimization を見るのがいいんでしょうか．

なお，こういった経緯があるのでアルゴリズム名も様々な名前で呼ばれています．自分の印象では連続系の人は Dinkelbach algorithm と呼ぶことが多く，離散系の人は Radzik の呼び名にしたがって Newton method や Discrete Newton method と呼ぶことが多いように思います．

実装例

ここでは離散分数計画における典型例である「最小比全域木問題」を実装します．これは無向グラフ G に2種類の非負重み a, b が与えられたとき w[a](T)/w[b](T) を最小化する全域木を求める問題です．

解法は上に述べたとおりで，パラメタ化された最小化問題 h(t) = min { w[a](T) - t w[b](T) | T \in Spanning Trees } を繰り返しとけば良いのですが，この問題は {0,1} 整数線型計画，すなわち w[a] - t w[b] = w[a - t b] という性質があるので，h(t) の評価は普通の最小全域木のルーチンで重みを a - t b に設定したものが使えます．

C++11 の勉強がてら，それっぽい実装ですが，まあ読めるとおもいます．コンパイルは g++ 4.7 で確認．

struct Graph {
  int n;
  struct Edge { int u, v, a, b; };
  vector<Edge> E;
  Graph(int n) : n(n) { }
  void addEdge(int u, int v, int a, int b) { E.push_back({u, v, a, b}); }

  double fcost, gcost;
  double minimumSpanningTree(double t) {
    fcost = gcost = 0;
    UnionFind unionFind(n);
    sort(E.begin(), E.end(), [t](Edge e, Edge f) { 
      return e.a - t*e.b < f.a - t*f.b; });
    for (Edge e : E) 
      if (unionFind.unite(e.u, e.v)) fcost += e.a, gcost += e.b;
    return fcost - t * gcost;
  }
  double fractionalMin() {
    for (double t = 0; ; t = fcost / gcost)
      if (fabs( minimumSpanningTree(t) ) < 1e-8) return t;
  }
};

参考文献

• T. Saaty and S. Gass (1954): "Parametric objective functions (Part 1)", Operations Research, vol. 2, pp. 316-319.
• J. R. Isbell and W. H. Marlow (1956): "Attrition games", Naval Research Logistic Quarterly, vol. 3, pp. 71-93.
• K. Ritter (1962): "Ein Verfahren zur Losung parameterabhangiger, nichtlinearer Maximum-Probleme", Unternehmensforschung, vol. 6, pp. 149-166.
• R. Jagannathan (1966): "On some propeties of progamming problems in parametric form pertaining to fractional programming", Management Science, vol. 12, pp. 609-615.
• W. Dinkelbach (1967): "On nonlinear fractional programming", Management Science, vol. 13, pp. 492-498.
• R. Cole (1987): "Slowing down sorting networks to obtain faster sorting algorithms", Journal of the ACM, vol. 37, no. 1, pp. 200-208.
• T. Radzik (1992): "Newton's Method for Fractional Combinatorial Optimization", Proceedings, 33rd Annual Symposium on Proceedings of Foundation of Computer Science, pp. 659-669.
• A. Jeflea (2003): "A parametric study for solving nonlinear fractional problems", Analele Stiintifice ale Universitatii Ovidus Constanta Seria: Mathematica, vol. 11, no. 2, pp. 72-92.
• R. van Oostrum and R. C. Veltkamp (2004): "Parametric search made practical", Computational Geometry: Theory and Applications - Special issue on the 18th annual symposium on computational geometry, vol. 28, no. 2-3, pp. 75-88.
• Y. Kawahara, K. Nagano, and Y. Okamoto (2011): "Submodular fractional programming for balanced clustering", Pattern Recognition Letters vol. 32, no. 2, pp. 235-243.

Shuffle Tree

10:59

Shuffle Treeは「確率的・自己組織化平衡二分探索木」です．Splay Tree, Scapegoat Tree, Treap の知識があると理解しやすいかもしれません．

テクニカルレポート [Patel 2009] が2009年に公開されていますが，仕事自体（＝提案）は2002年だそうです．なかなか面白いので是非論文誌などに出して欲しいのですが，経過を見るに，期待は薄そうです．

• [Patel 2009] M. Patel: A randomized self-adjusting binary search tree, technical report, http://www.clockandflame.com/media/ShuffleTree_Aug2009.pdf, 2009.

内容

二分探索木で挿入・削除・検索などを行う際には上から再帰的な探索を行いますが，Shuffle Treeではその際に「深い」ノードを見つけたら，その部分を「再平衡化」します（cf. Scapegoat Tree: rebuiliding）．

深いノードの条件は，基本的には「深さが log n を超えたら」ですが，乱択化しておきます．探索開始ノードで c = [0, n] の乱数を発生させ，降りるたびに c = (c-1)/n と更新， c = 0 以降のノードは「深い」と判定して，「再平衡化」の対象にします．※きちんと確認していませんが，多分期待値 O(n) の乱数を生成して log(n) で消化すれば同等の結果になるはずです．

ここで，c = 0 とか c = n にすれば乱択にしなくてもいいんじゃないかと思うかもしれませんが，それぞれ以下の問題があってダメで，乱択化することに真に意味があります．

• 毎回 c = n とすると浅いところに log n サイズのパス構造ができたとき解消されずに残る
• 毎回 c = 0 とすると劣化版 Splay Tree （二段回転なし）になるが．これはかなり遅い

再平衡化の手法は，木の回転操作に基づきます．あるノードから再帰的に左（or 右）に探索した後に，その方向のノードを上に持ち上げるように二分探索木の回転を行います．たとえば普通の2分探索木の挿入は

  insert(u, key) =
    if !u: return new Node(key)
    b = u.key < key
    u.ch[b] = insert(u.ch[b], key)
    return u

と実装できますが，Shuffle Treeでは，これを

  insert(u, key, c) =
    if !u: return new Node(key)
    b = u.key < key
    u.ch[b] = insert(u.ch[b], key, (c-1)/2)
    return rotate(u, b, c)

と書き換えます．ここで rotate(u, b, c) は c = 0 かつ回転可能のときにのみ u.ch[b] を上に持ち上げる木の回転操作です．

元の論文には「木を探索する一般的な関数 traverse を書いて，挿入・削除・検索全部の操作を traverse で実装するのだ」と書いてありますが，使いやすい traverse が思いつかなかったので，以下の実装では，上のように，それぞれを書き換えています（※実測したところ，検索は回転しないほうが早かったので，非回転にしています）．使いやすい traverse がかけた方は，教えてくれると嬉しいです．

パフォーマンスなど

• 実装量は Treap (insert-remove) より10行ほど長い
• 全体的に，挿入は Treap より少し遅い
• 全体的に，検索・削除は Treap より少し早い
• ソート例では Splay Tree 一回り遅い
• ランダム例では Splay Tree より一回り早い

実装量がTreapより長いのは，むしろTreapが異常です．Shuffle Treeの実装量は非平衡の二分探索木実装量＋回転実装量なので，これを下回るほうが凄いのです．

その他面白い点としては，ノードに全く平衡性情報を覚えさせなくて良い点があります（cf. 赤黒木：色情報, Treap：優先度）．これはSplay Treeの自己組織化性と同じ性質です．得意なデータも Splay Tree とおなじく「局所性のあるデータ（近いデータへのアクセス頻度が高い）」ですが，Shuffle Tree のほうが全体としてパフォーマンスのバラつきが少ない傾向にあります．

実装例

template <class T>
struct ShuffleTree {
  struct Node {
    T key;
    Node *ch[2];
    Node(T key) : key(key) { ch[0] = ch[1] = 0; }
  } *root;
  int size;
  ShuffleTree() : root(0), size(0) { }

  Node *rotate(Node *n, int b, int c) {
    Node *m = n->ch[b];
    if (c > 0 || !m) return n;
    n->ch[b] = m->ch[!b];
    m->ch[!b] = n;
    return m;
  }
  void insert(T key) { root = insert(root, key, rand() % (1 + size++)); }
  Node *insert(Node *n, T key, int c) {
    if (!n) return new Node(key);
    int b = (n->key < key);
    n->ch[b] = insert(n->ch[b], key, (c-1)/2);
    return rotate(n, b, c);
  }
  Node* removeMin(Node *n, T &value) {
    if (!n->ch[0]) { value = n->key; return n->ch[1]; }
    n->ch[0] = removeMin(n->ch[0], value);
    return n;
  }
  void remove(T key) { root = remove(root, key, rand() % (1 + size--)); }
  Node *remove(Node *n, T key, int c) {
    if (!n) return n;
    if (n->key == key) {
      if (!n->ch[1]) return n->ch[0];
      n->ch[1] = removeMin(n->ch[1], n->key); 
      return n;
    }
    int b = (n->key < key);
    n->ch[b] = remove(n->ch[b], key, (c-1)/2);
    return rotate(n, b, c);
  }
  Node *find(T key) { return find(root, key); }
  Node *find(Node *n, T key) {
    if (!n || n->key == key) return n;
    return find(n->ch[n->key < key], key);
  }
};

Double-Ended Priority Queue

11:43

普通の優先度付きキュー（ヒープ）は

• min
• push
• pop

の３つの操作が可能なデータ構造ですが，これに

• max

も可能になったデータ構造が Double-Ended Priority Queue です．

競技プログラミング的には std::multiset（平衡二分木）で代替するのが良いケースがほとんどだと思いますが，平衡木と完全二分木を配列で実装したヒープはキャッシュ効率で大きく差がつくので，計算量がシビアな問題ではこれが原因でTLEになることも十分考えられます（厳密には比較回数の下限にも定数倍で差があるのですが，実測ではそんなものよりキャッシュ効率が支配的です）．

DEPQの代表的な実装としては

• dual heap (aka. twin heap)
• min-max heap
• interval heap

の３つがあると思っています．以下これらを説明します．ちなみに，自分が一番好きなのはInterval Heapです．考え方がシンプルだし，この中では実測でも一番高速です．なお，最近でもいくつか新しいDEPQ実装は提案されていますが，それらの論文中に実測比較が無く，自分でも実装していないので，詳しいパフォーマンスの事情は知りません．

Dual Heap

一番基本的なDEPQの実装です．最小ヒープと最大ヒープを用意して以下のようにします．

• push(key): 両方のヒープに key を挿入
• deleteMin: minHeap を pop, 対応する要素を maxHeap から remove
• 他も同様

これは非常に単純なアイデアですが，字面の印象よりもずっと実装は大変です．試しに書いてみてください．

この手の同じデータ構造を並べてリンクでつなぐタイプの手法を Correspondence-based structure などといいますが，平衡二分木（Splay Tree）にすらパフォーマンス（比較回数と実行回数の両面）で勝てないと言われています（Chong-Sahni 2000）．

メリットはどんなデータ構造でも「とりあえず使える」ということで，例えばMeldable Double-Ended Priority Queueを作りたいけどアルゴリズムを知らない，といったケースではこの技法で実装するのが一つの戦略になります．

Min-Max Heap

Min-Max Heap（Atkinson-Sack-Santoro-Strothotte(1986)）は「偶数段目はMin-Heap，奇数段目はMax-Heap」であるヒープです．より正確には，以下の２条件を満たすものです：

• 偶数段目のノードの値は，その子孫よりも小さい（Min Heap property）
• 奇数段目のノードの値は，その子孫よりも大きい（Min Heap property）

push

要素を挿入する場合，普通と同じように完全二分木の末尾に挿入してheapUpするのですが，heapUpの前に一段回処理が必要です．自分が偶数段目（min level）にいる場合の処理は以下のようになります：

• 自分の親（max level）と比較して，自分のほうが小さい場合は min level の中で heapUp する
• 自分のほうが大きい場合は親と swap し，親から初めて max level の中で heapUp する

この処理でMin-Max Heap propertyが満たされることは，少し考えればわかります．

deleteMin

根と完全二分木の末尾を交換し，根をheapDownします．heapDownの処理は以下のようになります：

• 自分の子・孫の中の最小要素を m としたとき，自分が m より小さければ交換．
• 交換後の要素とその親でヒープ条件が崩れていたら修正．

deleteMax

根の子（深さ１のノード）のうち最大のものと末尾を交換し，heapDownします．

以下に実装例を示します．

// 51 lines
int log2(int n) { 
  int k = 0;
  for (; n; n >>= 1) ++k;
  return k-1;
}
const int MAXSIZE = 1000;
struct MinMaxHeap {
  int h[MAXSIZE];
  int size;
  void heapDown(int i, int b) {
    for (int m = i; ; i = m) {
      for (int c = i*2; c < min(size,i*2+2); ++c)
        if ((h[m] > h[c])^b) m = c;
      for (int g = i*4; g < min(size,i*4+4); ++g)
        if ((h[m] > h[g])^b) m = g;
      if (m == i) return;
      swap(h[m], h[i]);
      if (m - i*4 >= 0) 
        if ((h[m/2] < h[m])^b) swap(h[m/2], h[m]);
    }
  }
  void heapUp(int i, int b) {
    if (i > 1 && (h[i] > h[i/2])^b) { 
      swap(h[i], h[i/2]);
      b ^= 1; i /= 2;
    }
    for (int g; g = i >> 2; i = g) {
      if ((h[g] <= h[i])^b) break;
      swap(h[i], h[g]);
    }
  }
  void push(int x) { 
    h[size++] = x;
    heapUp(size-1,log2(size-1)%2); 
  }
  void deleteMin() { 
    h[1] = h[--size];
    heapDown(1, 0);
  }
  void deleteMax() {
    if (size == 2) return deleteMin();
    int m = 2;
    if (size >= 3 && h[3] > h[2]) m = 3;
    h[m] = h[--size];
    heapDown(m, 1);
  }
  int  findMin()   { return h[1]; }
  int  findMax()   { return size>3?max(h[2],h[3]):h[1+(size==3)]; }
  bool empty()     { return size == 1; }
  MinMaxHeap() : size(1) { }
};

Interval Heap

DEPQ実装で自分が一番好きなのがInterval Heap（Leeuwen-Wood, 1993）です．アイデアは非常に簡単で，以下のようなものです：

• 完全二分木の各ノードは「区間」に対応する（２つずつ値を格納する，末尾の要素はシングルトン [x,x] ）．
• [l,r] が [L,R] の子のとき [l,r] ⊆ [L,R]（Interval Heap Property）

最小値は根の左端点, 最大値は根の右端点を返せばOKです．

heapUp, heapDownも簡単で，変更のあった箇所と上下の包含関係をチェックして満たされていなければ満たされるように端点をスワップするのみです．

以下に実装例を示しますが，整理できていない段階で書いたものなので随分ヤヤコシくなっており，いずれ書き直す予定です．

// 54 lines
const int MAXSIZE = 10100100;
struct IntervalHeap {
  int h[MAXSIZE];
  int size, r;
  int &hof(int b, int i) { return h[b+2*i]; }
  void heapUp(int i, int b) {
    for (int p; p = i >> 1; i = p) {
      if (b ^ (hof(b,p) <= hof(b,i))) break;
      swap(hof(b,i), hof(b,p));
    }
  }
  void heapDown(int i, int b) {
    for (int c; (c = i << 1) < size+r; i = c) {
      if (c+1 < size+r && b ^ (hof(b,c) > hof(b,c+1))) ++c;
      if (b ^ (h[b+2*i] <= h[b+2*c])) break;
      swap(hof(b,i), hof(b,c));
      if (hof(0,c) > hof(1,c)) swap(hof(0,c), hof(1,c));
    }
  }
  void deleteMin() {
    if (r ^= 1) --size; else hof(1,size) = hof(0,size);
    hof(0,1) = hof(r,size);
    heapDown(1, 0);
  }
  void deleteMax() {
    if (size == r) return deleteMin();
    if (r ^= 1) --size; else hof(1,size) = hof(0,size);
    hof(1,1) = hof(r,size);
    heapDown(1, 1);
  }
  void push(int key) {
    hof(r,size) = key;
    if (r == 1) { // insert to large
      if (hof(0,size) > hof(1,size)) {
        swap(hof(0,size), hof(1,size));
        heapUp(size, 0); // min mode
      } else {
        heapUp(size, 1); // max mode
      }
    } else { // insert to large
      hof(1,size) = hof(0,size) = key;
      if (hof(0,size) < hof(0,size/2)) heapUp(size, 0);
      if (hof(1,size) > hof(1,size/2)) {
        heapUp(size, 1); 
        hof(0,size) = hof(1,size);
      }
    }
    if (!(r ^= 1)) ++size;
  }
  int findMin() { return hof(0,1); }
  int findMax() { return (size==r) ? hof(0,1) : hof(1,1); }
  bool empty() { return size+r == 1; }
  IntervalHeap() : size(1), r(0) { }
};

Pairing Heap

11:35

マージ可能ヒープの中で，現在「実用的に最速」と言われているものが Pairing heap です．

マージ可能ヒープなら，実装量的にはSkew Heapや乱択ヒープが非常に短いし，ヒープを取り替えたらTLEがACになるようなことはめったにないので，競技プログラミング的には微妙なネタですが，競技に直接使えるメジャーテクニックなんてだいたい語り尽くされているので，ここではそんなことは気にせずにやっていきます．

歴史的背景

Pairing Heap提案の2年前，Fredman-Tarjan(1984) がFibonacci Heapを提案しました．これにより最短路やらなんやらの計算量が全部落ちて凄い，という盛り上がりが起きたのですが，すぐに実用的には使えないと言われてしまいます（提案してすぐに「実用的には使えない」と言われるのはとんでもないことで，どれだけFibonacci Heapに注目が集まったかを表しています）．

一方，ほぼ同時期にSleator-Tarjan(1985) がLink Cut Treeの実装として，自己組織化平衡二分探索木のSplay Treeを提案しました．これは（当時知られていた）平衡木よりもずっとシンプルで，しかも実用的にも十分早いというものでした．

そこで，この２つのアイデアを組み合わせてFibonacci Heapを自己組織化版にする，という発想で作られたのが，Fredman-Sedgewick-Sleator-Tarjan(1986)のPairing Heapです．自己組織化アルゴリズムは実装がシンプルなわりに実用的に早いので，Fibonacci Heapの木の形状をamortizedで達成できれば実用的により高速になるだろうという発想です．

実装

コードを見ながら説明したほうがよいので，まずはコードを示します．Haskellにて．

data Heap a = Empty | Node a [Heap a] deriving (Show, Eq)

merge Empty h = h
merge h Empty = h
merge (Node x xs) (Node y ys)
  | x < y     = Node x (Node y ys : xs)
  | otherwise = Node y (Node x xs : ys)

mergeList []       = Empty
mergeList [x]      = x
mergeList (x:y:xs) = merge (merge x y) (mergeList xs)

push x h = merge (Node x []) h

pop (Node x xs) = mergeList xs

top (Node x xs) = x

ヒープのデータ構造は，各ノードがキーと子のリストをもつ可変長多分木で実装します．Binomial Heapの系列ですね．

マージ可能ヒープの基本は「mergeを作って，他の操作はmergeで実装」なので，Pairing Heapもそれに準じますが，木と木のmergeと，リストのmergeを区別します．木同士のmergeをmerge．リストのmergeをmergeListとします．

木同士のmergeは以下のように実装されます：

• キー同士を比較し，大きい方を小さい方の子リストに追加．

まあ，普通です．たとえば，要素を延々とpushしていくのは，線形リストの末尾に延々と要素を追加するのと同じ事になります．

続いて，mergeListは以下のようになります：

• リストの先頭から2つずつ要素を取り出し，それらをmergeする（リスト長1/2倍）．
• リストの後ろから順番にmergeして1つの木にする．

１段階目で２つずつ組化することからPairing Heapの名前がきています．１段階目と２段階目でmergeの方向が逆なのは少し重要で，この方向が一番計算時間が短くなります．理論保証するのは大変ですが，実測すると２倍程度の劣化があることがわかります．

C++での実装

上のアルゴリズムをC++に移植する場合に，２つの問題があります：

• 子リストを std::vector などで実装すると，あっさりメモリ溢れする
• mergeList を再帰実装すると，あっさりスタック溢れする

どちらも連続で insert された後のことを考えると，わかると思います．これらの問題が無ければPairing Heapは実装量的にもSkew Heapに匹敵する良データ構造になったのですが……．

どちらの問題も，解決自体は単純です．１つ目の問題は子をきちんとリストでもてばよいだけ．ただし，std::listが頭がオカシイレベルで遅いので，自分で単連結リストを作る必要があります．具体的には，ノードを以下のように定義します．

struct Node {
  int key;
  Node *head;
  Node *next;
};

headは子リストの先頭，nextは自分が含まれているリストにおける自分の次の要素を指します．

２つ目の問題の対応も簡単で，再帰をループに展開するだけです．実際の実装例を以下に示します．

struct PairingHeap {
  struct Node {
    int key;
    Node *head, *next;
    Node(int key) : key(key), head(0), next(0) { }
  } *r;

  Node *merge(Node *i, Node *j) {
    if (!i || !j) return i ? i : j;
    if (i->key > j->key) swap(i, j);
    j->next = i->head; i->head = j;
    return i;
  }
  Node *mergeList(Node *s) {
    Node n(0);
    while (s) {
      Node *a = s, *b = 0;
      s = s->next; a->next = 0;
      if (s) { b = s; s = s->next; b->next = 0; }
      a = merge(a, b); a->next = n.next;
      n.next = a;
    } // s == 0
    while (n.next) {
      Node *j = n.next;
      n.next = n.next->next;
      s = merge(j ,s);
    }
    return s;
  }
  int  top()         { return r->key; }
  void pop()         { r = mergeList(r->head); }
  void push(int key) { r = merge(r, new Node(key)); }
  bool empty()       { return !r; }
  PairingHeap() : r(0) { }
};

実行速度

普段ならここで幾つかの実装と実行速度比較をするところですが，今回はそれは各人のライブラリのマージ可能ヒープと比較してもらうことにします（自分で確認したほうが，わかりやすいと思いますし）．今回は既存の実装例・実行速度比較を紹介します．

C++ の boost に， boost::heap というプロジェクトがあります（http://www.boost.org/doc/libs/1_51_0/doc/html/heap.html）．ここにはいくつかのヒープの実装が置かれており，現在のラインナップは

• binomial heap
• d-ary heap
• fibonacci heap
• pairing heap
• skew heap

となっています．なぜかPairing Heapの実装が２の問題（再帰でスタック溢れ）を対処していない実装になっていますが，それでもそれなりに効率的には実装されています．

これらの実行時間（比較回数）比較が http://boost.2283326.n4.nabble.com/heap-A-performance-test-and-some-questions-tt3575491.html#a3581828 にあるので，見てみると良いかと思います．本当の見所は「効率的に実装されたFibonacci Heapの性能」という，世にも珍しいものだと思っていますが．

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だけメモしておいても十分事足りるかと思います．