Fibonacci ヒープを実装しました

参考にしたのは，いつもの CS166 のスライド．

「できる操作」「その操作でできるうれしいこと」「実装」の順に書きます（予定）．

英語が楽に読める人は上のスライドを見るとよさそうです．

できる操作

まず，これはヒープなので，優先度つきキューの基本三演算ができます．

• push(k, v)
  • 優先度が k である要素 v を追加する
  • \(O(1)\) time
• top()
  • 優先度が最も高い要素を返す
  • \(O(1)\) time
• pop()
  • 優先度が最も高い要素を取り除く
  • amortized \(O(\log n)\) time

便宜上，「要素 v を入れて v の値で比較する」ではなく「優先度 k を持つ要素 v を入れて k の値で比較する」ということにしておきます．

で，これに加えて以下の操作ができます．

• meld(h)
  • 別の Fibonacci ヒープ h と統合する
  • \(O(1)\) time
• prioritize(e, k')
  • // e は要素 (k, v) への参照
  • e の優先度を元々の k より高い k' に変更する
  • amortized \(O(1)\) time (!!)

優先度の操作がならし定数時間でできてうれしいです．

操作の呼び方について

prioritize は decrease-key と書かれたりすることが多いのかなぁ．min や max などに対して中立な表現を探していました．

prioritizeどうだろう

— 衛藤渚 (@eto_nagisa) October 26, 2019

ありがとうございました．

うれしいこと

例として，Dijkstra 法の計算量のオーダー改善を挙げます．

Dijkstra 法は単一始点の最短距離を求めるもので，概ね次のようなアルゴリズムです． 一回は見たことがある人を対象として書いているつもりです．

dijkstra(g, s) {
  // g はグラフを表す隣接リスト．頂点数 n で辺数 m とする
  // s からの最短距離を求める

  vector<weight_type> dp(n, inf);
  dp[s] = 0;
  priority_queue<pair<weight_type, index_type>> pq;
  pq.push({0, s});  // s への距離 が 0
  while (!pq.empty()) {
    weight_type w;
    index_type v;
    tie(w, v) = pq.top();  // 最も近い頂点 v とそれへの距離 w
    pq.pop();
    if (dp[v] < w) continue;  // w より短く行けるならパス
    for (auto e: g[v]) {
      weight_type = nw = dp[v] + e.cost
      index_type nv = e.target
      if (nw < dp[nv]) {
        // 短く行けることがわかったら更新してヒープにも入れる
        dp[nv] = nw
        pq.push({nw, nv});  // ***
      }
    }
  }
  return dp;
}

で，*** の部分なんですが，実際に要素 (nw, nv) を入れ直す必要はなくて，すでに入っている要素 nv の優先度を nw にすることができればよさそうですよね．

s 以外の要素も，優先度を inf として最初からヒープに入れておくと考えます．

そうすると，結局，ヒープに対する操作は次のようになります．

• \(n\) 回の push
  • \(O(n)\) time
• \(n\) 回の pop
  • 最初の追加の後，新たに頂点がヒープに入ることはないため
  • O(n log m) time
  • \(O(n \log n)\) time
• \(m\) 回の prioritize
  • 各頂点から出る辺に対して，一回ずつしか見られないため
  • \(O(m)\) time

また，グラフの性質から m = O(n²)で， O(log m) = O(log n) ですので， 全体で \(O(m + n \log n)\) time となります． そもそも \(n\) 個しかヒープ中に入らないので \(\log m\) とかは関係なかったですね．

実装

tl; dr：ポインタでたくさんがんばります．きらいな人はつらそう．

以下のような流れで説明します．

• ベースとなる二項ヒープについての話
  • binomial heap．二分ヒープ (binary heap) とは別
  • ポインタによる木に基づくヒープ
  • 上に書いた Fibonacci heap の操作のうち，prioritize 以外ができる
• prioritize を行う方法
• データ構造の表現方法
  • 愚直なポインタ木だと amortized \(O(1)\) time を実現できないので工夫する必要がある

二項ヒープ

2 べきの要素数のヒープたちを連結リスト（linked list）で持ちます． list<heap> のようなものをイメージするとよくて，これ全体で二項ヒープです．

連結リストで管理された各ヒープを「子ヒープ」，全体がなすヒープを「親ヒープ」と呼ぶことにします．

親ヒープが持っているのは次の二つです．

• 子ヒープを要素として持つ連結リスト
• 全体で最も優先度が高い要素を管理する子ヒープへのポインタ

子ヒープをどう表現するかの説明は後回しにしますが，子ヒープは以下の操作ができるとします．

• meld(h) in \(O(1)\) time
  • 別の子ヒープ h と統合する
  • これにより要素数が 2 倍になる
  • Note: \(2^i + 2^i = 2^{i+1}\)
• top() in \(O(1)\) time
  • 最も優先度の高い要素を返す
• pop() in \(O(i)\) time
  • // 対象の子ヒープの要素数を \(2^i\) とする
  • 最も優先度の高い要素を取り除く
  • 要素数 \(2^0\), \(2^1\), ..., \(2^{i-1}\) の子ヒープたちに分解する
  • Note: \(1 + 2^0 + 2^1 + \dots + 2^{i-1} = 2^i\)

さて，これのもとで，親ヒープに対する処理たちを考えます．

push

要素数 1 の子ヒープをつくり，リストの末尾に追加します．\(O(1)\) time でできますね． 必要に応じて，最も優先度が高い要素を持つ子ヒープへのポインタを書き換えます．

Note: \(1 = 2^0\)．

meld

リスト同士をくっつけるだけ．当然 \(O(1)\) time でできます．ポインタ書き換えも適宜．

top

最も優先度が高い要素を持つ子ヒープへのポインタを持っていますね． それが指すヒープに top を聞けばよいです．\(O(1)\) time．

pop

ぐわー．

まず，最も優先度が高い要素を持つ子ヒープをバラします． 上で述べた子ヒープの pop でできあがる子ヒープたちをリストに入れます．

それまでの push や，今の処理などの結果，要素数が少ない子ヒープがたくさんあるのがイメージできてもらえるとうれしいです．

この状態で次に優先度が高い要素を探そうとすると大変そうなので，できるだけ要素数を減らしておきたいです． 全体で要素数が \(n\) のとき，これらをまとめ上げて \(\log_2 n+O(1)\) 個の子ヒープにするのを目標にします．

要素数 \(2^i\) のヒープふたつを meld すると要素数 \(2^{i+1}\) の子ヒープがひとつできるのがポイントです．

連結リストで管理している子ヒープを順に見て，同じ要素数の子ヒープが現れたらそれらを meld していきます．

その結果，\(n\) を二進数で表したときに \(i\) bit 目が立っていれば \(2^i\) サイズの子ヒープができあがります．

\(n\) は \(\log_2 n+O(1)\) bit で表せるので，子ヒープの個数もこれで抑えられることになります．

あとは，このまとめられた子ヒープの中で最も優先度の高いものを持つものを探せばよいです．

計算量の直感的な話

\(n\) 回 push した後に pop を行うと \(O(n)\) time かかってしまいますが，続けてもう一度 pop した場合，子ヒープの個数が減っているので \(O(\log n)\) time で処理が終わります．

ならし計算量をお勉強しましょう．

子ヒープの表現

子ヒープは多分木で表します．

同じ要素数の木をくっつけるときは，片方の根を他方の根につなげます．優先度の高い方を新たな根にします．

要素数が大きくなる順に，次のような木になります．

子ヒープの形

二項ヒープでは，この木を left-child/right-sibling (LCRS) 表現で表します． 「\(i\) 番目の子は \(i\) 個目のポインタ」のようなものではなく，「子へのポインタ」と「次の兄弟へのポインタ」を持ち，二分木のように扱います．

これにより，たぶん根を pop するのを効率的に行えるんだったと思います．

Fibonacci ヒープの木では別の表現方法をするので，ここでは深入りしません．

ともあれ，ここまでの話をちゃんと読むと二項ヒープが実装できるはずです．

「子ヒープを論理的にはこういう木の形で持つ」というのと「その木をこういう方法で表現する」という二つのレイヤーがあることに注意するといいかもしれません． 上の図は前者に対応していて，LCRS 表現のリンクの張り方ではないです．

prioritize

さて，上で作った木はヒープ条件を満たしています．つまり，親はどの子よりも優先度が高いです．

あるノードの優先度を高めても，親の優先度を超えなければヒープ条件が崩れないのでそれで終了です．

親の優先度より高くなってしまった場合はどうしましょうか．

A Crazy Idea

そのノードを親から切り離してしまいましょう． そうすればヒープ条件が崩れませんね（おいおい）．

そのノードを根とする子ヒープをリストに追加しましょう．必要に応じてポインタも書き換えます．

これだけで終了です，だと話は単純なんですが，そうもいきません．

細かい子ヒープを簡単にたくさん作ることができてしまい，pop のならし計算量が悪化してしまいます．

対策

細かい子ヒープを簡単に作らせなければいいんですよね． より詳しくは，各子ヒープが次のことを満たしてくれるといいです*1．

• 根ノードが持つ子の数に対して，全体のノード数が指数的に大きい

ということで，根を除く各ノードについて，切り離せる子ノードはたかだか一つという制約をつけてみましょう．

一度子ノードが切り離されたとき，そのことを覚えておくようにします（damaged と呼んでいます）．

damaged なノード v から再び子を切り離したくなったときは，どうしましょうか． v と v の親も切り離すことにしましょう．これにより，連鎖的に切り離しが発生することはありますが，問題ないらしいので問題ないです．

というわけで，prioritize ができました．

木の表現方法

Fibonacci ヒープで扱う木では，以下のような操作が要求されます．

• ノード u の子にノード v を追加する
• ノード v の親にアクセスする
  • damaged のフラグを書き換え，および値の比較が必要
• ノード v とその親 u を切り離す

これらが定数時間で行えないと話になりません．

たとえば，子のポインタの配列を持ち，\(i\) 番目の要素は \(i\) 番目の子を指し，それとは別で親へのポインタも持ってみるとしましょう．次のような図に対応します．

木のよくない表現方法

子を切り離す際に，子の個数に対して線形とか対数時間かかってしまいますね．

そこで，次のような表現方法を採用します．

木のよい表現方法

兄弟間は circularly doubly-linked list になっています． 子から親へはリンクが張られていますが，親からは代表の子にのみリンクが張られています．

これにより，親からのリンクを張り替えるのは，代表の子を切り離すときのみでよくなります．

他の操作も定数時間で行えますね．

注意すること

要素数 1 の circularly doubly-linked list をつくるとき，隣の子を NULL とするか自分自身とするかはどちらでもいいと思いますが，要素数 2 の状態から雑に削除を行うと後者の状態になるので，意図しない状態になっていないか気をつけましょう．

説明おわり

にゃー．

ならし計算量などの議論をちゃんとする方がよかったのかもなんですが，割愛してしまいました． 冒頭にリンクを置いた元スライドには書かれているので，それを見るといいかもしれません．

割愛したせいでなぜ Fibonacci ヒープなのかの説明ができていませんね．解析の途中で Fibonacci 数列が出てきます．

実測

よく，フィボナッチヒープは遅いと言われますね．

実際に計測して確かめてみましょう．Dijkstra を貼ると解ける問題として，次の問題を使います．

atcoder.jp

グラフは \(3\le n \le 1002\), \(m = n(n-1)\) になっていると思います．

実行時間は次の通りです．

実行時間

まず，1268 ms のやつは Fibonacci ヒープで普通の Dijkstra を実装したものです．宝の持ち腐れですね，定数倍も重くてつらそうです．

114 ms のやつが STL の std::priority_queue です．二分ヒープを配列で表現するものが使われていて，やはり速いですね．

41 ms のやつは自前ヒープで，二分ヒープの配列表現に加え，\(O(\log n)\) time で prioritize をサポートしたものです．速いですね．

さて，37 ms のやつが Fibonacci ヒープによる \(O(m+n\log n)\) の Dijkstra です．一番速くて素敵ですね．うれしい．

上から，
- STL の二分ヒープの Dijkstra
- Fibonacci ヒープの愚直 Dijkstra, O(E log V)
- Fibonacci ヒープの Dijkstra, O(E+V log V)

オーダー改善しても定数倍で負ける流れだと思ってたのでうれしい pic.twitter.com/iiB4RZGKLU

— えびちゃん (@rsk0315_h4x) October 28, 2019

ついでに承認欲求も満たされてわーいわーいです．

提出コード

今回はこれでおわりです．おつかれさまでした．

*1:これでいい理由については，ならし計算量の議論をちゃんとする必要があります．