ZDDでマインスイーパーを解く話(卒研お勉強枠➂(最終回))

<< 前回

次回 >>(最終回のため、なし)

 

 

はじめに

今度こそ某有名パズルゲームを解く話です。

いよいよ明かされる某有名パズルゲームの正体とは…

 

ズバリ、マインスイーパです!!!

 

まあ、半年間ずっと𝕏(旧Twitter)でマインスイーパーソルバーの話をポスト(旧ツイート)してたので、何も隠せてないのですが…

 

この記事でやること

※本記事において、式番号をつけた等式は全てZDDの新たな演算として導入したものになります。詳細な説明は省きますので、詳しい実装を知りたい方はGitHubを参照してください。

 

マインスイーパーに関する諸定義

説明に使う用語

  • 盤面上の開示された数字をヒント数字と呼ぶ。
  • 爆弾がないと確定したマスにつける目印をと呼ぶ。
  • ヒント数字も旗もないマスを未開放マスと呼ぶ。
  • 未開放マスのうち、8方向の隣接したマスにヒント数字が1つ以上存在するマスを隣接マスと呼び、それ以外の未開放マスを非隣接マスと呼ぶ。

細かい設定

  • 爆弾総数の情報はプレイヤーに開示されているものとします。

 

マインスイーパーのZDDによる解法

ZDDによる爆弾配置の表現

隣接マスを台集合とし、条件を満たすような隣接マス内の各爆弾配置に対して、爆弾のあるマスの集合を組合せ集合とします。以下、特に断りのない限り台集合を  U、爆弾配置を列挙したZDDを  \mathcal{P} と表記します。

例えば、次の盤面(盤面1)では  \mathcal{P} = \{\{ a, d, e\}, \{c, d\}\} となります。さらに、この時点でマス  b に爆弾がないこととマス  d に爆弾があることが確定します。

盤面1 (爆弾総数3)

より一般的に書くと、 \mathcal{P}に対して

\begin{equation}\label{1} \mathcal{P}.\text{items}() := \bigcup_{P \in \mathcal{P}} P \tag{1} \end{equation}

\begin{equation}\label{2} \mathcal{P}.\text{flags}() := \bigcap_{P \in \mathcal{P}} P \tag{2} \end{equation}

とおくと

  • 爆弾がないと確定する隣接マスの集合:

\begin{equation} U \setminus \mathcal{P}.\text{items}() \end{equation}

  • 爆弾があると確定する隣接マスの集合:

\begin{equation} \mathcal{P}.\text{flags}() \end{equation}

と表せます。よって、 \mathcal{P} を作ることさえできればマインスイーパーを解き進めることができそうですね。

 

式 \ref{1}、\ref{2} はメモ化再帰により実装できます。

ここで、返り値の集合を表すデータ構造として前回の記事で導入した共有リストを利用すれば、メモリの消費量を少し抑えられると期待できます*1

 

新たな隣接マスの追加

マスを開くと、非隣接マスの一部が隣接マスに変化するため、それらを台集合に追加する操作が必要になります。

そこで、暫定のZDDとして、  \mathcal{P} の根に、新たなマスをラベルに持つ節点の0-枝と1-枝を突き刺します。こうすれば、「このマスには爆弾があってもなくてもよい」ということを暫定的に表現できますね。

新たな隣接マスをラベルに持つ節点の0-枝と1-枝を、根に突き刺す

ただし、ZDDでは節点がアイテムの添え字に対してトポロジカル順序を満たす必要があるため、このマスの添え字は(これまでに付した添え字の最小値)- 1 とします。

 

そして、この操作で得られた暫定のZDDに対し、新たなヒント数字と爆弾総数による条件を改めて適用することで、正しいZDDを構築することができます。

ここからは、それらに対応したZDDの新たな演算を2つ導入します。

 

ヒント数字を適用する演算

 U: 台集合, \mathcal{P} \subseteq 2^U, E \subseteq U, h \in \mathbb{N}_{\ge 0} に対して

\begin{equation} \mathcal{P}.\text{hint}(E, h) := \{P \in \mathcal{P} \mathrel{|} |P \cap E| = h\} \tag{3} \end{equation}

と定義します。

ここで、 h をヒント数字、 E をヒント数字が影響する隣接マスの集合とすれば、暫定のZDD  \mathcal{P}からヒント数字の条件を満たす組合せ集合を抜き出す演算になりますね。

 

実はこの演算、メモ化再帰で処理出来ちゃいます。びっくり。

引数の  E はまさに共有リストの出番ですね。

 

爆弾総数を適用する演算

 U:台集合, \mathcal{P} \in 2^U, b \in \mathbb{N}_{\ge 0}, m \in \mathbb{N}_{\ge 0} に対して

\begin{equation} \label{4} \mathcal{P}.\text{limit}(b, m) := \{P \in \mathcal{P} \mathrel{|}  0 \le b-|P| \le m\} \tag{4} \end{equation}

と定義します。

ここで、 b を爆弾総数、 m を非隣接マスの個数とすれば、暫定のZDD  \mathcal{P}から爆弾総数の条件を満たす組合せ集合を抜き出す演算になりますね。

この演算もメモ化再帰で処理出来ちゃいます。やったぜ。

 

ここまでに見た2つの演算hint, limitを適用することで、暫定のZDDから正しいZDDを得ることができます。

 

爆弾総数のもう一つの使い方

次の盤面(盤面2)で、爆弾がないと確定するマスはどこでしょう?

盤面2(爆弾総数3)

(しばらくシンキングタイム用の空白)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

…少し考えてもらえればわかると思いますが、右端の列の3マスには爆弾がないと確定します。

左と中の2列で爆弾を3つ使い切るため、右列に爆弾はない

より一般的には、爆弾総数を  b、非隣接マスの個数を m としたとき、非隣接マスに爆弾がないと確定することと以下は同値です:

\begin{equation} \forall P \in \mathcal{P}; b-|P| = 0 \end{equation}

この判定は式 \ref{4} のソースコードをちょちょいと弄ればできます*2

 

危険度の計算

爆弾がないと確定するマスがない場合は、危険度(そのマスに爆弾があるような爆弾配置の総数)が最も低いマスを開くことにします。

暇な人は、次の盤面(盤面3)で実際に数えてみるとよいでしょう。シンプルな盤面なのでまだマシですが、それでも人力でやるにはうんざりする作業であると実感できるかと思います。

盤面3(爆弾総数2)a: 1, b: 3, c: 3, d: 1, 非隣接マス: 2

 

この計算を、ZDDを使ってうまいことやります。

 

隣接マスの危険度

 \text{dp}_{Total}[e][k] :=  \mathcal{P} のうち、隣接マス  e が含まれるサイズが  k の組合せ集合の個数

と定義します。組合せ集合はZDDの根から1-葉までのパスに対応しているので、競プロerならお馴染みのDAG上のトップダウンボトムアップ2方向のDPで求められますね。

サイズが  k の組合せ集合(= 爆弾が  k 個あるような隣接マス内の爆弾配置)に対しては、非隣接マス内の爆弾配置の総数  \binom{m}{b-k} b: 爆弾総数, m: 非隣接マスの個数)による重みがかかるので、非隣接マス  e の危険度は次で求められます:

\begin{equation} \sum^{b}_{k=1} \text{dp}_{Total}[e][k] \times \binom{m}{b-k} \end{equation}

非隣接マスの危険度

 \text{dp}[k] :=  \mathcal{P} のうち、サイズが  k の組合せ集合の個数*3

とします。どの非隣接マスについても、危険度の計算式は次の通りです:

\begin{equation} \sum^{b-1}_{k=0}  \text{dp}[k] \times  \binom{m - 1}{b - k - 1} \end{equation}

 

旗を立てる

爆弾があると確定したマスについては、盤面上に旗を立てて爆弾総数を1減らすことで対応します。また、ヒント数字は周囲の旗の個数分値を減らして処理するようにします。旗を立てたことで隣接マスの定義から外れるため台集合から取り除く必要がありますが、これはonset演算で処理できます。

この操作は行わなくても理論上は解くことができますが、持っているZDDの節点数を大きく減らせるため、解答速度がかなり速くなります。

 

不要節点の削除

ZDDや共有リストはハッシュテーブル上に永続的に節点を保存するため、マインスイーパーを解き進めていくと、爆弾の有無が確定したマスの節点やその祖先の節点がどんどん増えていきます。これらは今後共有されることのない不要な節点なので、削除してしまって問題ないはずです

そこで、逆辺を追加し、ラベルごとに節点を管理することで、開くor旗を立てると同時に、節点の削除を行えるようにしました*4

 

不要節点の削除は、盤面が大きければ大きいほど効果が発揮されます。

200×200, 爆弾割合10% での、実行時間に対するメモリ消費量の推移

全体フローチャート



 

実験結果

盤面サイズと爆弾総数割合ごとに50戦し、勝数と解答時間を集計しました。

ソースコード

  • 実装言語: Python 3.11.4 (tags/v3.11.4:d2340ef, Jun  7 2023, 05:45:37) [MSC v.1934 64 bit (AMD64)] on win32
  • OS: Windows11
  • CPU: Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz   2.30 GHz
  • メモリ: 8GB

勝数

 

20×20

40×40

60×60

80×80

10%

47

47

49

46

20%

36

17

13

8

解答時間

「勝った試合のみ」と「全試合」それぞれで解答時間の最小値、中央値、最大値を記載。(単位: sec)

 

 

20×20

40×40

60×60

80×80

 

 

10%

最小

0.012

0.012

0.064

0.064

0.156

0.156

0.351

0.351

中央

0.016

0.016

0.076

0.076

0.191

0.191

0.424

0.424

最大

0.042

0.042

0.110

0.110

0.290

0.290

0.757

0.757

20%

最小

0.020

0.002

0.139

0.004

0.858

0.005

3.357

0.006

中央

0.038

0.035

0.246

0.243

1.396

1.400

7.947

5.275

最大

0.128

0.128

0.707

1.121

17.471

23.505

38.111

74.303

ちなみに:『Microsoft Minesweeper』の「超上級」は 16×30, 爆弾総数99 なので、この表での 20×20, 20% に相当。

 

余談

恐らく私の大学最後のブログ執筆になるので、学部4年の思い出の記録として研究の背景をまったりと書きました。

研究テーマを決めたきっかけの話

夏休みに入る直前くらいのこと。友人とパズルゲームの話をしていて、自分もなんか買ってみるかーとSteamを漁ってみたらこんなの↓を見つけました。

store.steampowered.com

マインスイーパーは昔から好きでよく遊んでいたので、面白そうだと思い即購入。

プレイしてみると、Windows標準搭載のマインスイーパーとは決定的に違う点がありました。

なんと、全ての問題が論理的に完全解答できるように作られており、開示された情報からは爆弾の有無が特定できないマスに対して操作(開く or 旗を立てる)をすると、自動的にそれを満たさない解を生成し、その時点でゲームを終了させる機能が実装されていたのです。

プレイ画面

未確定のマスをクリックするとゲームが終了する

「なぜそうなる?」をクリックすると、実際に爆弾配置を見ることができる

これにそこそこ感動し、おそらく内蔵されているのであろうマインスイーパーのソルバーを、自分でも作ってみたいと思いました。

 

ZDDについて勉強していると、パズルを解くのに使えるという情報をいくつか目にします。そこで、これマインスイーパーも解けるんじゃね…?と軽く考察を進めてみると、なんか解けちゃいました。

 

というわけで、夏休みの間ひたすら実装して、夏休み明けの初回ゼミで成果を発表したら、教授っちから学会で発表しようと提案されたという流れです。貴重な経験でした。

 

ちなみに 参考文献[1] はソルバーを作ってから見つけました。(アホ)

良い子のみんなは先に先行研究を調べてから卒研に着手しましょう。

 

参考文献

  1. 鈴木 浩史, 孫 浩, 湊 真一. BDD/ZDDを用いたマインスイーパーの爆弾配置パタンの列挙. The 30th Annual Conference of the Japanese Society for Artificial Intelligence, 1D5-OS-02b-3in2, 2016.
  2. 湊 真一. Zero-suppressed BDDs for Set Manipulation in Combinational Problems. 30th ACM/IEEE Design Automation Conference,1993.
  3. 3. 湊 真一. Zero-suppressed BDDs and their applications. Int. J. Softw. Tools Technol. Transf. (STTT), vol. 3, no. 2, pp. 156-170, Springer 2001

 

<< 前回

次回 >>(最終回のため、なし)

*1:これはかなり例外的な使い方で、最悪計算量の変化は(おそらく)ありません。要検証です。

*2:私の実装では演算と判定を同時に行っています。

*3: \text{dp}_{Total}の計算過程で行うボトムアップ方向のDPに関して、根に残った結果にあたります。

*4:より厳密な話をすると、ガベージコレクション君に仕事をしてもらうため、参照元を全て断ち切る必要があります。つまり、指定した節点とその祖先だけでなく、メモされた再帰関数のポインタなど色んなものを消さないといけません。実装超面倒くさかった…。

ZDDによるN-Queen problemの解法を効率化する話(卒研お勉強枠②)

<< 前回

次回 >>

 

 

はじめに

すっごく久しぶりの投稿です。本当はもっと早く投稿するつもりだったのですが、教授っちの指示で卒研お勉強枠】から【卒研】を通り越して【学会発表】に飛躍したので、投稿を自粛してました。びっくり。

 

さて、某有名パズルゲームをZDDで解くと前回予告していたのですが、卒研やってるうちに謎に成果が発生してしまったので、2部に分けることになりました。まずは前半として、前回扱ったN-Queen ProblemをZDDでさらに効率的に解く方法を紹介します。

しれっと言ってますが、これがとんでもなく応用の利きそうなアイディアとなります。

 

この記事でやること

  • 補助データ構造: 共有リスト の導入
  • ZDDを用いたN-Queen Problemの解法の効率化

ソースコードGitHubで公開してます。

 

前回のおさらいと効率化アイディア

N-Queen Problem (N=8)

N-Queen Problemは、N×Nのチェス盤にN個のクイーンを相互に攻撃が当たらないように置く方法は何通りか?という組合せ問題です。

詳細な解法は前回の記事にて説明しているので省略しますが、上の行から順に演算を組み合わせてZDDを構築していくことで解くことができます。

上の行から順にZDDを構築する ( \mathcal{P}: 暫定のZDD)

この解法では、クイーンの置き場所を決め打った際に何度もoffset演算を呼び出すことで☒印のマスを含まない組合せ集合を抜き出しています。offset演算の具体的な手続きを思い出すと、offsetしたいマスよりも添え字が小さいマスを表す節点を何度も生成し、それらの節点に対し次のoffsetで再び再帰関数が呼び出されることになります。これ、かなり無駄な動きをしている気がしませんか?

実は、offsetの引数を集合に拡張した演算offset_sを考えると、次のようなメモ化再帰で手続きを記述できます。

    @cache
    def offset_s(self, V:frozenset):
        """
            selfにおいて、各v∊Vを含まない組合せ集合全体をなす部分zddを得る。
            例:
                self := {ac, ab, ad, c}
                のとき
                self.offset({b, c}) = {ad}
        """
        if len(V) == 0 or self is self.zdd.zero_terminal or self is self.zdd.one_terminal:
            return self
        v = min(V, key=lambda v: self.zdd.idx[v])
        if self.top == v:
            return self.zero.offset_s(V-frozenset([v]))
        elif self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self.offset_s(V-frozenset([v]))
        else:
            return self.zdd._get_node(self.top, self.zero.offset_s(V), self.one.offset_s(V))

メモ化する関係上、引数の集合を表すデータ構造としてPythonのfrozensetを用いています。これでは結局、各関数内で時間・空間両方で集合のサイズ分の計算量をかけてしまうため、あまり解決にはなりません*1。この問題を解決する都合の良いデータ構造があればいいのですが…。

 

共有リスト

前述の再帰関数offset_sを実現するうえで、引数の集合を表すデータ構造が満たすべき要件は以下の通りです:

  • 一致判定(メモ化するため)
  • 添え字が最小のアイテムを除いた新たな集合を生成(再帰呼び出し時に使用)

これらが \mathcal{O}(1)でできれば革命がおきますね。

 

革命がおきました。

実は意外と簡単で、ZDDから0-枝と0-葉を省くことで得られる単方向連結リストを考えれば、上記要件を満たします。これを、共有リストと呼ぶことにしましょう。

共有リスト

基本的な実装は、ほとんどZDDと変わりません。

ZDD同様ハッシュテーブル上に節点を保存し、節点生成時にget_node関数を挟むことで共有規則を適用します。

 

また詳細は省きますが、共有リストは和集合や共通部分を取るといった基本的な集合演算もZDDと同様の再帰関数で記述できます。(ただし、メモ化は必要ありません。)

 

その他の効率化ポイント

マスの添え字を逆順に振る

マスの添え字を逆順に振る

図の通りです。

実は、前回N=11までしか解けなかったのが、これだけでN=12も解けるようになりました。

 

左右反転解を省略する

Nが偶数の場合は、1行目のみ左半分にしかクイーンを置けないとすることで左右反転解を省略できます。

Nが奇数の場合は、1行目で真ん中に置き、2行目で左半分に置くパターンを考える必要があります。

 

実験結果

ソースコード

  • 実装言語: Java SE 11*2
  • OS: Windows11
  • CPU: Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz   2.30 GHz
  • メモリ: 8GB

     

   

offset

offset_s

N

節点数

解答時間(ms)

節点数

解答時間(ms)

8

92

4,596

9

1,091

4

9

352

18,069

89

4,091

27

10

724

67,534

719

14,735

53

11

2,680

294,397

10,471

61,830

302

12

14,200

1,265,628

312,421

272,413

1,419

13

73,712

-

-

1,343,397

7,650

14

365,596

-

-

6,715,477

71,288

速い!

 

おわりに

共有リストを導入したことで、集合を引数に取るZDDの演算を効率よくメモ化できるようになりました。今回はoffset演算を拡張しましたが、前回の記事紹介したonset演算やchange演算も同様に拡張できます。

次回は某有名パズルゲームを扱いますが、そこで新たな集合を引数に取る演算が登場します。共有リスト、結構便利。

 

参考文献

  1. 湊 真一. Zero-suppressed BDDs for Set Manipulation in Combinational Problems. 30th ACM/IEEE Design Automation Conference,1993.
  2. 湊 真一. Zero-suppressed BDDs and their applications. Int. J. Softw. Tools Technol. Transf. (STTT), vol. 3, no. 2, pp. 156-170, Springer 2001

 

<< 前回

次回 >>

*1:実際のところは、ある程度改善されるようです。

*2:いつもはPythonで書いているのですが、この実験に関してはメモリがカツカツでPythonだと実行時間のブレが大きかったので、Javaでの実験結果を載せています。

ZDDと8-Queen problemの話【卒研お勉強枠①】

 

次回 >>

 

 

はじめに

競プロをのほほんとやりながら大学生活を過ごしてきましたが、実は大学を卒業するには卒論とかいうものを書かないといけないらしいです。困りました…。

一応ざっくりとやりたいことは決まってて、それをやる武器としてどうやら「ZDD」とかいうデータ構造が使えそう(な気がする)(使ってみたらどうなるかな?くらいの気持ち)なので、しばらくそれのお勉強をします。

せっかくなので、お勉強したことをブログにまとめていきたいと思います。やる気が続けば。

 

この記事でやること

・ZDDの解説と、各種演算のPython3.11による実装

・ZDDを用いて8-Queen problemを解く。

 

ZDD(Zero-suppressed Binary Decision Diagram, ゼロサプレス型二分決定グラフ)

ZDDとは?

ZDDは、与えられた集合(台集合)に対しその部分集合族を表現するデータ構造です。

特に、台集合に対しある種の条件を満たす部分集合のことを組合せ集合と呼び、ZDDはこの組合せ集合全体を表現するのによく用いられます。例えば「ナップサック問題の重量制限に収まる商品の集合全体」とか、「グラフ上のs-tパス全体」とかですね。

どうやって表現するの?

(各画像は [3] より引用)

台集合を  V := \{v_0, \dots, v_{n-1}\}、組合せ集合全体を  \mathcal{C} とします。また、組合せ集合を書くときは、含まれる要素を並べて略記(例:  \{v_0, v_1, v_2\} \rightarrow v_{0}v_{1}v_{2})し、特に空の組合せ集合を  \lambda で表します。

まず、以下の手順で Binary Decision Tree (二分決定木) を構築します。(長いので、クソデカ木と呼称することにします):

  • 高さ  n の完全二分外向木(根→葉方向の有向木)を用意し、深さ  i \in \{0, \dots, n-1\} の各ノードには  v_i をラベル付けしておく。また、各ノードから出る2本の枝をそれぞれ 0-枝1-枝 と呼ぶ。
  • 木なので、葉を1つ選ぶと根からその葉までのパスが一意に定まる。このとき「 v_i のノードから出るときに1-枝を通る」⇔「 v_i を含む」とすることで  V の部分集合に対応付けし、各葉に「対応する部分集合が  \mathcal{C} に属するかどうか」のbool値を記録する。

クソデカ木。台集合の要素が  v_i ではなく  e_i になっているけど許して。

クソデカ木は明らかに  \mathcal{C} を一意に表現できていますが、  \mathcal{C} の実態に依らず  2^n 程度のノードを必要となり効率的ではありません。発想は悪くないのですが、実用上  \mathcal{C} 2^V に比べて疎であることが多いので、実態に応じたノード数で  \mathcal{C} を表現できれば嬉しいです。

ZDDは、クソデカ木に対し以下の2種の操作を可能な限り繰り返すことでノードを削減したものとなります*1

  1. 冗長接点の削除
    • 1-枝が0の値を持つ葉を指している場合に、この節点を取り除き、0-枝の行き先に直結させる。
  2. 等価接点の共有
    • 等価な節点(ラベルが同じで、0-枝同士,1-枝同士の行き先が同じ)を共有する。
    • ※特に、葉は False/True それぞれ1つにまとめ、それぞれ0-葉 (あるいは \bot), 1-葉 (あるいは  \top)と表します 。

ノード削減のイメージ図

実際にクソデカ木からZDDを構築する様子

ただ、この手順でZDDを構築しようとするのは、最初にクソデカ木を陽に作ることになるので実質不可能です。悲しいね。

ではZDDをどう扱うかというと、部分集合族に対する様々な演算(和集合をとる など)をZDD上での操作で表現するというのが基本的な用法になります。

 

ここからは、具体的な実装を見ながらZDDの各種演算を解説していきます。

 

ZDDのプログラム上での表現方法と実装準備

まず実装上の前提として、ZDDクラスDiagramクラスを分けて作ります。

ZDDクラスには基礎的な組合せ集合のDiagramを作る関数を用意し、Diagramクラスにメインとなる演算を書いていきます。

Diagramと言っても、インスタンスそのものはノードを表現します。このクラス名は、各ノードを根とする部分グラフによって小さな部分集合族が表現されていることを踏まえています。

また、ZDDインスタンス上のハッシュテーブルにDiagramインスタンスを記録していく(正確にはDiagramインスタンス生成する関数 _get_node() をメモ化する)ことで、全体として永続データ構造のように扱います。

ZDDクラスのコンストラク
  •  \text{self.Diagram} は扱うDiagramクラスを表します。*2
  •  \text{self.zero_terminal} が0-葉、 \text{self.one_terminal} が1-葉を表し、これらは特別に属性として記録しておきます。
  •  \text{self.idx} は 台集合の各要素の添え字を表し、葉のラベル(None)には最大の添え字を番兵として付しておきます。
class ZDD:

    def __init__(self, ground_set, Diagram=Obsolete_Diagram) -> None:
        self.idx = {v: i for i, v in enumerate(ground_set)}
        self.idx[None] = len(ground_set)
        self.Diagram = Diagram
        self.zero_terminal = self.Diagram(zdd=self)
        self.one_terminal = self.Diagram(zdd=self)
Diagramクラスのコンストラクタ、ハッシュ関数
  •  \text{self.zdd} は自身が属するZDDインスタンスを表します。
  •  \text{self.zero}, \text{self.one} はそれぞれ 0-枝、1-枝 の先にあるDiagramインスタンスを表します。
  •  \text{self.top} はノードのラベルを表します。
  • ハッシュ関数は、2つの葉については乱数で、そうでない場合はラベルと左右のDiagramインスタンスハッシュ値から計算します。これは、特に【等価節点の共有】に用います。
  • Diagramのインスタンスに対し破壊的な処理を行うことはないので、コンストラクタ内で計算して  \text{self.hash} に記録してそれを  \text{self.__hash__()} で呼ぶようにしています。
  • ※ コードではObsolete_Diagramクラスとしています。
class Obsolete_Diagram:

    def __init__(self, zdd, top=None, zero=None, one=None) -> None:
        self.zdd = zdd
        self.top = top
        self.zero = zero
        self.one = one
        if top is None:
            self.hash = randint(0, 10**9)
        else:
            self.hash = (top, zero.hash, one.hash).__hash__()
ZDDを図示する関数 Diagram.show(self, filename)
  • Diagramクラス で宣言します。
  • pyvisを別途インストールし、from pyvis.network import Network する必要があります。
  • filename は .html (HTML拡張子) である必要があります。
    def show(self, filename):
        """
            selfによるzddを図示
            filemame: html拡張子のファイル名
        """
        nodes = {self: 0}
        net = Network(directed=True)
        dq = deque()
        if self is self.zdd.zero_terminal:
            net.add_node(nodes[self], shape="box", label="⏊")
        elif self is self.zdd.one_terminal:
            net.add_node(nodes[self], shape="box", label="⏉")
        else:
            net.add_node(nodes[self], shape="circle", label=self.top.__repr__())
            dq.appendleft(self)
        while dq:
            p = dq.pop()
            if p.zero not in nodes:
                nodes[p.zero] = len(nodes)
                if p.zero is self.zdd.zero_terminal:
                    net.add_node(nodes[p.zero], shape="box", label="⏊")
                elif p.zero is self.zdd.one_terminal:
                    net.add_node(nodes[p.zero], shape="box", label="⏉")
                else:
                    net.add_node(nodes[p.zero], shape="circle", label=p.zero.top.__repr__())
                    dq.appendleft(p.zero)
            net.add_edge(nodes[p], nodes[p.zero], color="red")
            if p.one not in nodes:
                nodes[p.one] = len(nodes)
                if p.one is self.zdd.zero_terminal:
                    net.add_node(nodes[p.one], shape="box", label="⏊")
                elif p.one is self.zdd.one_terminal:
                    net.add_node(nodes[p.one], shape="box", label="⏉")
                else:
                    net.add_node(nodes[p.one], shape="circle", label=p.one.top.__repr__())
                    dq.appendleft(p.one)
            net.add_edge(nodes[p], nodes[p.one], color="blue")
        net.show(filename, notebook=False)
  • 適当なインスタンスで実行すると、例えば下図のように表示されます。

 \{abcd, acd\} のZDD
  • 図の見方は以下の通りです:
    • 赤辺:0-枝
    • 青辺:1-枝
    •  \bot:0-葉
    •  \top:1-葉
インスタンスごとに関数をメモ化するデコレーター
  • グローバル に宣言します。
  • Pythonでは functools.cache を使えば関数をメモ化できますが、インスタンスごとにテーブルを分けずに記録してしまうため、インスタンスをdelする際のメモリの解放が困難になります。(する方法があれば教えていただけると幸いです!)
  • そこで、インスタンスごとにテーブルを作ってそこにメモするようなデコレーターを作っておきます。使い方は functools.cache と全く同じです。
def cache(func):
    """
        インスタンスメソッドをメモ化するデコレーター
    """
    def _wrapper(self, *args):
        if not hasattr(self, "_cache_table"):
            self._cache_table = {}
        key = (func, args)
        if key not in self._cache_table:
            self._cache_table[key] = func(self, *args)
        return self._cache_table[key]
    return _wrapper
Diagramインスタンスの生成に用いる関数 ZDD._get_node(self, v, zero, one)
  • ZDDクラス で宣言します。
  • Diagramインスタンス生成時に使い、2つの削除規則を漏れなく適用します。
    1. 【冗長接点の削除】は if 分岐で適用されます。
    2. 【等価節点の共有】は @cache によるメモ化で適用されます。
  •  v, zero, one は生成するDiagramのラベル、0-枝の先、1-枝の先 を表します。
    @cache
    def _get_node(self, v, zero, one):
        if one is self.zero_terminal:
            return zero
        return self.Diagram(self, v, zero, one)
アルゴリズムの説明用の記号・記法
  • Diagramインスタンス  \mathcal{P} は、部分集合族を表す記号としてそのまま利用する。
  • 二項演算  \triangleleft v \in V とどの要素にも  v が含まれないような部分集合族  \mathcal{P} \subset 2^V の間に定義し、

    \begin{align} \mathcal{P} \triangleleft v := \{C \cup \{v\} \mathrel{|} C \in \mathcal{P} \} \end{align}

    と定める。
  • 定義より

\begin{align} \mathcal{P}\text{.zero} & = \{C \mathrel{|} C \in \mathcal{P}, \mathcal{P}\text{.top} \notin C\} \\ \mathcal{P}\text{.one} & = \{C \setminus \{\mathcal{P}\text{.top}\} \mathrel{|} C \in \mathcal{P}, \mathcal{P}\text{.top} \in C\} \end{align}

が成り立ちます。

  • 特に

\begin{align} \mathcal{P} = \mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) \end{align}

が成り立ちます。逆に言えば、Diagramインスタンス  \mathcal{X}, \mathcal{Y} と添え字が  \mathcal{X}\text{.top}, \mathcal{Y}\text{.top} より小さい  v \in V に対し

\begin{align} \mathcal{Z} = \mathcal{X} \cup (\mathcal{Y} \triangleleft v) \end{align}

が成り立つならば、 \mathcal{Z}

\begin{cases} \mathcal{Z}\text{.top} & := & v \\ \mathcal{Z}\text{.zero} & := & \mathcal{X} \\ \mathcal{Z}\text{.one} & := & \mathcal{Y} \end{cases} とするDiagramインスタンスとして表現できます。

 

準備が完了したので、ここから各種演算の実装に入っていきます。

 

ZDDの各種演算

※ 説明文中の小文字zddは、ZDDクラスのことではなく、部分集合族を表現する本来の意味でのZDDを表します。

 \mathcal{P}, \mathcal{Q} はDiagramインスタンスとします。

 

実装する演算は以下になります:(具体例はコードにコメント文で入れてあります。)

  1.  \text{ZDD.empty_set}(\text{self})
    • 空のzdd( \{\})を得る。
  2.  \text{ZDD.unit_set}(\text{self})
    • 空の組合せ集合のみを持つzdd( \{ \lambda \})を得る。
  3.  \text{ZDD.single_literal_set}(\text{self}, v)
    • 組合せ集合  v のみを持つzdd( \{v\})を得る。
  4.  \text{Diagram.offset}(\text{self}, v)
    •  \text{self} において、 v を含まない組合せ集合全体をなす部分zddを得る。
  5.  \text{Diagram.onset}(\text{self}, v)
    •  \text{self} において、 v を含む各組合せ集合から  v を除いて取り出したzddを得る。
  6.  \text{Diagram.change}(\text{self}, v)
    •  \text{self} の各組合せ集合について、 v が含まれて いる/いない を反転させたzddを得る。
  7.  \mathcal{P} |  \mathcal{Q}
    • 部分集合族  \mathcal{P}, \mathcal{Q} に対し、 \mathcal{P} \cup \mathcal{Q} を表すzddを得る。
  8.  \mathcal{P} &  \mathcal{Q}
    • 部分集合族  \mathcal{P}, \mathcal{Q} に対し、 \mathcal{P} \cap \mathcal{Q} を表すzddを得る。
  9.  \mathcal{P} + \mathcal{Q}
    • 部分集合族  \mathcal{P}, \mathcal{Q} に対し、 \mathcal{P} \cup \mathcal{Q} を表すzddを得る。( \mathcal{P} |  \mathcal{Q} と等価)
  10.  \mathcal{P} - \mathcal{Q}
    • 部分集合族  \mathcal{P}, \mathcal{Q} に対し、 \mathcal{P} \setminus \mathcal{Q} を表すzddを得る。
  11.  \mathcal{P} * \mathcal{Q}
    •  \mathcal{P} \mathcal{Q} から組合せ集合を1つずつ取り、それらの和集合をとって得られる組合せ集合」全体を表すzddを得る。
  12.  \mathcal{P} // \mathcal{Q}
    •  \mathcal{Q} \neq \{\} とする。
    •  \mathcal{Q} の各要素  C に対し「 \mathcal{P} C を含む各組合せ集合から  C を除いて取り出した集合全体」を取ったときの、それらの共通部分を表すzddを得る。
  13.  \mathcal{P} \% \mathcal{Q}
    •  \mathcal{Q} \neq \{\} とする。
    •  \mathcal{P} - \mathcal{Q} * (\mathcal{P} // \mathcal{Q}) を表すzddを得る。
1. ZDD.empty_set(self)
  • ZDDクラス で宣言します。
  • 0-葉がそのまま求めたいzddを表します。
    def empty_set(self):
        """
            空のzdd ({}) を得る。
        """
        return self.zero_terminal
2. ZDD.unit_set(self)
  • ZDDクラス で宣言します。
  • 1-葉がそのまま求めたいzddを表します。
  • unitというのは、11. で定義される積の単位元をなすことから名づけられています。
    def unit_set(self):
        """
            空の組合せ集合のみを持つzdd ({λ}) を得る。
        """
        return self.one_terminal
3. ZDD.single_literal_set(self, v)
  • ZDDクラス で宣言します。
  • 下図が求めたいzddです。

 \{v\} の zdd



    def single_literal_set(self, v):
        """
            組合せ集合vのみを持つzdd ({v}) を得る。
        """
        assert v in self.idx
        return self._get_node(v, self.zero_terminal, self.one_terminal)
4. Diagram.offset(self, v)
  • Diagramクラス で宣言します。
  •  \mathcal{P} := \text{self} とおきます。
  •  \mathcal{P}\text{.top} = v の場合:

\begin{align} \mathcal{P}\text{.offset}(v) = \mathcal{P}\text{.zero} \end{align}

  •  \mathcal{P}\text{.top} の添え字が  v の添え字より大きい場合:

\begin{align} \mathcal{P}\text{.offset}(v) = \mathcal{P} \end{align}

  •  \mathcal{P}\text{.top} の添え字が  v の添え字より小さい場合:

\begin{align} & \mathcal{P}\text{.offset}(v) \\ & := \{C \in \mathcal{P} \mathrel{|} v \notin C\} \\ & = \{C \in \mathcal{P}\text{.zero} \mathrel{|} v \notin C\} \cup \{C \in (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) \mathrel{|} v \notin C\} \\ & = \{C \in \mathcal{P}\text{.zero} \mathrel{|} v \notin C\} \cup (\{C \in \mathcal{P}\text{.one} \mathrel{|} v \notin C\} \triangleleft \mathcal{P}\text{.top}) \\ & = \mathcal{P}\text{.zero.offset}(v) \cup (\mathcal{P}\text{.one.offset}(v) \triangleleft \mathcal{P}\text{.top}) \end{align}

    @cache
    def offset(self, v):
        """
            selfにおいて、vを含まない組合せ集合全体をなす部分zddを得る。
            例:
                self := {abc, ab, c}
                のとき
                self.offset(c) = {ab}
        """
        assert v in self.zdd.idx
        if self.top == v:
            return self.zero
        if self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self
        return self.zdd._get_node(self.top, self.zero.offset(v), self.one.offset(v))
5. Diagram.onset(self, v)
  • Diagramクラス で宣言します。
  •  \mathcal{P} := \text{self} とおきます。
  •  \mathcal{P}\text{.top} = v の場合:

\begin{align} \mathcal{P}\text{.onset}(v) = \mathcal{P}\text{.one} \end{align}

  •  \mathcal{P}\text{.top} の添え字が  v の添え字より大きい場合:

\begin{align} \mathcal{P}\text{.onset}(v) = \{\} \end{align}

  •  \mathcal{P}\text{.top} の添え字が  v の添え字より小さい場合:

\begin{align} & \mathcal{P}\text{.onset}(v) \\ & :=  \{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}, v \in C\} \\ & = \{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}\text{.zero}, v \in C\} \cup \{ C \setminus \{v\} \mathrel{|} C \in (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}), v \in C\} \\ & =  \{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}\text{.zero}, v \in C\} \cup (\{ C \setminus \{v\} \mathrel{|} C \in \mathcal{P}\text{.one}, v \in C\}  \triangleleft \mathcal{P}\text{.top}) \\ & = \mathcal{P}\text{.zero.onset}(v) \cup (\mathcal{P}\text{.one.onset}(v) \triangleleft \mathcal{P}\text{.top}) \end{align}

    @cache
    def onset(self, v):
        """
            selfにおいて、vを含む各組合せ集合からvを除いて取り出したzddを得る。
            例:
                self := {abc, ab, c}
                のとき
                self.onset(c) = {ab, λ}
        """
        assert v in self.zdd.idx
        if self.zdd.idx[self.top] == self.zdd.idx[v]:
            return self.one
        if self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self.zdd.zero_terminal
        return self.zdd._get_node(self.top, self.zero.onset(v), self.one.onset(v))
6. Diagram.change(self, v)
  • Diagramクラス で宣言します。
  •  \mathcal{P} := \text{self} とおきます。
  •  \mathcal{P}\text{.top} = v の場合:

\begin{align} \mathcal{P}\text{.change}(v) = \mathcal{P}\text{.one} \cup (\mathcal{P}\text{.zero} \triangleleft v) \end{align}

  •  \mathcal{P}\text{.top} の添え字が  v の添え字より大きい場合:

\begin{align} \mathcal{P}\text{.change}(v) = \{\} \cup (\mathcal{P} \triangleleft v) \end{align}

  •  \mathcal{P}\text{.top} の添え字が  v の添え字より小さい場合:

\begin{align} & \mathcal{P}\text{.change}(v) \\ & := \{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}, v \in C\} \cup \{C \cup \{v\} \mathrel{|} C \in \mathcal{P}, v \notin C\} \\ & = (\{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}\text{.zero}, v \in C\} \cup \{C \setminus \{v\} \mathrel{|} C \in (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}), v \in C\}) \\ & \cup (\{C \cup \{v\} \mathrel{|} C \in \mathcal{P}\text{.zero}, v \notin C\} \cup \{C \cup \{v\} \mathrel{|} C \in (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}), v \notin C\}) \\ & = (\{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}\text{.zero}, v \in C\} \cup \{C \cup \{v\} \mathrel{|} C \in \mathcal{P}\text{.zero}, v \notin C\}) \\ & \cup ( (\{C \setminus \{v\} \mathrel{|} C \in \mathcal{P}\text{.one}, v \in C\} \cup \{C \cup \{v\} \mathrel{|} C \in \mathcal{P}\text{.one}, v \notin C\}) \triangleleft \mathcal{P}\text{.top}) \\ & = \mathcal{P}\text{.zero.change}(v) \cup (\mathcal{P}\text{.one.change}(v) \triangleleft \mathcal{P}\text{.top}) \end{align}

    @cache
    def change(self, v):
        """
            selfの各組合せ集合について、vが含まれて いる/いない を反転させたzddを得る。
            例:
                self := {abc, ab, c}
                のとき
                self.change(b) = {ac, a, bc}
        """
        assert v in self.zdd.idx
        if self.top == v:
            return self.zdd._get_node(v, self.one, self.zero)
        if self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self.zdd._get_node(v, self.zdd.zero_terminal, self)
        return self.zdd._get_node(self.top, self.zero.change(v), self.one.change(v))
7. P | Q
  • Diagramクラス で宣言します。
  • 一方が空ならもう一方を、同じならどちらか一方を返せばよいです。
  •  \mathcal{P}\text{.top}, \mathcal{Q}\text{.top} の添え字の大小を比較して、
  •  \mathcal{P} の方が小さい場合: 

\begin{align} & \mathcal{P} \cup \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) \cup \mathcal{Q} \\ & =  ( \mathcal{P}\text{.zero} \cup \mathcal{Q} ) \cup ( \mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) \end{align} 

  •  \mathcal{Q} の方が小さい場合: 

\begin{align} & \mathcal{P} \cup \mathcal{Q} \\ & = \mathcal{P} \cup (\mathcal{Q}\text{.zero} \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{Q}\text{.top}) ) \\ & = (\mathcal{P} \cup \mathcal{Q}\text{.zero}) \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{Q}\text{.top}) \end{align} 

  • 等しい場合:

\begin{align} & \mathcal{P} \cup \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) \cup (\mathcal{Q}\text{.zero} \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) \\ & = (\mathcal{P}\text{.zero} \cup \mathcal{Q}\text{.zero}) \cup ( (\mathcal{P}\text{.one} \cup \mathcal{Q}\text{.one}) \triangleleft \mathcal{P}\text{.top}) \end{align}

    @cache
    def __or__(self, other):
        """
            self ∪ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self | other = {abc, ab, bc, c, λ}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self is self.zdd.zero_terminal:
            return other
        if other is self.zdd.zero_terminal or self is other:
            return self
        if self.zdd.idx[self.top] < self.zdd.idx[other.top]:
            return self.zdd._get_node(self.top, self.zero | other, self.one)
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return self.zdd._get_node(other.top, self | other.zero, other.one)
        return self.zdd._get_node(self.top, self.zero | other.zero, self.one | other.one)
8. P & Q
  • Diagramクラス で宣言します。
  • 一方が空なら0-葉を、同じならどちらか一方を返せばよいです。
  •  \mathcal{P}\text{.top}, \mathcal{Q}\text{.top} の添え字の大小を比較して、
  •  \mathcal{P} の方が小さい場合: 

\begin{align} & \mathcal{P} \cap \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) \cap \mathcal{Q} \\ & = \mathcal{P}\text{.zero} \cap \mathcal{Q} \end{align}

  •  \mathcal{Q} の方が小さい場合:

\begin{align} & \mathcal{P} \cap \mathcal{Q} \\ & = \mathcal{P} \cap (\mathcal{Q}\text{.zero} \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{Q}\text{.top}) )  \\ & = \mathcal{P} \cap \mathcal{Q}\text{.zero} \end{align}

  • 等しい場合:

\begin{align} & \mathcal{P} \cap \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) \cap (\mathcal{Q}\text{.zero} \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{P}\text{.top}) )  \\ & = (\mathcal{P}\text{.zero} \cap \mathcal{Q}\text{.zero}) \cup ( (\mathcal{P}\text{.one} \cap \mathcal{Q}\text{.one}) \triangleleft \mathcal{P}\text{.top})\end{align}

    @cache
    def __and__(self, other):
        """
            self ∩ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self & other = {ab}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self is self.zdd.zero_terminal or other is self.zdd.zero_terminal:
            return self.zdd.zero_terminal
        if self is other:
            return self
        if self.zdd.idx[self.top] < self.zdd.idx[other.top]:
            return self.zero & other
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return self & other.zero
        return self.zdd._get_node(self.top, self.zero & other.zero, self.one & other.one)
9. P + Q
  • Diagramクラス で宣言します。
  •  \mathcal{P} | \mathcal{Q} と等価なので、そのまま呼び出します。
    def __add__(self, other):
        """
            (self | other と等価。)
            self ∪ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self + other = {abc, ab, bc, c, λ}
        """
        assert isinstance(other, Obsolete_Diagram)
        return self.__or__(other)
10. P - Q
  • Diagramクラス で宣言します。
  •  \mathcal{P} = \{\} または  \mathcal{P} = \mathcal{Q} なら  \{\} を返せばよいです。
  •  \mathcal{Q} = \{\} なら  \mathcal{P} をそのまま返せばよいです。
  •  \mathcal{P}\text{.top}, \mathcal{Q}\text{.top} の添え字の大小を比較して、
  •  \mathcal{P} の方が小さい場合: 

\begin{align} & \mathcal{P} - \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) - \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} - \mathcal{Q}) \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) \end{align}

  •  \mathcal{Q} の方が小さい場合:

\begin{align} & \mathcal{P} - \mathcal{Q} \\ & = \mathcal{P} - (\mathcal{Q}\text{.zero} \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{Q}\text{.top}) ) \\ & = \mathcal{P} - \mathcal{Q}\text{.zero} \end{align}

  • 等しい場合:

\begin{align} & \mathcal{P} - \mathcal{Q} \\ & = (\mathcal{P}\text{.zero} \cup (\mathcal{P}\text{.one} \triangleleft \mathcal{P}\text{.top}) ) - (\mathcal{Q}\text{.zero} \cup (\mathcal{Q}\text{.one} \triangleleft \mathcal{P}\text{.top}) \\ & = (\mathcal{P}\text{.zero} - \mathcal{Q}\text{.zero}) \cup ( (\mathcal{P}\text{.one} - \mathcal{Q}\text{.one}) \triangleleft \mathcal{P}\text{.top}) \end{align}

    @cache
    def __sub__(self, other):
        """
            self \ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self - other = {abc, c}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self is self.zdd.zero_terminal or self is other:
            return self.zdd.zero_terminal
        if other is self.zdd.zero_terminal:
            return self
        if self.zdd.idx[self.top] < self.zdd.idx[other.top]:
            return self.zdd._get_node(self.top, self.zero - other, self.one)
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return self - other.zero
        return self.zdd._get_node(self.top, self.zero - other.zero, self.one - other.one)
11. P * Q
  • Diagramクラス で宣言します。
  • 明らかに可換なので、一般性を失わず  \mathcal{P}\text{.top} の添え字は  \mathcal{Q}\text{.top} の添え字以下であると仮定します。
  •  \mathcal{Q} = \{\} なら  \{\} を返せばよいです。
  •  \mathcal{Q} = \{ \lambda \} なら  \mathcal{P} を返せばよいです。
  •  v := \mathcal{P}\text{.top} とし、

\begin{align} \mathcal{P}_0 & :=  \mathcal{P}\text{.offset}(v) \\ \mathcal{P}_1 & :=  \mathcal{P}\text{.onset}(v) \\ \mathcal{Q}_0 & :=  \mathcal{Q}\text{.offset}(v) \\ \mathcal{Q}_1 &  := \mathcal{Q}\text{.onset}(v) \end{align}

とおきます。

\begin{align} & \mathcal{P} * \mathcal{Q} \\ & = (\mathcal{P}_0 \cup (\mathcal{P}_1 \triangleleft v) ) * (\mathcal{Q}_0 \cup (\mathcal{Q}_1 \triangleleft v) ) \\ & = (\mathcal{P}_0 * \mathcal{Q}_0) \cup ( ( (\mathcal{P}_0 * \mathcal{Q}_1) \cup (\mathcal{P}_1 * \mathcal{Q}_0) \cup (\mathcal{P}_1 * \mathcal{Q}_1) ) \triangleleft v) \end{align}

    @cache
    def __mul__(self, other):
        """
        「selfとotherから1つずつ取り、それらの和集合を取って得られる組合せ集合」全体を表すzddを得る。
        例:
            self := {ab, b, c}
            other := {ab, λ}
            のとき
            self * other = {ab, abc, b, c}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return other.__mul__(self)
        if other is self.zdd.zero_terminal:
            return self.zdd.zero_terminal
        if other is self.zdd.one_terminal:
            return self
        v = self.top
        p0, p1 = self.offset(v), self.onset(v)
        q0, q1 = other.offset(v), other.onset(v)
        return self.zdd._get_node(v, p0*q0,  p1*q1 + p1*q0 + p0*q1)
12. P // Q
  • Diagramクラス で宣言します。
  •  \mathcal{Q} = \{\lambda \} (乗法単位元)なら、 \mathcal{P} をそのまま返せばよいです。
  • (そうでなくて、) \mathcal{P} = \{\} または  \{\lambda \} なら、 \{\} を返せばよいです。
  •  \mathcal{P} = \mathcal{Q} なら  \{ \lambda \} を返せばよいです。
  •  v := \mathcal{Q}\text{.top} とし、

\begin{align} \mathcal{P}_0 & :=  \mathcal{P}\text{.offset}(v) \\ \mathcal{P}_1 & :=  \mathcal{P}\text{.onset}(v) \\ \mathcal{Q}_0 & :=  \mathcal{Q}\text{.offset}(v) \\ \mathcal{Q}_1 &  := \mathcal{Q}\text{.onset}(v) \end{align}

とおきます。

  •  \mathcal{Q}_0 = \{\} の場合:

\begin{align} \mathcal{P} // \mathcal{Q} =  \mathcal{P}_1 // \mathcal{Q}_1 \end{align}

  •  \mathcal{Q}_0 \neq \{\} の場合:

\begin{align} \mathcal{P} // \mathcal{Q} = (\mathcal{P}_0 // \mathcal{Q}_0) \cap (\mathcal{P}_1 // \mathcal{Q}_1) \end{align}

    @cache
    def __floordiv__(self, other):
        """
            other ≠ {} (空集合)とする。
            otherの各要素Cに対し「selfのCを含む各組合せ集合からCを除いて取り出した集合全体」を取ったときの、
            それらの共通部分を表すzddを得る。
            例:
                self := {abc, abd, ac, cd}
                other := {ab, c}
                のとき
                self//other
                    = {c, d} ∩ {ab, a, d}
                    = {d}
        """
        assert isinstance(other, Obsolete_Diagram) and other is not self.zdd.zero_terminal
        if other is self.zdd.one_terminal:
            return self
        if self is self.zdd.zero_terminal or self is self.zdd.one_terminal:
            return self.zdd.zero_terminal
        if self is other:
            return self.zdd.one_terminal
        v = other.top
        p0, p1 = self.offset(v), self.onset(v)
        q0, q1 = other.offset(v), other.onset(v)
        r = p1//q1
        if r is not self.zdd.zero_terminal and q0 is not self.zdd.zero_terminal:
            r = r & p0//q0
        return r
13. P % Q
  • Diagramクラス で宣言します。
  •  -, *, // は既に実装したので、定義通りに計算するだけです。
    def __mod__(self, other):
        """
            other ≠ {} (空集合)とする。
            self % other := self - other*(self/other)
            例:
                self := {abc, abd, ac, cd}
                other := {ab, c}
                のとき
                self % other
                    = {abc, abd, ac, cd} - {ab, c}*{d}
                    = {abc, abd, ac, cd} - {abd, cd}
                    = {abc, ac}
        """
        assert isinstance(other, Obsolete_Diagram) and other is not self.zdd.zero_terminal
        return self - other*(self//other)

 

これにて、演算の実装は完了です。これらをうまく組み合わせて、様々な部分集合族に対応するZDDをボトムアップに構築することができます。

 

組合せ集合の総数を求める

ZDDに含まれる組合せ集合は、根から1-葉までのパスに1対1対応します。

ということは、組合せ集合の総数を求めたければパスの総数を求めればよいということになります。DAGなので、メモ化再帰(またはDP)でよいですね。

    @cache
    def count(self):
        """
            selfに含まれる組合せ集合の個数を得る。
            例:
                self := {abc, ab, c}
                のとき
                self.count() = 3
        """
        if self is self.zdd.zero_terminal:
            return 0
        if self is self.zdd.one_terminal:
            return 1
        return self.zero.count() + self.one.count()

 

ここまで実装まとめ

import sys
sys.setrecursionlimit(1000000000) # ごめんなさい
from pyvis.network import Network
from collections import deque
from random import randint

class Obsolete_Diagram:
    def __init__(self, zdd, top=None, zero=None, one=None) -> None:
        self.zdd = zdd
        self.top = top
        self.zero = zero
        self.one = one
        if top is None:
            self.hash = randint(0, 10**9)
        else:
            self.hash = (top, zero.hash, one.hash).__hash__()
    
    def __hash__(self) -> int:
        return self.hash
    
    def show(self, filename):
        """
            selfによるzddを図示
            filemame: html拡張子のファイル名
        """
        nodes = {self: 0}
        net = Network(directed=True)
        dq = deque()
        if self is self.zdd.zero_terminal:
            net.add_node(nodes[self], shape="box", label="⏊")
        elif self is self.zdd.one_terminal:
            net.add_node(nodes[self], shape="box", label="⏉")
        else:
            net.add_node(nodes[self], shape="circle", label=self.top.__repr__())
            dq.appendleft(self)
        while dq:
            p = dq.pop()
            if p.zero not in nodes:
                nodes[p.zero] = len(nodes)
                if p.zero is self.zdd.zero_terminal:
                    net.add_node(nodes[p.zero], shape="box", label="⏊")
                elif p.zero is self.zdd.one_terminal:
                    net.add_node(nodes[p.zero], shape="box", label="⏉")
                else:
                    net.add_node(nodes[p.zero], shape="circle", label=p.zero.top.__repr__())
                    dq.appendleft(p.zero)
            net.add_edge(nodes[p], nodes[p.zero], color="red")
            if p.one not in nodes:
                nodes[p.one] = len(nodes)
                if p.one is self.zdd.zero_terminal:
                    net.add_node(nodes[p.one], shape="box", label="⏊")
                elif p.one is self.zdd.one_terminal:
                    net.add_node(nodes[p.one], shape="box", label="⏉")
                else:
                    net.add_node(nodes[p.one], shape="circle", label=p.one.top.__repr__())
                    dq.appendleft(p.one)
            net.add_edge(nodes[p], nodes[p.one], color="blue")
        net.show(filename, notebook=False)

    @cache
    def offset(self, v):
        """
            selfにおいて、vを含まない組合せ集合全体をなす部分zddを得る。
            例:
                self := {abc, ab, c}
                のとき
                self.offset(c) = {ab}
        """
        assert v in self.zdd.idx
        if self.top == v:
            return self.zero
        if self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self
        return self.zdd._get_node(self.top, self.zero.offset(v), self.one.offset(v))
    
    @cache
    def onset(self, v):
        """
            selfにおいて、vを含む各組合せ集合からvを除いて取り出したzddを得る。
            例:
                self := {abc, ab, c}
                のとき
                self.onset(c) = {ab, λ}
        """
        assert v in self.zdd.idx
        if self.zdd.idx[self.top] == self.zdd.idx[v]:
            return self.one
        if self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self.zdd.zero_terminal
        return self.zdd._get_node(self.top, self.zero.onset(v), self.one.onset(v))
    
    @cache
    def change(self, v):
        """
            selfの各組合せ集合について、vが含まれて いる/いない を反転させたzddを得る。
            例:
                self := {abc, ab, c}
                のとき
                self.change(b) = {ac, a, bc}
        """
        assert v in self.zdd.idx
        if self.top == v:
            return self.zdd._get_node(v, self.one, self.zero)
        if self.zdd.idx[self.top] > self.zdd.idx[v]:
            return self.zdd._get_node(v, self.zdd.zero_terminal, self)
        return self.zdd._get_node(self.top, self.zero.change(v), self.one.change(v))
    
    @cache
    def __or__(self, other):
        """
            self ∪ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self | other = {abc, ab, bc, c, λ}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self is self.zdd.zero_terminal:
            return other
        if other is self.zdd.zero_terminal or self is other:
            return self
        if self.zdd.idx[self.top] < self.zdd.idx[other.top]:
            return self.zdd._get_node(self.top, self.zero | other, self.one)
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return self.zdd._get_node(other.top, self | other.zero, other.one)
        return self.zdd._get_node(self.top, self.zero | other.zero, self.one | other.one)
    
    @cache
    def __and__(self, other):
        """
            self ∩ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self & other = {ab}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self is self.zdd.zero_terminal or other is self.zdd.zero_terminal:
            return self.zdd.zero_terminal
        if self is other:
            return self
        if self.zdd.idx[self.top] < self.zdd.idx[other.top]:
            return self.zero & other
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return self & other.zero
        return self.zdd._get_node(self.top, self.zero & other.zero, self.one & other.one)
    
    def __add__(self, other):
        """
            (self | other と等価。)
            self ∪ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self + other = {abc, ab, bc, c, λ}
        """
        assert isinstance(other, Obsolete_Diagram)
        return self.__or__(other)

    @cache
    def __sub__(self, other):
        """
            self \ other を表すzddを得る。
            例:
                self := {abc, ab, c}
                other := {ab, bc, λ}
                のとき
                self - other = {abc, c}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self is self.zdd.zero_terminal or self is other:
            return self.zdd.zero_terminal
        if other is self.zdd.zero_terminal:
            return self
        if self.zdd.idx[self.top] < self.zdd.idx[other.top]:
            return self.zdd._get_node(self.top, self.zero - other, self.one)
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return self - other.zero
        return self.zdd._get_node(self.top, self.zero - other.zero, self.one - other.one)
    
    @cache
    def __mul__(self, other):
        """
        「selfとotherから1つずつ取り、それらの和集合を取って得られる組合せ集合」全体を表すzddを得る。
        例:
            self := {ab, b, c}
            other := {ab, λ}
            のとき
            self * other = {ab, abc, b, c}
        """
        assert isinstance(other, Obsolete_Diagram)
        if self.zdd.idx[self.top] > self.zdd.idx[other.top]:
            return other.__mul__(self)
        if other is self.zdd.zero_terminal:
            return self.zdd.zero_terminal
        if other is self.zdd.one_terminal:
            return self
        v = self.top
        p0, p1 = self.offset(v), self.onset(v)
        q0, q1 = other.offset(v), other.onset(v)
        return self.zdd._get_node(v, p0*q0,  p1*q1 + p1*q0 + p0*q1)

    @cache
    def __floordiv__(self, other):
        """
            other ≠ {} (空集合)とする。
            otherの各要素Cに対し「selfのCを含む各組合せ集合からCを除いて取り出した集合全体」を取ったときの、
            それらの共通部分を表すzddを得る。
            例:
                self := {abc, abd, ac, cd}
                other := {ab, c}
                のとき
                self//other
                    = {c, d} ∩ {ab, a, d}
                    = {d}
        """
        assert isinstance(other, Obsolete_Diagram) and other is not self.zdd.zero_terminal
        if other is self.zdd.one_terminal:
            return self
        if self is self.zdd.zero_terminal or self is self.zdd.one_terminal:
            return self.zdd.zero_terminal
        if self is other:
            return self.zdd.one_terminal
        v = other.top
        p0, p1 = self.offset(v), self.onset(v)
        q0, q1 = other.offset(v), other.onset(v)
        r = p1//q1
        if r is not self.zdd.zero_terminal and q0 is not self.zdd.zero_terminal:
            r = r & p0//q0
        return r
    
    def __mod__(self, other):
        """
            other ≠ {} (空集合)とする。
            self % other := self - other*(self/other)
            例:
                self := {abc, abd, ac, cd}
                other := {ab, c}
                のとき
                self % other
                    = {abc, abd, ac, cd} - {ab, c}*{d}
                    = {abc, abd, ac, cd} - {abd, cd}
                    = {abc, ac}
        """
        assert isinstance(other, Obsolete_Diagram) and other is not self.zdd.zero_terminal
        return self - other*(self//other)

    @cache
    def count(self):
        """
            selfに含まれる組合せ集合の個数を得る。
            例:
                self := {abc, ab, c}
                のとき
                self.count() = 3
        """
        if self is self.zdd.zero_terminal:
            return 0
        if self is self.zdd.one_terminal:
            return 1
        return self.zero.count() + self.one.count()

class ZDD:

    def __init__(self, ground_set, Diagram=Obsolete_Diagram) -> None:
        self.idx = {v: i for i, v in enumerate(ground_set)}
        self.idx[None] = len(ground_set)
        self.Diagram = Diagram
        self.zero_terminal = self.Diagram(zdd=self)
        self.one_terminal = self.Diagram(zdd=self)

    @cache
    def _get_node(self, v, zero, one):
        if one is self.zero_terminal:
            return zero
        return self.Diagram(self, v, zero, one)
    
    def empty_set(self):
        """
            空のzdd ({}) を得る。
        """
        return self.zero_terminal
    
    def unit_set(self):
        """
            空の組合せ集合のみを持つzdd ({λ}) を得る。
        """
        return self.one_terminal

    def single_literal_set(self, v):
        """
            組合せ集合vのみを持つzdd ({v}) を得る。
        """
        assert v in self.idx
        return self._get_node(v, self.zero_terminal, self.one_terminal)

 

8-Queen problem

最後に、これらのZDD演算を利用して解ける有名な数え上げの問題を1つ紹介します。

エイト・クイーン - Wikipedia

 

8×8のチェス盤に、相互に攻撃が当たらないようクイーンを8個配置する方法は何通り?と言う問題です。

配置方法の1例(エイト・クイーン - Wikipedia)より

この問題は「各マスを台集合とし、クイーンを配置するマスを部分集合として取り出す」と捉えることで、解全体をZDDを用いて表現することができます。

 

8-Queen problem の解全体を表すZDDの構築方法

各行にはクイーンを1つ配置するので、上の行から順にどの列にクイーンを配置するかを決めていきながら、その行までで考えられる配置全体を表すZDDを構築していきます。マス  (i, j) grid_{ij} 表し、この全体を台集合とします。

 

 i-1 行目までのZDDを  \mathcal{S}_{i-1} として、新たにマス  (i, j) にクイーンを配置した場合の暫定ZDDを作るとします。

このマスにクイーンを配置できるのは上と斜め上の各マスにクイーンが置かれていないときに限るので、 \mathcal{S}_i からそれらのマスを含まない組合せ集合を取り出し、それらにそれぞれ  grid_{ij} を追加すればよいです。

つまり、それぞれのマスについて  \text{offset}() したあとに、 \text{change}(grid_{ij}) すればいいわけですね。

のマスを含まない組合せ集合を  S_{i-1} から抜き出す。

 \mathcal{S}_i は、各  j について上記手順で得られるZDDの和を取れば求まります。

実装例

8×8のチェス盤に限る必要はないので、一般のn×nのチェス盤で解けるようにしました。solve() を show = True で実行すればZDDを図示してくれます。

class QueenProblem:

    def solve(n: int, show: bool = False) -> int:
        """
            n×n チェス盤でのクイーン問題を解く。
            show: zddを表示するか -> 
        """
        assert type(n) is int and 1 <= n

        def grid(i, j):
            return i*n + j
        
        zdd = ZDD([grid(i, j) for i in range(n) for j in range(n)])
        s = zdd.empty_set()
        for j in range(n):
            s += zdd.single_literal_set(j)
        for i in range(1, n):
            t = zdd.empty_set()
            for j in range(n):
                u = s
                for k in range(i):
                    u = u.offset(grid(i-k-1, j))
                for k in range(min(i, j)):
                    u = u.offset(grid(i-k-1, j-k-1))
                for k in range(min(i, n-j-1)):
                    u = u.offset(grid(i-k-1, j+k+1))
                t += u.change(grid(i, j))
            s = t
        if show:
            s.show(f"Queen_problem{n}.html")
        return s.count()

n = 8 での ZDD

nごとの実行時間

せっかくなので、私の環境でどれくらいの速度で n-Queen problem が解けたかの記録を残しておきます。各nごとに10回の平均を取ります。

from time import perf_counter

def timer_return_time(func):
    """
        関数の実行時間を計測するwrapper
        実行時間(ms)を返り値に追加する
    """
    def wrapper(*arg, **darg):
        begin = perf_counter()
        ret = func(*arg, **darg)
        end = perf_counter()
        # print(f'{func.__name__}: {(end-begin)*1000}ms')
        return ret, (end-begin)*1000
    return wrapper

class QueenProblem:

    @timer_return_time
    def solve(n: int, show: bool = False) -> int:
        (中略)
        return s.count()

if __name__ == "__main__":
    repeat = 10
    for n in range(11, 12):
        t = 0
        for _ in range(repeat):
            ans, time = QueenProblem.solve(n, False)
            t += time
        print(f"{n}×{n}: {ans}個, {t/repeat}ms")

実行環境:

  • Python 3.11.4 (tags/v3.11.4:d2340ef, Jun  7 2023, 05:45:37) [MSC v.1934 64 bit (AMD64)] on win32
  • windows 11
  • Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz   2.30 GHz
  • 8.00 GB

実行結果:

n = 1 ~ 11 での実行結果。n = 12 はしばらく待っても終わりませんでした。

※ ちなみにバックトラックした方が普通に速いらしいです。

Nクイーン【前編】 | TECH PROjin

 

おわりに

お読みいただきありがとうございました。

卒研お勉強枠次回は、ZDDを使って某有名パズルゲームを解くお話になる予定です。気長にお待ちください。

 

次回 >>

 

参考文献

[1] Zero-suppressed BDDs and their applications (湊真一)

[2] 超高速グラフ列挙アルゴリズム|森北出版株式会社 (ERATO 湊離散構造処理系プロジェクト 著, 湊真一 編)

[3] ZDD (Zero-suppressed binary Decision Diagram) 超高速グラフ列挙アルゴリズム 著者: 湊 真一 (石井 健太)

[4] エイト・クイーン - Wikipedia

[5] Nクイーン【前編】 | TECH PROjin

*1:正確には削減していないものもZDDと呼び、完全に削減したものを既約なZDDと呼ぶ方が主流だそうです。また、既約なZDDは削除の順序に依らず一意に定まるそうですが、これの証明は調べてもよくわからなかったです…。

*2:Diagramクラスはどんどん継承して拡張する想定の設計となっています。

マトロイドによるクラスカル法の正当性の証明

 

はじめに

クラスカル法を知っていますか?私は知ってます。

クラスカル法 - Wikipedia

最小全域木問題の解法として知られるアルゴリズムで、重みの昇順に閉路ができないよう辺を追加していくだけで最小全域木になりますよ~という、貪欲法で最適解が得られる稀有(?)な例です。基本的に、貪欲法は「局所最適解に嵌ってしまう」という弱点があるのですが、クラスカル法は帰納法wikiのやつ)によって正当性が証明できます。

 

しかし、貪欲法が使えるかどうかを問題ごとに一々考えるのは大変です。そこで考えられた貪欲法が使える指標の一つが「マトロイド」と呼ばれる構造です。

 

もし「全域木がマトロイドだからクラスカル法は正しい」と雑に説明できてしまうと、とても楽ですよね。というわけで、これをやります。*1

 

この記事でやること

  • 離散最適化問題において貪欲法が使える「マトロイド」という構造について、クラスカル法の正当性の証明を目標に解説します。
  • 記事の前半ではマトロイドの定義と貪欲法との関係とその証明、後半ではクラスカル法にマトロイドを適用する方法の説明となります。これら説明を最短で行うため、具体例をほぼほぼ出さずに定義→命題→証明の流れで突き進みます。
  • あくまでもクラスカル法に関する話を目標に定めているため、マトロイドについてのほんの一部分にしか触れておりません(もちろん証明は厳密に書いているつもりです)。ですので、本記事を読んでマトロイドに興味を持たれた方は是非ご自身で色々調べてみることをお勧めします。

 

注意点

  • ちょっぴり線形代数の知識がいります。(線形独立とか次元とかが何かわかっていれば十分のはず…。)
  • グラフに単純性を仮定しません。(すなわち、多重辺・自己ループの存在を認めます。)

 

(前半)マトロイドと貪欲法

マトロイドの定義

【定義 1.1.1】*2
【定義 1.1.1】(独立性システム、マトロイド)
 E: 非空な有限集合, \mathcal{F}: E の部分集合族 とする。
 \mathcal{F} が次の2条件を満たすとき、 \mathcal{F} Eの独立性システムと言う: 
  1.  \emptyset \in \mathcal{F}
  2.  X \in \mathcal{F}, Y \subseteq X \Rightarrow Y \in \mathcal{F}

さらに、次の条件も満たすとき、 \mathcal{F} E のマトロイドと言う: 

  1.  X, Y \in \mathcal{F}, |X| \gt |Y| \Rightarrow \exists e \in X \setminus Y : Y \cup e \in \mathcal{F}

※以下、マトロイドに対しては  \mathcal{F} のかわりに主に  \mathcal{I} を使うこととします。

 

マトロイドに関する用語

【定義 1.2.1】
【定義 1.2.1】(独立集合、従属集合)
 E: 非空な有限集合, \mathcal{I}: E のマトロイド とする。
 X \subseteq E について
  •  X: (\mathcal{I}の)独立集合 : \underset{\text{def}}{\Leftrightarrow} X \in \mathcal{I}
  •  X: (\mathcal{I}の)従属集合 : \underset{\text{def}}{\Leftrightarrow} X \notin \mathcal{I}
【定義 1.2.2】
【定義 1.2.2】(基、サーキット)
 E: 非空な有限集合, \mathcal{I}: E のマトロイド とする。
 B \in \mathcal{I}, C \in 2^{E} \setminus \mathcal{I} について
  •  B: (\mathcal{I}の)基: \underset{\text{def}}{\Leftrightarrow} B: \mathcal{I}において極大
  •  C: (\mathcal{I}の)サーキット: \underset{\text{def}}{\Leftrightarrow} 2^{E}\backslash \mathcal{I}において極小
【定義 1.2.3】
【定義 1.2.3】(集合の基)
 E: 非空な有限集合, \mathcal{I}: E のマトロイド, X \subseteq E とする。
このとき、 \exists B_X \subseteq X: B_X \in \mathcal{I} かつ \forall e \in X \setminus B_X, B_X \cup \{e\} \notin \mathcal{I}
この B_X Xの基と呼ぶ。

※【定義 1.2.2】の「基」とは異なる概念であることに注意してください。

存在性は簡単に示せます。

 \mathcal{F} := \{Y \subseteq X \mathrel{|} Y \in \mathcal{I} \} としたときに  \emptyset \in \mathcal{F} より  \mathcal{F} \neq \emptyset であることに注意して、  B_X として  \mathcal{F} の中で最大のものを取ればよいです。

特に、 X \in \mathcal{I} のとき  B_X = X となります。

 

【定義 1.2.4】
【定義 1.2.4】(階数関数)
 E: 非空な有限集合, \mathcal{I}: E のマトロイド とする。
関数  r: 2^E \to \mathbb{N}_{\ge 0}; X \mapsto r(X)
\begin{align} r(X) := \text{max} \{ |Y| \mathrel{|} Y \subseteq X, Y \in \mathcal{I} \} \end{align}と定め、これを( \mathcal{I}の)階数関数と呼ぶ。
つまり、 r(X) Xの基のサイズのことです。集合の基のサイズは一意であるため、well-definedです。

 

【命題 1.2.5】*3
【命題 1.2.5】
 E: 非空な有限集合, \mathcal{I}: E のマトロイド, r: 階数関数 について、以下が成立:
  1.  X \subseteq E \Rightarrow r(X) \le |X|
  2.  X \in \mathcal{I} \Leftrightarrow r(X) = |X|
  3.  X \subseteq Y \subseteq E \Rightarrow r(X) \le r(Y)

直感的に明らかだと思うので、証明略*4

 

マトロイドと貪欲法

はじめにざっくりと「貪欲法が使える」と述べましたが、これを正しく定式化しておきます。

【定義 1.3.1】
【定義 1.3.1】(最大独立集合問題)
 E: 非空な有限集合, \mathcal{F}: E の独立性システム とする。
与えられた重み関数  w: E \to \mathbb{R}_{\ge 0} に対し、
 \sum_{e \in X} w(e) が最大となる  X \in \mathcal{F} を見つける最適化問題を最大独立集合問題と呼ぶ。

※例えば「01ナップサック問題」や「最大マッチング問題」などは、最大独立集合問題として表現できます。( \mathcal{F} をうまく定義してあげればよいです。)

【定義 1.3.2】
【定義 1.3.2】(最小基問題)
 E: 非空な有限集合, \mathcal{F}: E の独立性システム とする。
与えられた重み関数  w: E \to \mathbb{R}_{\ge 0} に対し、
 \sum_{e \in X} w(e) が最小となる極大な  X \in \mathcal{F} を見つける最適化問題を最小基問題と呼ぶ。

※「最小全域木問題」は最小基問題として表現できます。

 

これら問題に対し、貪欲法を疑似コードで表現します。

【定義 1.3.3】
【定義 1.3.3】(最大独立集合問題に対する貪欲法)
 E := \{ e_{1}, \dots, e_{n} \}, \mathcal{F}: E の独立性システム, w: E \to \mathbb{R}_{\ge 0}: 重み関数
による最大独立集合問題について、(必要があれば適切に並び替えて) w(e_{1}) \ge \dots \ge w(e_{n}) としたとき、貪欲法を以下のアルゴリズムで定義する:
  1.  X \gets \emptyset
  2.  \text{for } i = 1, \dots, n:
    1.  \text{if } X \cup e_{i} \in \mathcal{F}: X \gets X \cup e_{i}
  3.  \text{return } X
【定義 1.3.4】
【定義 1.3.4】(最小基問題に対する貪欲法)
 E := \{ e_{1}, \dots, e_{n} \}, \mathcal{F}: E の独立性システム, w: E \to \mathbb{R}_{\ge 0}: 重み関数
による最小基問題について、(必要があれば適切に並び替えて) w(e_{1}) \le \dots \le w(e_{n}) としたとき、貪欲法を以下のアルゴリズムで定義する:
  1.  X \gets \emptyset
  2.  \text{for } i = 1, \dots, n:
    1.  \text{if } X \cup e_{i} \in \mathcal{F}: X \gets X \cup e_{i}
  3.  \text{return } X

※ここで出力される  X が極大であることは【定義 1.1.1-(Ⅱ)】から示せます。

 

勿論、これら貪欲法は一般にはうまくいきません。
しかし、ここで以下の定理が成り立ちます。

 

【定理 1.3.5】
【定理 1.3.5】(マトロイドの最大独立集合問題に対する貪欲法)
 E := \{ e_{1}, \dots, e_{n} \}, \mathcal{F}: E の独立性システム とする。

 \mathcal{F}: マトロイド
 \Leftrightarrow 任意の重み関数  w: E \to \mathbb{R}_{\ge 0} に対し、【定義 1.3.3】の貪欲法が最大独立集合問題の最適解を出力する。

【定理 1.3.6】
【定理 1.3.6】(マトロイドの最小基問題に対する貪欲法)
 E := \{ e_{1}, \dots, e_{n} \}, \mathcal{F}: E の独立性システム とする。

 \mathcal{F}: マトロイド
 \Rightarrow 任意の重み関数  w: E \to \mathbb{R}_{\ge 0} に対し、【定義 1.3.4】の貪欲法が最小基問題の最適解を出力する。

※【定理 1.3.6】については、一般に逆が成り立つのか私自身わかっておりません。勉強不足で申し訳ありません…

 

早速証明していきましょう!

 

【定理 1.3.5】の証明

  •  \Rightarrowの証明:( \mathcal{F}: Eのマトロイド を仮定)
補題 1.4.1】

補題 1.4.1】

貪欲法の出力を B、最適解の1つを B^{*}とする。

また、 i = 1, \dots , n に対し

 B_i := B \cap \{e_1, \dots , e_i \}, B^{*}_i := B^{*} \cap  \{e_1, \dots , e_i \} とする。このとき、 r: 階数関数 として

  1.  |B_i| = r(\{e_1, \dots , e_i \})

  2.  |B^{*}_i| \le r(\{e_1, \dots , e_i \})

  3.  |B_i| \ge |B^{*}_i|

証明:

 B, B^{*}: \mathcal{F}の基 であることに注意する。

i:

  •  B_i: \{e_1, \dots , e_i \}の基 を示せば、【命題1.2.5-(ii)】より示される。
    •  B_i \in \mathcal{F} である。
    •  B_iアルゴリズム2行目のforループ終了時における Xに等しいため、
    •  \forall e_j \in \{e_1, \dots , e_i \} \setminus B_i, B_i \cup \{e_j\} \notin \mathcal{F}
    •  \therefore B_i: \{e_1, \dots , e_i \}の基 \because【定義 1.2.3】)   \square

ii: 

  •  B^*_i \subseteq B^* より  B^*_i \in \mathcal{F} \because【定義 1.1.1-(Ⅱ)】)であるため

\begin{align} |B^*_i| &= r(B^*_i) &\because \text{【命題 1.2.5-(ii)】} \\ &\le r(\{e_1, \dots , e_i \}) &\because \text{【命題 1.2.5-(i)】} \square \end{align}

iii:

  • i, ii と【命題 1.2.5-(iii)】より  \square

 

 B_0 = B^*_0 := \emptyset, w(e_{n+1}) := 0 とすると

\begin{align} \sum_{e \in B} w(e) &= \sum_{i=1}^{n} (|B_i|-|B_{i-1}|)w(e_i) \\ &= \sum_{i=1}^{n} |B_i|w(e_i) - \sum_{i=1}^{n} |B_{i-1}|w(e_i) \\ &= \sum_{i=1}^{n} |B_i|w(e_i) - \sum_{i=0}^{n-1} |B_i|w(e_{i+1}) \\ &= \sum_{i=1}^{n} |B_i|w(e_i) - \sum_{i=1}^{n} |B_i|w(e_{i+1}) \\ &= \sum_{i=1}^{n} |B_i|(w(e_i)-w(e_{i+1})) \\ &\ge \sum_{i=1}^{n} |B^*_i|(w(e_i)-w(e_{i+1})) \tag{※} \\ &= \sum_{i=1}^{n} |B^*_i|w(e_i) - \sum_{i=1}^{n} |B^*_i|w(e_{i+1}) \\ &= \sum_{i=1}^{n} |B^*_i|w(e_i) - \sum_{i=0}^{n-1} |B^*_i|w(e_{i+1}) \\ &= \sum_{i=1}^{n} |B^*_i|w(e_i) - \sum_{i=1}^{n} |B^*_{i-1}|w(e_i) \\ &= \sum_{i=1}^{n} (|B^*_i|-|B^*_{i-1}|)w(e_i) \\ &= \sum_{e \in B^*} w(e) \end{align}

(※)では【補題 1.4.1】と  w(e_i)-w(e_{i+1}) \ge 0 を用いた。

以上より、貪欲法による解 Bは最適解の1つである。  \square

 

  •  \Leftarrowの証明:(対偶を示す。 \mathcal{F}: Eのマトロイドでない を仮定)
    •  \mathcal{F}は独立性システムであるため、仮定より【定義 1.1.1-(Ⅲ)】を満たさない、すなわち  \exists X, Y \in \mathcal{F}, |X| \gt |Y| : \forall e \in X \setminus Y, Y \cup \{e\} \notin \mathcal{F}
    • ここで、  k := |Y| とおき、次のように重み関数 w: E \to \mathbb{R}_{\ge 0} を定義:

\begin{align} w(e) := \begin{cases} k+2 & (e \in Y) \\ k+1 & (e \in X \setminus Y) \\ 0 & (\text{otherwise}) \end{cases} \end{align}

    •  B を貪欲法による解とすると、 \sum_{e \in B} w(e) = k(k+2)
      •  \because貪欲法のアルゴリズムの動作中、まず  B  e \in Y を満たす  e が全て追加される。
      • このとき、 X, Y に関する条件より  e \in X \setminus Y を満たす  e は追加されない。
      • その後、いくつかの  w(e) = 0 である  e が追加される。
      • よって、 \sum_{e \in B} w(e) = |Y|(k+2) = k(k+2)  \square
    • 一方、 X \in \mathcal{F} であって

\begin{align} \sum_{e \in X} w(e) &= |X|(k+2) \\ &\ge (|Y|+1)(k+1) \\ &= (k+1)(k+1) \\ &\gt k(k+2) \end{align}

    • よって、この重み関数  w に対して貪欲法は最適解を出力しない  \square

 

【定理 1.3.6】の証明

  •  \mathcal{F}: Eのマトロイド を仮定。
  •  W := \sum_{e \in E}w(e) とおき、関数  \tilde{w}: E \to \mathbb{R}_{\ge 0}; e \mapsto \tilde{w}(e) \tilde{w}(e) := W-w(e) で定める。
補題 1.4.2】

補題 1.4.2】

 w による最小基問題の最適値を  x \tilde{w} による最大独立集合問題の最適値を  y とし、 \mathcal{F} の階数関数を  r とする。このとき、

 x = Wr(E)-y

証明:

\begin{align} x &= \underset{B: E\text{の基}}{\text{min}} \{ \sum_{e \in B}w(e) \} \\ &= \underset{B: E\text{の基}}{\text{min}} \{ \sum_{e \in B}(W-\tilde{w}(e)) \} \\ &= \underset{B: E\text{の基}}{\text{min}} \{ W|B| - \sum_{e \in B}(\tilde{w}(e)) \} \\ &= Wr(E) - \underset{B: E\text{の基}}{\text{max}} \{ \sum_{e \in B}(\tilde{w}(e)) \} \\ &= Wr(E) - \underset{X \in \mathcal{F}}{\text{max}} \{ \sum_{e \in B}(\tilde{w}(e)) \} \\ &= Wr(E)-y \quad \square \end{align}

 

  •  \tilde{w}(e_1) \ge \dots \ge \tilde{w}(e_n) に注意すると、 w による最小基問題に対しての貪欲法による解  B は、  \tilde{w} による最大独立集合問題の最適解の1つとなっている。
  • すなわち、【補題 1.4.2】において  y = \sum_{e \in B}\tilde{w}(e) と表せる。
  • ここで

\begin{align} \sum_{e \in B}\tilde{w}(e) &= \sum_{e \in B}(W-w(e)) \\ &= Wr(E) - \sum_{e \in B}w(e) \end{align}

  • なので、【補題 1.4.2】の  y にこれを代入すれば

\begin{align} x &= Wr(E)-(Wr(E) - \sum_{e \in B}w(e)) \\ &= \sum_{e \in B}w(e) \end{align}

  • したがって、貪欲法による解  B w による最小基問題の最適解でもある。  \square

 

長い証明を経て、ついに前半戦は終了です!

【定義 1.1.1】のたったの3条件から、貪欲法が使えるという凄まじい結論が得られるの、鮮やか過ぎると思いません?私は思いました。

 

さて、ここからはグラフとマトロイドの関係に迫る後半戦です。

 

(後半)マトロイドのクラスカル法への適用

ベクトルマトロイド

非常に重要なマトロイドです。グラフとマトロイドを繋ぐ要となります。

そのまま競プロでも使います。

【命題 2.1.1】
【命題 2.1.1】(ベクトルマトロイド)
 \boldsymbol{a}_{1}, \dots , \boldsymbol{a}_{n}: (適当な体上で定義された)m次元ベクトル,
 E := \{ \boldsymbol{a}_{1}, \dots , \boldsymbol{a}_{m} \} とする。
 E の部分集合族  \mathcal{I} を次で定義する:
 \mathcal{I} := \{ X \subset E \mathrel{|} X: 線形独立 \}
このとき  \mathcal{I} はマトロイドであり、ベクトルマトロイドと呼ぶ。 

証明:

以下の補題を用いる: (線形代数の知識なので証明は略*5

補題

 A, B: m次元ベクトルの集合 について

  1.  A: 線形独立 \Leftrightarrow \text{dim} \langle A \rangle = |A|
  2.  A \subseteq B \Rightarrow \text{dim} \langle A \rangle \le \text{dim} \langle B \rangle
  3.  A \subseteq B のとき、  \langle A \rangle = \langle B \rangle \Leftrightarrow \text{dim} \langle A \rangle = \text{dim} \langle B \rangle
  • マトロイドの3条件を確認する。
    1. ok
    2.  Y \subset X \in \mathcal{I} とし、 Y: 線形独立 を示す。
      •  X = \{\boldsymbol{a}_1, \dots, \boldsymbol{a}_k\}, Y = \{\boldsymbol{b}_{i_1}, \dots, \boldsymbol{b}_{i_l}\} \quad (\{i_1, \dots i_l\} \subset \{1, \dots, k\}) とし、 \Delta := \{1, \dots, k\} \setminus \{i_1, \dots, i_l\} とおく。

\begin{align} \sum_{j=1}^{l} \lambda_{i_j} \boldsymbol{a}_{i_j} = \boldsymbol{0} &\Rightarrow \sum_{j=1}^{l} \lambda_{i_j} \boldsymbol{a}_{i_j} + \sum_{i \in \Delta } 0 \cdot \boldsymbol{a}_i = \boldsymbol{0} \\ &\Rightarrow \sum_{i=1}^{k} \lambda_i \boldsymbol{a}_i = \boldsymbol{0} \quad (\lambda_i := 0 \text{ for } i \in \Delta ) \\ &\Rightarrow \lambda_i = 0 \quad \text{for } i \in \{1, \dots, k\} (\because X: \text{線形独立}) \\ &\Rightarrow \lambda_{i_j} = 0 \quad \text{for } j \in \{1, \dots, l\} \end{align}

      • ok
    1.  X, Y \in \mathcal{I}, |X| \gt |Y| とし、 Y: 線形独立 を示す。
      •  X, Y: 線形独立 より、 \text{dim} \langle X \rangle = |X|, \text{dim} \langle Y \rangle = |Y| \because補題-(i)】)
      • 以下、背理法で示す。 \forall x \in X \setminus Y, Y \cup \{x\}: 線形従属 と仮定。
      • このとき、 \forall x \in X, \langle Y \cup \{x\} \rangle = \langle Y \rangle …(※)
        •  \because)
        •  x \in X \cap Y の場合:
          •  Y \cup \{x\} = Y より  \square
        •  x \notin X \cap Y の場合:

\begin{align} |Y| &= \text{dim} \langle Y \rangle \\ &\le \text{dim} \langle Y \cup \{ x \} \rangle &\because 【補 題\text{-(ii)}】 \\ &\lt | Y \cup \{ x \} | &\because 【補 題\text{-(i)}】\\ &= |Y|+1 \end{align}

      • (※)を繰り返し適用すれば  \langle Y \cup X \rangle = \langle Y \rangle
      • 補題-(iii)】より

\begin{align} \text{dim} \langle Y \cup X \rangle &= \text{dim} \langle Y \rangle \\ &= |Y|\end{align}

      • 一方、 \text{dim} \langle X \rangle \le \text{dim} \langle Y \cup X \rangle \quad \because補題-(ii)】
      •  \therefore |X| \le |Y| となり  |X| \gt |Y| に矛盾  \square
      • ok

 

グラフのマトロイド

【定義 2.2.1】
【定義 2.2.1】(接続行列)
 G: 無向グラフ とする。
 e \in E(G) に対し、 n 次元列ベクトル  B(G)_e の第  v \in V(G) 成分  B(G)_{ve} を次で定める:
\begin{align} B(G)_{ve} := \begin{cases} 0 &\text{if } v \notin e \\ 1 &\text{else if } v \in e \text{かつ} e \neq \{v, v\} \\ 2 &\text{otherwise} \end{cases} \end{align}
この列ベクトルを横に並べた  V(G) \times E(G) 行列  B(G) を、 G の接続行列という。

特に、 B(G) v 行目の総和は  v の次数を表します。(重要)

【定義 2.2.2】
【定義 2.2.2】(閉路マトロイド)
 G: 無向グラフ, B(G)_e: Gの接続行列の第e列ベクトル とする。
 \mathcal{I} \subseteq 2^{E(G)} を次で定める:
\begin{align} \mathrel{I} := \{X \subseteq E(G) \mathrel{|} \{ &B(G)_e \mathrel{|} e \in X\} が \\ &\mathbb{Z}_2 上の |V(G)| 次元ベクトル空間において線形独立 \} \end{align}このとき、 \mathcal{I} E(G) のマトロイドであり、Gの閉路マトロイドと呼ぶ。

これがマトロイドであることは、【命題 2.1.1】より明らかです。ありがとうベクトルマトロイド…!!

【命題 2.2.3】
【命題 2.2.3】(閉路マトロイドと森)
 G: 無向グラフ, \mathcal{I}: Gの閉路マトロイド とする。このとき
 X \in \mathcal{I} \Leftrightarrow G[X]は森

 G[X] X による  G の誘導部分グラフを表します。*6

証明:

  •  \Rightarrow の証明:(対偶を示す。 X \subseteq E(G) について、 G[X]: 閉路を持つ と仮定)
    • 閉路に含まれる辺集合を  C とおくと  \sum_{e \in C} B(G)_e = \boldsymbol{0}*7
    • よって、以下のように  \sum_{e \in X} \lambda_e B(G)_e = \boldsymbol{0} の非自明解  \lambda_e が得られる:
    • \begin{align} \lambda_e := \begin{cases} 1 & (\text{if } e \in C) \\ 0 & (\text{otherwise}) \end{cases} \end{align}
    •  \therefore X \notin \mathcal{I} ~~ \square
  •  \Leftarrow の証明:(対偶を示す。 X \notin \mathcal{I} と仮定。)
    •  X: 自己ループを持つ なら、明らかに  G[X] は森でない。以下、 X: 自己ループを持たない とする。
    •  \lambda_e \sum_{e \in X} \lambda_e B(G) = \boldsymbol{0} の非自明解とし、 Y := \{e \in X \mathrel{|} \lambda_e = 1\} とする。( Y \neq \emptyset
    • このとき、 Y は閉路を持つ。
      •  \because )
      •  \boldsymbol{0} = \sum_{e \in X} \lambda_e B(G)_e = \sum_{e \in Y} \lambda_e B(G)_e より、各  v \in V(G[Y]) について  v G[Y] における次数は偶数であり、 Y に誘導されていることに注意すると次数は2以上であるとわかる。
      • ここで、 G[Y] で最大のパス  P とその端点の1つ  w をとれば、 w の次数が2以上であることから  \exists u \in V(P) : \{w, u\} \in E(G[Y]) \setminus E(P) であり、この  P+\{w, u\} Y の持つ閉路である。 \square
    • したがって、 X は閉路を持つため、 G[X] は森でない。 \square

 

この命題から、森の辺集合全体は閉路マトロイドであることが得られます。また、閉路の辺集合は閉路マトロイドに対するサーキットであることもわかります。

特に、グラフが連結だとすると極大な森は全域木です。つまり、全域木は閉路マトロイドの基に対応します。

 

以上により、【定理 1.3.6】と合わせることでついに結論が得られます!

【定理 2.3.1】
【定理 2.3.1】(クラスカル法の正当性)
最小全域木問題に対し、クラスカル法は最適解を出力する。

 

ちなみにおまけとして、【定理 1.3.5】と合わせれば最大全域木問題も貪欲法で解けるという結論が得られます。知っていましたか?

最大全域木問題を使って解ける問題

 

おわりに

最後まで読んでいただいた方、お疲れさまでした!

本当は具体例を書いて定義の理解が正しいかなどを細かく確かめられるようにした方が絶対いいと思うのですが、めんどくさかったのでサボりました。気が向いたら追記するかもしれません。

 

参考文献

マトロイド - Wikipedia

離散最適化基礎論 (2015年度後学期) 組合せ最適化におけるマトロイドの役割

組合せ最適化 理論とアルゴリズム 第6版(B. コルテ, J. フィーゲン)

入門線形代数(三宅敏恒)

 

 

 

*1:ちなみに、この説明は雑過ぎて不正確です。

*2:マトロイドには複数の同値な定義があり、マトロイドの公理系と呼ばれています。特に、ここに書いた定義は「独立集合による定義」になります。他には「サーキットによる定義」、「階数関数による定義」など、色々あるようです。

*3:マトロイドの公理系の観点から、本来は 劣モジュラ性などの性質=階数関数による定義 が成立することを示しておくべきなのですが、クラスカル法の正当性の証明には使わないため省いております。

*4: Xの基を B_Xとして包含関係を見れば示せます。

*5:私の手元にある教科書では、入門線形代数-三宅敏恒 の 定理4.4.4 から得られます。

*6:「誘導部分グラフ」は通常頂点集合により誘導することがほとんどですが、ここでは辺集合で誘導しています。

*7: C: 自己ループ の場合も、法2で考えると  \boldsymbol{0} です。

永続セグ木と区間mexクエリの話

この記事でやること

  • 全永続セグ木のPython (PyPy3) による非再帰実装
  • 区間mexクエリのオフライン処理方法の解説

全永続セグ木って?

セグ木を全永続化したものです。

1点更新をする際にバージョンを名付けておき、以降の1点更新や区間演算を指定したバージョンのセグ木上で行うことができます。計算量はいずれも元のセグ木と変わらず  O(\text{log}N) 。すごい。*1

実現方法は、こちらのブログ の解説が非常に簡潔でわかりやすいと思いました。(丸投げ)

根を選ぶだけでバージョンが指定できるため、そこからは元のセグ木と同様に二分探索もできます。

全永続セグ木の実装例

PyPy提出を想定して非再帰で書いてます。木dpを非再帰で書く方法 に近いイメージで、dfsで通った頂点を \text{tps}に積んでおき、帰りがけは \text{tps}の逆順にまとめて処理しています。

class PersistentSegTree:
    class Node:
        def __init__(self, L, R) -> None:
            self.L = L
            self.R = R
            self.x = None
            self.left = None
            self.right = None
            self._calc = None
    def __init__(self, array, op, e, original_v=0) -> None:
        self.length = len(array)
        self.size = 1<<(self.length-1).bit_length()
        self.op = op
        self.e = e
        self.original_v = original_v
        self.version = {original_v: self.Node(0, self.size)}
        stack = [self.version[self.original_v]]
        tps = []
        while stack:
            node = stack.pop()
            if node.R-node.L > 1:
                tps.append(node)
                node.left = self.Node(node.L, (node.L+node.R)>>1)
                stack.append(node.left)
                node.right = self.Node((node.L+node.R)>>1, node.R)
                stack.append(node.right)
            else:
                node.x = array[node.L] if node.L < self.length else self.e
        for node in reversed(tps):
            node.x = self.op(node.left.x, node.right.x)
    def set(self, old_v, new_v, i, x):
        assert old_v in self.version and new_v not in self.version and 0 <= i < self.length
        self.version[new_v] = self.Node(0, self.size)
        old_node, new_node = self.version[old_v], self.version[new_v]
        tps = []
        while new_node.R-new_node.L > 1:
            tps.append(new_node)
            m = (new_node.L+new_node.R)>>1
            if new_node.L <= i < m:
                new_node.left = self.Node(new_node.L, m)
                new_node.right = old_node.right
                old_node, new_node = old_node.left, new_node.left
            else:
                new_node.left = old_node.left
                new_node.right = self.Node(m, new_node.R)
                old_node, new_node = old_node.right, new_node.right
        new_node.x = x
        for node in reversed(tps):
            node.x = self.op(node.left.x, node.right.x)
    def get(self, v, i):
        assert v in self.version and 0 <= i < self.length
        node = self.version[v]
        while node.R-node.L > 1:
            if node.L <= i < (node.L+node.R)>>1:
                node = node.left
            else:
                node = node.right
        return node.x
    def all_prod(self, v):
        assert v in self.version
        return self.version[v].x
    def prod(self, v, l, r):
        assert v in self.version
        if not 0 <= l < r <= self.length: return self.e
        stack = [(self.version[v], l, r)]
        tps = []
        while stack:
            node, l, r = stack.pop()
            m = (node.L+node.R)>>1
            if node.L == l and r == node.R:
                node._calc = node.x
            elif r <= m:
                tps.append((node, 0))
                stack.append((node.left, l, r))
            elif m <= l:
                tps.append((node, 1))
                stack.append((node.right, l, r))
            else:
                tps.append((node, 2))
                stack.append((node.left, l, m))
                stack.append((node.right, m, r))
        for node, look in reversed(tps):
            if look == 0:
                node._calc = node.left._calc
            elif look == 1:
                node._calc = node.right._calc
            else:
                node._calc = self.op(node.left._calc, node.right._calc)
        return self.version[v]._calc
    def max_right(self, v, l, f):
        assert v in self.version and 0 <= l < self.length and f(self.e)
        tps = []
        node = self.version[v]
        while node.L != l:
            if l < (node.L+node.R)>>1:
                tps.append(node.right)
                node = node.left
            else:
                node = node.right
        tps.append(node)
        x = self.e
        for node in reversed(tps):
            if not f(self.op(x, node.x)):
                while node.R-node.L > 1:
                    if f(self.op(x, node.left.x)):
                        x = self.op(x, node.left.x)
                        node = node.right
                    else:
                        node = node.left
                return node.L
            x = self.op(x, node.x)
        return self.length
    def min_left(self, v, r, f):
        assert v in self.version and 0 < r <= self.length and f(self.e)
        tps = []
        node = self.version[v]
        while node.R != r:
            if (self.L+self.R)>>1 < r:
                tps.append(node.left)
                node = node.right
            else:
                node = node.left
        tps.append(node)
        x = self.e
        for node in reversed(tps):
            if not f(self.op(node.x, x)):
                while node.R-node.L > 1:
                    if f(self.op(node.right.x, x)):
                        x = self.op(node.right.x, x)
                        node = node.left
                    else:
                        node = node.right
                return node.R
            x = self.op(node.x, x)
        return 0

 

区間mexクエリの処理方法

 \text{mex}とは、非負整数の集合を引数に取る関数で、「集合に含まれない最小の非負整数」を返します。

競プロでは、特にGrundy数の定義に登場することでよく知られています。

 

以下、長さ Nの非負整数列 Aに対し、 \text{mex}(A[l, r)) を求めるクエリを考えてみます。モノイドでないので、そのままセグ木に乗ったりはしません。悲しい。

オフラインの場合

 0 \le i \lt N について、 A[0, r)における i の最大のindex(ないなら-∞)を S_r[i] と表します。このとき、

 i \in A[l, r) \Leftrightarrow S_r [i] \ge l

であることに注意すると、 \text{min}(S_r [0, i)) \lt l を満たす最小の非負整数 i \text{mex}(A[l, r)) であるとわかります。

これは、 S_rをRmQセグ木で持つことで二分探索により求められる*2ので、クエリを rについて昇順に見ていきながら適切にセグ木を更新していくことで、全体 O((N+Q)\text{log}N)で解けました。

オンラインの場合

オフラインの場合における S_rを永続セグ木のバージョン rとするだけです。計算量は変わらず O(N+Q)\text{log}N

実装例

class RangeMexQuery:
    def __init__(self, array) -> None:
        self.pst = PersistentSegTree([-inf]*len(array), min, inf, 0)
        for r, a in enumerate(array):
            if a < len(array): self.pst.set(r, r+1, a, r)
    def prod(self, l, r):
        assert 0 <= l < r <= self.pst.length
        return self.pst.max_right(r, 0, lambda x: x >= l)

 

終わりに

そのうち『永続セグ木と区間K-thクエリの話』も書きたい…。

 

参考文献

熨斗袋さんのツイート

永続データ構造 - Wikipedia

永続セグメント木・永続遅延セグメント木 – 37zigenのHP

非再帰セグ木上の任意始点にぶたん - えびちゃんの日記

非再帰BFSでトポソから木DPをする - Qiita

Grundy数(Nim数, Nimber)の理論

*1:空間計算量は  O(N+Q\text{log}N) に悪化します。

*2:非再帰セグ木上の任意始点にぶたん - えびちゃんの日記 の表記を借りると、述語を p(x) := x \ge l とすればよいです。

彩色多項式の話【ABC294-Ex問題】

 

この記事でやること

ABC294-Ex問題 にてグラフ彩色の通り数に関する問題が出題され、その際彩色多項式なるものを勉強したので備忘録として残しておきます。

 

彩色多項式とは

グラフ  G k-彩色の通り数を  k に関する多項式で表したものです。

これを  P_{G}(k) と表すことにします。

 

木の彩色多項式

 n 頂点の木  T に対し、彩色多項式

 P_{T}(k) = k(k-1)^{n-1}

と表せます。これは、適当な根を取って自由に色を決め、そこから親と異なる色を順に塗っていけることからすぐにわかります。

 

彩色多項式の漸化式

彩色多項式  P_{G}(k) と辺  e \in V(G) について、以下の漸化式が成立します:

 P_{G}(k) = P_{G-e}(k) - P_{G \backslash e}(k)

 G-e G から  e を除いたグラフ、 G \backslash e G e について縮約したグラフ

 

証明

グラフ  G-e k-彩色の通り数を考える。 e = (u, v) としたとき、

  1.  u, v を異なる色に彩色する場合:
    •  G k-彩色と1対1対応するため、 P_{G}(k) 通り。
  2.  u, v を同じ色に彩色する場合:
    •  G \backslash e k-彩色と1対1対応するため、 P_{G \backslash e}(k) 通り。

以上より、 P_{G-e}(k) = P_{G}(k) + P_{G \backslash e}(k) であり、漸化式が得られる。 \square

 

単純連結グラフの彩色多項式

与えられたグラフに対し、閉路に含まれる1辺を  e として漸化式を適用することを木になるまで繰り返すことで得られます。*1縮約の際に多重辺ができる場合は、重複を除いて単純性が崩れないように注意しましょう。

 

ACコード

Submission #39887705 - AtCoder Beginner Contest 294

 

参考文献

https://ocw.hokudai.ac.jp/wp-content/uploads/2016/01/GraphTheory-2005-Note-10.pdf

*1:連結性を保つために閉路上の辺を選んでいます。また、この解法の計算量についてはあまりよくわかっていません。最悪ケースをかなり緩く見積もるモデルでの測定値が10^9程度の状態数となっており、実際はもっと少ないはずなので十分間に合うだろう推測しております。

Dinic法の計算量を真面目に調べた話

 

この記事でやること

  • Dinic法の仕組みの概略解説
  • 一般のグラフに対するDinic法の計算量と、実用上の高速化が期待できる基本的なポイントの考察
  • 二部グラフの最大マッチング問題でのDinic法の計算量の考察

 

はじめに(茶番)

AtCoder Beginner Contest 285 にて、二部グラフのマッチングに関する問題が出題されました。

atcoder.jp

自力では解法が思いつけなかったので(精進不足)解説を開いてみると、どうやら「2」と書かれた頂点を全て含むようなマッチングが存在するか?という問題に帰着されるらしい。なるほど確かに。

あとは最小流量制限付きのフローが存在するかを判定すればいいので、昔書いたDinic法のライブラリを引っ張り出してきて…グラフを構築して…ソイヤ!!

あああああ!!!!!

いやいや待て待て、さすがにちょっと適当過ぎたか…。よく考えてみればここらへんとかこんな感じに書けばいいか…あーここは書き直した方がいいなぁ…。こんな感じか?ダメか…こうか?こうすればいいのか?

ああああああああああああああ!!!!!!!!!!

というとても悲しい事件が起きたので、Dinic法と大真面目に格闘して知ったこと・考えたことを備忘録として残しておきます。

 

注意点

  • 計算量解析については参照元が英語なので、私のガバガバ英訳&補完となっていることをご了承ください。間違っていたらごめんなさい。(どなたか査読して頂けると非常に助かります。)
  • 一応Python3.8.3によるソースコードを掲載しておりますが、かなり汚いです。なので、具体的な実装例を参照したい方は他の方のソースコードを漁ってみることをお勧めします。

 

Dinic法の概要

全体像としては、下記の記事が非常にわかりやすいです。

Dinic法 - tkw’s diary

概略

  • Ford-Fulkerson法によると、「残余ネットワーク上で s-tパスを見つけて流す」を繰り返せば最大流が求まるので、「いい感じの順番で s-tパスを見つけて流す」ことでループ回数を抑えたい。
  • その手順は、
    1. 残余ネットワーク上でBFSにより各頂点の sからの(辺の重みを全て1とした場合での)最短距離 \text{dist}を記録する。
    2.  \text{cap} \gt 0 かつ  \text{dist} [ u ] \lt \text{dist} [ v ] 」を満たす辺 (u, v)のみを見るDFSで s-tパスを見つけ流せるだけ流すことを、 s-tパスが見つからなくなるまで繰り返す。
    3. BFSで s-tパスが見つからなくなるまで1, 2を繰り返す。

という感じです。

また、2のDFSを繰り返すパートでは配列

 \text{search_from}[v] := 頂点v から出る辺について, 探索中の辺の\text{index}

を管理することで、(BFS基準での同一パート内で)次回以降のDFSで既に探索済みの辺は参照しないようにします。具体的な処理としては、「DFSで戻ってきたとき、かつ頂点 t にまだ到達できていないとき」にインクリメントすればよいです。

 

実装例

(Python3.8.3)

from collections import deque

class Dinic:

    def __init__(self, n, inf=float("inf")): # n: 頂点数, inf: minの単位元、十分大の整数にした方が速いかも
        self.V = n
        self.inf = inf
        self.G = [[] for _ in range(n)]
        self.dist = []
        self.search_from = []
    
    def add_edge(self, from_, to, cap): # 容量capの辺(from_, to)をグラフに追加
        self.G[from_].append([to, cap, len(self.G[to])])
        self.G[to].append([from_, 0, len(self.G[from_])-1])

    def _bfs(self, s, t):
        self.dist = [-1]*self.V
        self.dist[s] = 0
        dq = deque([s])
        while dq:
            now = dq.pop()
            for edge in self.G[now]:
                next, cap, _ = edge
                if cap > 0 and self.dist[next] < 0:
                    self.dist[next] = self.dist[now]+1
                    if next == t: return True
                    dq.appendleft(next)
        return False

    def _dfs(self, s, t):
        stack = [t]
        goal = False
        flow = self.inf
        while stack:
            now = stack.pop()
            if now >= 0:
                if now == s:
                    goal = True
                    continue
                i = self.search_from[now]
                while i < len(self.G[now]):
                    self.search_from[now] = i
                    pre, _, rev = self.G[now][i]
                    _, cap, _ = self.G[pre][rev]
                    if cap > 0 and self.dist[now] > self.dist[pre] >= 0:
                        stack.append(~now)
                        stack.append(pre)
                        break
                    i += 1
            else:
                now = ~now
                if goal:
                    i = self.search_from[now]
                    pre, _, rev = self.G[now][i]
                    flow = min(flow, self.G[pre][rev][1])
                    continue
                else:
                    i = self.search_from[now]+1
                    while i < len(self.G[now]):
                        self.search_from[now] = i
                        pre, _, rev = self.G[now][i]
                        _, cap, _ = self.G[pre][rev]
                        if cap > 0 and self.dist[now] > self.dist[pre] >= 0:
                            stack.append(~now)
                            stack.append(pre)
                            break
                        i += 1
        if not goal: return 0
        now = t
        while now != s:
            i = self.search_from[now]
            pre, _, rev = self.G[now][i]
            self.G[now][i][1] += flow
            self.G[pre][rev][1] -= flow
            now = pre
        return flow

    def max_flow(self, s, t): # sからtへの最大流を求める
        flow = 0
        while True:
            if not self._bfs(s, t): return flow
            self.search_from = [0]*self.V
            f = self._dfs(s, t)
            while f > 0:
                flow += f
                f = self._dfs(s, t)

使用例↓

Submission #38339906 - AtCoder Beginner Contest 010

 

計算量評価

一般のグラフでのDinic法の計算量は  O(V^{2}E) となることが知られています:

  •  \because )
  •  \text{dist}[t]はBFS毎に狭義単調増加する(※証明は後述)ので、1, 2のループは高々 V-1回しか発生しません。
  • BFSは1回あたり  O(V+E)
  • DFSでは1回あたり、主に下記の2種類の(BFS基準の同一ループ内で)探索されなくなる辺が存在します:
    •  s-tパスに含まれず、行き止まりに到達する辺
    •  s-tパス上に含まれる辺のうち、 \text{cap} が最小の辺*1
  • よって、BFS基準の同一ループ内でのDFSによる計算量は「全ての辺が探索されなくなるまでの計算量」で評価することにします。
    • 前者の辺は一度使用すると \text{search_from} により以降は参照されなくなるため、1本あたりの計算量は O(1) とみなせます。*2
    • 後者の辺は1本あたりの計算量を O(\text{dist}[t]) = O(V) とみなせます。
  • 辺は E 本なので、合わせて  O(VE) となります。
  • 以上より、全体で  O(V^{2}E)

 

  • (※の証明)
    • 現在のBFSパートでの残余ネットワークにおいて、その前のBFSパートでの \text{dist}の振る舞いを考えます。
    • 現在の残余ネットワークの任意の辺 (i, j)で、 \text{dist}[j]-\text{dist}[i] \le 1が成立します。・・・①
      •  \because ) 
      •  \text{dist}[j]-\text{dist}[i] \gt 1 を満たす辺 (i, j)が存在すると仮定すると、BFSの性質よりその辺は前BFSパートでは存在しない、つまり、その後のDFSパートで逆辺にフローを流したことになります。
      • DFS時のルールにより \text{dist}[j]+1 = \text{dist}[i]
      •  \therefore \text{dist}[j] - \text{dist}[i] = -1 となり、仮定に矛盾  \square
    • 現在のBFSパートで見つけた s-tパスにおいて、 \text{dist}[s] = 0であることと①より、このパスの長さの下限は \text{dist}[t]とわかります。
    • また、このパスの中には \text{dist}[j] - \text{dist}[i] \lt 1を満たす辺 (i, j)が少なくとも1つ含まれます。
      •  \because )  
      • 含まれない、すなわち任意の辺 (i, j) \text{dist}[j] - \text{dist}[i] = 1 と仮定すると、これは前DFSパートで見つかるはずの増加パスとなり現在のBFSパートに移っていることに矛盾  \square
    • よって、このパスの長さの下限は \text{dist}[t]+1となります。  \square

辺1本あたりの \text{dist}の増分は高々+1、かつ少なくとも1本は増分+0以下で、 \text{dist}[t]まで増やす必要がある

実用上の高速化ポイント

(言語仕様などによる)実装上の最適化ではなく、アルゴリズム的に実用上の高速化が期待できるポイントを挙げます。

BFSは頂点t に到達した時点で打ち切る

dijkstra等の最短路問題でもよく言われるやつです。DFSで参照する辺の条件から、頂点 t 以降に訪れる頂点は必ず行き止まりとなります。*3

実は頂点s からDFSするより頂点t からDFSした方が速い(らしい?)

 s からのDFSでは「(BFS基準の同一ループ内で)初めから行き止まりの頂点」にも訪れてしまいますが、 t からのDFSの場合はそのような頂点には訪れることがありません。*4元から行き止まりが少ない場合は逆辺を参照する回数が増える分ほんの少し遅くなるかもしれませんが、トータルで見れば概ね実行時間は改善されます。(要検証)

Submission #38178001 - AtCoder Beginner Contest 285

 s からDFS, TLE)

Submission #38177995 - AtCoder Beginner Contest 285

 t からDFS, AC [2754ms])

 

二部マッチングにおけるDinic法の計算量

二部グラフにおける最大マッチングは、下図のようなグラフを構築することで最大流問題に帰着できます。

辺の \text{cap}は全て1

このグラフにおいては、Dinic法の計算量が O( (V+E)\sqrt{V})となることが知られています。*5正確には、下記の条件を満たすグラフではこの計算量でDinic法が動作します:

  •  s, t を除くすべての頂点について、流せる流量が高々1。すなわち、「出る辺の \text{cap}の和」と「入る辺の \text{cap}の和」のいずれか一方は高々1。
  • この条件を満たすグラフは以下の2つの性質を持ちます:
    1. 増加パスに流せる流量は明らかに1。*6
    2. 増加パスにフローを流しても、残余ネットワークは常に条件を満たす。
      •  \because ) 
      • 残余ネットワークが条件を満たすとし、新たに増加パスにフローを流したとします。この流量は1です。
      • パス上の s, t を除く各頂点について、流れた流量1のフローにより、入る辺の \text{cap}1つが出る辺(逆辺)へと変わり、出る辺の \text{cap}1つが入る辺(逆辺)へと変わります。
      • このとき、 \text{cap}の収支は変わっていないため、条件を満たしたままとなります。 \square
  • この条件を満たすグラフをunit networkと言うらしい…?*7

以下、計算量の証明:

  •  \because )
  • DFSパートで s-t パスを見つけた際、先述の条件よりパスに含まれる辺のうち少なくとも \text{dist}[t]/2本の \text{cap}は1なので、それらは全て以降のDFSで探索されなくなります。
  • 言い換えると、DFSパートで見つかる s-t パスの総距離は高々 2E となります。つまり、DFSパート1回あたりの計算量は(行き止まりの分も合わせて) O(E)となります。

 s-t パス1本で探索されなくなる辺が \text{dist}[t]/2本となる例
  • 更に、BFSパートは O(\sqrt{V}) 回しか発生しないことが示せます:
    •  \because )  
    • 初めて \text{dist}[t] \sqrt{V}以上となったときのBFSパートを考えます。
    •  \text{dist}の狭義単調増加性より、それまでにBFSパートは O(\sqrt{V})回しか発生しません。
    • このBFSパートに至るまでに流した流量を f、最終的に得られる最大フローを F とします。このとき、現在の残余ネットワークにおける最大流問題の解、すなわち追加で流せる流量は F-f となります。
    • 更に、先述のグラフの性質から s-t の点素なパスが F-f 本あると言い換えられます。
    • これらパスの長さは \text{dist}[t] なので、パスに使用する頂点数に関する不等式  (F-f)(\text{dist}[t]-1) \le V が得られます。
    • ここで、 \text{dist}[t] \ge \sqrt{V} より F-f = O(\sqrt{V})
    • ゆえに、BFSパート1回につき少なくともフローを1は流せることからこれ以降のBFSパートは O(\sqrt{V}) 回しか発生ないとわかります。
    • 以上を合わせると、BFSパートは O(\sqrt{V})回しか発生しません。 \square
  • 以上より、全体で O( (V+E)\sqrt{V}) \square 

 

おわりに

二部マッチングに対するDinic法は O(E\sqrt{V})で動作するという話は結構色々なところで見かけましたが、その証明が書かれた文献は全然見つからなかったのでかなり苦労しました。誰か書いておいてくださいよー。

ちなみにDinic法の特殊なグラフでの計算量ってなんか色々あるみたいですね。疲れたので調べるつもりはありませんが…(今のところは)

 

参考文献

Dinic法 - tkw’s diary

Dinic 法とその時間計算量 - みさわめも

https://people.orie.cornell.edu/dpw/orie633/LectureNotes/lecture9.pdf

Maximum flow - Dinic's algorithm - Algorithms for Competitive Programming

*1:フローを流すことでこの辺が探索されなくなる代わりに逆辺が生じますが、 \text{dist}が反転しているためDFSには影響しません。

*2:行き止まりになる前に「 s-tパス上に含まれる辺のうち、 \text{cap} が最小でない辺」として複数回使われる可能性はありますが、これは後者の辺の計算量に含まれます。

*3:BFSの性質上、頂点 t 以降に訪れる任意の頂点の \text{dist} \text{dist}[t] 以上であるため。

*4:一度以上DFSすることにより新たに生じる行き止まりには到達する可能性があります。

*5: V=O(E)とみなして O(E\sqrt{V})と書く場合も多いです。

*6:厳密には s, t \text{cap} \gt 1 の辺で接続しているという例外はありますが、計算量への影響は明らかにないので無視します。

*7:参考文献[4] ではそう呼ばれていましたが、広く使われている表現なのかはわかりません。