A. Lucas, Front. Phys. (2014) 掲載例題の実装と解説 ー 有向帰還頂点集合問題¶
本サンプルコードでは、論文 A. Lucas, "Ising formulations of many NP problems", Front. Phys. (2014) で紹介されている『有向帰還頂点集合問題』に Fixstars Amplify を用いて取り組みます。同論文に紹介されている他の NP 完全・NP 困難な問題も以下で解説しています（カッコ内は論文内で問題に対応する節番号）。
- グラフの分割問題（2.2節）
- 最大クリーク問題（2.3節）
- 厳密被覆問題（4.1節）
- 集合パッキング問題（4.2節）
- 最小頂点被覆問題（4.3節）
- 充足可能性問題（SAT）（4.4節）
- 最小極大マッチング問題（4.5節）
- グラフ彩色問題（6.1節）
- クリーク被覆問題（6.2節）
- 整数長ジョブスケジューリング問題（6.3節）
- ハミルトン閉路問題（7.1節）
- 有向帰還頂点集合問題（8.3節）
- 最小帰還辺集合問題（8.5節）
- グラフ同型性判定問題（9節）
有向帰還頂点集合問題¶
有向グラフ $G$ が与えられたとき、$G$ の頂点の部分集合 $F$ であって、どの $G$ の閉路も $F$ の頂点を $1$ つ以上通るものを 有向帰還頂点集合 といいます。 言い換えると、$G$ の $F$ に含まれない頂点を出発して、$F$ に含まれない頂点のみを通ってもとの頂点に戻ってくることができないとき、$F$ は $G$ の有向帰還頂点集合になります。
例えば、下図のグラフに対して、オレンジ色で示した頂点の部分集合は、有向帰還頂点集合の 1 つとなります。
有向帰還頂点集合問題 とは、有向グラフ $G$ に対して、$G$ の帰還頂点集合のうち最小の要素数のものを求める問題です。本サンプルプログラムの定式化は A. Lucas, Front. Phys. (2014) の 8.3 節のものに沿って行います。
問題の作成¶
本サンプルプログラムで取り組む帰還頂点集合問題として、NetworkX を用いて有向グラフ $G$ を作成します。
import networkx as nx
N = 8 # グラフ頂点の数
G = nx.DiGraph()
G.add_nodes_from(range(N))
edge_list = [
(0, 1),
(1, 2),
(2, 3),
(3, 4),
(4, 5),
(5, 6),
(6, 7),
(7, 0),
(4, 2),
(7, 1),
(7, 5),
]
G.add_edges_from(edge_list)
node_labels = {
0: "a",
1: "b",
2: "c",
3: "d",
4: "e",
5: "f",
6: "g",
7: "h",
}
pos = nx.circular_layout(G) # レイアウトを保存しておく
nx.draw_networkx(G, node_size=600, font_color="w", labels=node_labels, pos=pos)
作成したグラフ $G$ から頂点 $e$, $f$ を取り除いたグラフは閉路を持たないことは簡単に分かります。つまり、$\{e, f\}$ は $G$ の帰還頂点集合です。 また、閉路 $c\rightarrow d\rightarrow e\rightarrow c$ と閉路 $f\rightarrow g\rightarrow h\rightarrow f$ は共通部分を持たないので、$G$ の帰還頂点集合の要素数は $2$ 以上です。 したがって、この問題に関しては、$G$ の帰還頂点集合の要素数の最小値は $2$ となります。
定式化¶
方針¶
以下、$G$ の頂点の個数を $N$ とします。
まず、$N$ 個のバイナリ変数 $y$ を各頂点と対応付けて、それぞれの頂点が帰還頂点集合 $F$ に含まれるかどうかを表すことにします。$F$ に含まれるなら $0$ で含まれないなら $1$ です。
次に、本問題の言い換えである「$F$ に含まれない頂点からなる $G$ の部分グラフ $H$ が閉路を持たない」という条件は、さらに「$H$ の頂点に番号をうまく付けると、 $H$ のすべての辺が番号の小さな頂点から大きな頂点に向かって出ているようにできる」という条件に言い換えることができます（証明：$\Rightarrow$ は簡単、$\Leftarrow$ はトポロジカルソート）。 この番号付けは $N\times N$ のバイナリ変数テーブル $x$ を用いて、頂点 $v$ の番号が $i$ であるとき $v$ 行 $i$ 列のバイナリ変数を $1$ とすることで表現できます。
たとえば、上で作成した問題は以下のようなグラフになっています。
このグラフの各頂点に以下のように色と番号を付けると、オレンジ色の2点 $e$, $f$ は帰還頂点集合となっていて、青い頂点同士を結ぶ辺は番号の小さな頂点から大きな頂点に向かって出ています。
このような帰還頂点集合の選び方および番号の付け方に対応するバイナリ変数 $y$, $x$ は以下の表のようになります。ただし、帰還頂点集合に含まれる頂点と対応する $x$ の行はすべて $0$ とすることにします。
| a | b | c | d | e | f | g | h | |
|---|---|---|---|---|---|---|---|---|
| $y$ | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 |
| $$x$$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| a | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| b | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
| c | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| d | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
| e | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| f | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| g | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| h | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
目的関数¶
帰還頂点集合の要素数ができるだけ少なくなればよいので、目的関数は $\displaystyle -\sum_{v=0}^{N-1} y_v$ となります。$y_v$ は 頂点 $v$ が帰還頂点集合に含まれるなら $0$、そうでないなら $1$ となることに注意してください。
制約条件¶
$y$ および $x$ が帰還頂点集合を表すためには、以下が必要です。
- 条件1： $F$ に含まれない頂点には番号が $1$ つ付けられている。つまり、 $x$ の $v$ 行目は、$v$ が帰還頂点集合に含まれるならばすべて $0$ であり、そうでないならば $1$ つだけが $1$ である。
- 条件2： $G$ の辺 $u\rightarrow v$ について、$u$ と $v$ がともに帰還頂点集合に含まれないならば、$u$ の番号は $v$ の番号よりも小さい。つまり、このとき自然数 $i \leq j$ に対して、$x_{u, j}$ と $x_{v, i}$ が両方 $1$ であってはならない（注意：$x_{u, \underline{i}}$ と $x_{v, \underline{j}}$ は両方 $1$ になり得る）。
条件 1 は
$$ \sum_{i=0}^{N-1} x_{v, i} = y_v \quad \text{for} \quad v \in \{0, 1, \ldots, N-1\} $$
と表すことができます。
また、条件 1 より、$u$, $v$ のどちらかが帰還頂点集合に含まれるならば $x_{u, j}$ と $x_{v, i}$ がともに $1$ になることはないので、条件 2 のうち「$u$ と $v$ がともに帰還頂点集合に含まれない」という条件は自然と考慮されます。したがって、条件 2 は
$$ x_{u, j} x_{v, i} = 0 \quad \text{for} \quad (u, v) \in E, \ 0 \leq i \leq j < N $$
と表すことができます。
逆に、バイナリ変数 $y$, $x$ が条件 1, 2 をみたしているとき、対応する $y$ が $y=0$ となっている頂点の集合は帰還頂点集合となるので、これらを制約条件として与えればよいです。
実装¶
上で作成した問題と定式化を使って、実際に問題を解いてみましょう。
まず、Fixstars Amplify SDK の BinarySymbolGenerator を使い、バイナリ変数 $y$ と $x$ を作成します。
from amplify import VariableGenerator
gen = VariableGenerator()
y = gen.array("Binary", shape=(N,))
x = gen.array("Binary", shape=(N, N))
次に、目的関数 $\displaystyle -\sum_v y_v$ を作成します。
cost = -y.sum()
条件 1 に対応する制約条件を作成します。条件 1 は、$F$ に含まれない各頂点に番号が付けられていることを表し、これは、前述の通り、$x$ の各行の和が $y$ の各要素に等しいと言い換えることができます。
まず、$x$ の各行の和と $y$ の各要素との差を表す一次元配列を作成します。次に、 equal_to 関数の axis
パラメータに空のタプルを指定することで、この一次元配列の各要素がすべて 0 に等しいという制約条件を一度に生成できます。
from amplify import equal_to
diff = x.sum(axis=1) - y
constraint1 = equal_to(diff, 0, axis=())
条件 2 に対応する制約条件を作成します。条件 2 は、$x_{u, j} x_{v, i} = 0 \ \bigl((u, v) \in E, \ 0 \leq i \leq j < N\bigr)$ という制約です。
from amplify import sum as amplify_sum
constraint2 = amplify_sum(
equal_to(x[u, j] * x[v, i], 0)
for u, v in G.edges
for i in range(N)
for j in range(i, N)
)
目的関数と制約条件を足し合わせて組合せ最適化モデルを作成します。
制約条件は目的関数に対するペナルティ関数としてイジングマシンに与えられるため、制約条件に対する重みとして、目的関数の取り得る値とおよそ同等の値またはそれより少々大きめの値を推定して決定します。今回は、制約の重みを $2$ とします。
model = cost + (constraint1 + constraint2) * 2
クライアントを設定し、Fixstars Amplify Annealing Engine (AE) で実行します。Amplify SDK
は制約条件をみたす解を自動でフィルターするので、result が空でなければ解が見つかったと分かります。
from amplify import FixstarsClient, solve
from datetime import timedelta
client = FixstarsClient()
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # ローカル環境等で使用する場合は、Fixstars Amplify AE のアクセストークンを入力してください。
client.parameters.timeout = timedelta(milliseconds=1000) # タイムアウトは 1000 ms
# 求解を実行
result = solve(model, client)
if len(result) == 0:
print("解が見つかりませんでした")
else:
print("解が見つかりました")
最後に、結果を可視化します。
import numpy as np
values = result.best.values
y_values = y.evaluate(values)
x_values = x.evaluate(values)
numbering = {v: "" for v in G.nodes}
numbering.update(dict(np.argwhere(x_values == 1)))
colors = ["C0" if v == 1 else "C1" for v in y_values]
nx.draw_networkx(
G, node_size=600, node_color=colors, font_color="w", labels=numbering, pos=pos
)
