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$ のすべての頂点を一回ずつ通ってもとに戻ってくるような閉路をハミルトン閉路といいます。 一般に、グラフのサイズが大きいとき、グラフにハミルトン経路が存在するかどうかを現実的な時間で判定することは困難です。
ここでは、Fixstars Amplify を用いて、ハミルトン閉路を探索するプログラムを作成します。本問題は A. Lucas, Front. Phys. (2014) の 7.1 節に対応します。
問題の作成¶
まず、NetworkX を用いて本サンプルプログラムで取り扱うグラフ $G$ を作成します。頂点の数は $N$ 個です。
import networkx as nx
N = 5 # グラフの頂点の数
G = nx.Graph()
G.add_nodes_from(range(N))
edge_list = [(0, 1), (0, 2), (1, 2), (1, 3), (1, 4), (2, 4), (3, 4)]
pos = nx.circular_layout(G) # レイアウトを保存しておく
G.add_edges_from(edge_list)
nx.draw_networkx(G, node_size=600, font_color="w", pos=pos)
定式化¶
決定変数¶
$N\times N$ 個のバイナリ決定変数 $q$ を用意し、どの頂点を何番目に通るかを表します。つまり、バイナリ決定変数のある成分 $ q_{k, i}$ は、頂点 $i$ を $k$ 番目に通る ($=1$) か否 ($=0$)か、と表します。例えば、バイナリ変数が以下のようになっているときは、閉路 $0 \rightarrow 1 \rightarrow 3 \rightarrow 4 \rightarrow 2 \rightarrow 0$ に対応します。
| 順番 \ 頂点番号 | 0 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|---|
| 1 番目 | 1 | 0 | 0 | 0 | 0 |
| 2 番目 | 0 | 1 | 0 | 0 | 0 |
| 3 番目 | 0 | 0 | 0 | 1 | 0 |
| 4 番目 | 0 | 0 | 0 | 0 | 1 |
| 5 番目 | 0 | 0 | 1 | 0 | 0 |
目的関数¶
ハミルトン閉路問題は条件をみたすものを見つける問題なので、目的関数は $0$（無し）となります。
制約条件¶
$q$ がハミルトン閉路を表すためには、以下が必要です。
-
$k$ 番目に通る頂点は必ず $1$ つである必要があります。これは、バイナリ変数 $q$ の各行に $1$ つだけ $1$ があると言い換えることができます。
-
各頂点はちょうど $1$ 回通る必要があります。これは、バイナリ変数 $q$ の各列に $1$ つだけ $1$ があると言い換えることができます。
-
辺が張られていない頂点間での移動はできません。つまり、頂点 $i$ と頂点 $j$ の間に辺が張られていないとき、$q_{k, i}$ と $q_{k+1, j}$ がともに $1$ であってはいけません。
条件 1～3 を数式で書き下すと、それぞれ以下のようになります。
\begin{align*} \sum_{i=0}^{N-1} q_{k, i} = 1 & \quad \text{for} \quad k \in \{0, 1, \ldots, N-1\} \\ \sum_{k=0}^{N-1} q_{k, i} = 1 & \quad \text{for} \quad i \in \{0, 1, \ldots, N-1\} \\ q_{k, i}q_{k+1, j} = 0 & \quad \text{for} \quad k \in \{0, 1, \ldots, N-1\}, (i, j) \notin E \end{align*}
ここで、$E$ は $G$ の辺集合を表します。
また、バイナリ変数 $q$ が条件 1~3 をすべてみたすとき、$q$ は $G$ のハミルトン閉路に対応します。
実装¶
上で作成した問題と定式化を使って、実際に問題を解いてみましょう。まず、Fixstars Amplify SDK の BinarySymbolGenerator を使って
$N\times N$ 個のバイナリ変数 $q$ を作成します。
from amplify import VariableGenerator
gen = VariableGenerator()
q = gen.array("Binary", shape=(N, N))
次に、条件 1 と 2 に対応する制約条件を作成します。これらは、$q$ のそれぞれの行と列にひとつだけ $1$ があるという条件でしたので、one_hot
を使って書くことができます。axis パラメータに 1 を指定すると二次元配列の各行に対する one-hot 制約を一度に生成でき、0 を指定すると各列に対する one-hot
制約を一度に生成できます。
from amplify import one_hot
row_constraints = one_hot(q, axis=1)
col_constraints = one_hot(q, axis=0)
構築したそれぞれの制約条件を表示し、正しく各行及び列に one_hot 条件が課されているか確認します。
row_constraints
col_constraints
次に、条件 3 に対応する制約条件を作成します。条件 3 は $q_{k, i}q_{k+1, j} = 0$ ($i$ と $j$ は辺で結ばれていない2頂点) という条件でした。 注意点として、$k=N-1$ のとき、 $q_{k+1, j}$ は $q_{0, j}$ を意味する必要があることに注意します。
from amplify import equal_to, sum as amplify_sum
edge_constraints = amplify_sum(
equal_to(q[k, i] * q[(k + 1) % N, j], 0) + equal_to(q[k, j] * q[(k + 1) % N, i], 0)
for (i, j) in nx.non_edges(G)
for k in range(N)
)
以上で必要な制約条件が揃いました。最後に、これらをまとめて組合せ最適化モデルを作成します。
from amplify import Model
model = Model(row_constraints + col_constraints + edge_constraints)
クライアントを設定して、Fixstars Amplify Annealing Engine で実行します。
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)
ハミルトン経路が見つかったかどうかを確認します。Amplify SDK は制約条件をみたす解を自動でフィルターするので、result
が空でなければハミルトン経路が見つかったと分かります。
if len(result) == 0:
print("ハミルトン経路が見つかりませんでした。")
else:
print("ハミルトン経路が見つかりました。")
最後に、結果を可視化します。ハミルトン経路をオレンジで示します。
import numpy as np
# デフォルトの辺のアトリビュートを設定
for edge in G.edges.values():
edge["color"] = "k"
edge["width"] = 1.0
# ハミルトン経路に含まれる辺のアトリビュートを設定
values = q.evaluate(result.best.values)
route = np.where(values == 1)[1]
for i, j in zip(route, np.roll(route, -1)):
G.edges[i, j]["color"] = "C1"
G.edges[i, j]["width"] = 2.0
# 描画
edge_color = [edge["color"] for edge in G.edges.values()]
edge_width = [edge["width"] for edge in G.edges.values()]
nx.draw_networkx(
G, node_size=600, font_color="w", pos=pos, edge_color=edge_color, width=edge_width
)
