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節）
最小帰還辺集合問題¶
有向グラフ GG が与えられたとき、GG の辺の部分集合 FF であって、どの GG の閉路も FF の辺を 11 つ以上含むものを 帰還辺集合 といいます。これを言い換えると、GG の FF に含まれない辺のみを通る閉路が存在しないとき、FF は GG の帰還辺集合になります。
例えば、下図のグラフにおいて、オレンジ色で示した辺は、帰還辺集合となります。
最小帰還辺集合問題 とは、有向グラフ GG に対して、GG の帰還辺集合のうち最小の要素数のものを求める問題です。
ここでは、Fixstars Amplify を用いて、GG の最小帰還辺集合を求めるプログラムを作成します。本サンプルプログラムの定式化は A. Lucas, Front. Phys. (2014) の 8.5 節のものに沿って行います。
問題の作成¶
本サンプルプログラムで取り組む最小帰還辺集合問題として、NetworkX を用いて有向グラフ GG を作成します。
import networkx as nx
N = 8 # グラフの頂点数
G = nx.DiGraph()
G.add_nodes_from(range(N))
edge_list = [
(0, 1),
(0, 6),
(1, 2),
(2, 3),
(3, 1),
(3, 4),
(3, 5),
(4, 5),
(5, 6),
(6, 7),
(7, 0),
]
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)
作成したグラフ GG から辺 b→cb→c 及び a→ga→g を取り除いたグラフは閉路を持たないことは簡単に分かります。つまり、{b→c{b→c, a→g}a→g} は GG の帰還辺集合です。
また、閉路 b→c→d→bb→c→d→b と閉路 a→g→h→aa→g→h→a は共通部分を持たないので、GG の帰還辺集合の要素数は必ず 22 以上です。
したがって、この問題に関しては、GG の帰還辺集合の要素数の最小値は 22 となります。
定式化¶
以下、GG の 頂点の数を NN、辺の数を MM とします。
方針¶
まず、MM 個のバイナリ変数 yy を各辺と対応付けて、それぞれの辺が帰還辺集合 FF に含まれるかどうかを表すことにします。FF に含まれるなら 00 で含まれないなら 11 です。
次に、「GG の FF に含まれない辺のみを通る閉路が存在しない」という条件は、 「GG の頂点に番号をうまく付けると、 GG の FF に含まれないすべての辺が番号の小さな頂点から大きな頂点に向かって出ているようにできる」という条件に言い換えることができます (証明：⇒⇒ は簡単、⇐⇐ はトポロジカルソート)。
この番号付けは N×NN×N のバイナリ変数テーブル xx を用いて、頂点 vv の番号を ii とするとき vv 行 ii 列のバイナリ変数を 11 とすることで表現できます。
上で作成したグラフを用いて、変数の対応付けの例を説明します。上で作成したグラフは以下のようになっています。
このグラフに対して、以下のように頂点に番号を付けると、オレンジ色の2辺は帰還辺集合となっていて（番号が減少する方向に辺が出ている為）、黒い辺はすべて番号の小さな頂点から大きな頂点に向かって出ています。
このような帰還辺集合の選び方および番号の付け方に対応するバイナリ変数 yy, xx は以下の表のようになります。
| 辺 | a→ba→b | a→ga→g | b→cb→c | c→dc→d | d→bd→b | d→ed→e | d→fd→f | e→fe→f | f→gf→g | g→hg→h | h→ah→a |
|---|---|---|---|---|---|---|---|---|---|---|---|
| y | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| x | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| a | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
| b | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |
| c | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| d | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| e | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| f | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
| g | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| h | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
目的関数¶
帰還辺集合の要素数ができるだけ少なくなればよいので、目的関数は −M−1∑e=0ye となります。ye は 辺 e が帰還辺集合 F に含まれるなら 0、そうでないなら 1 となることに注意してください。
制約条件¶
y および x が帰還辺集合を表すためには、以下が必要です。
- 条件 1： G の各頂点には、0 以上 N 未満の番号がが 1 つ付けられている。 つまり、x の各行には、ちょうど 1 つだけ 1 がある。
- 条件 2： G の辺 u→v について、u→v が帰還辺集合 F に含まれないならば、u の番号は v の番号よりも小さい。
条件 1 は
N−1∑i=0xv,i=1forv∈{0,1,…,N−1}
と表すことができます。
また、条件 2 は、「辺 u→v が F に含まれないならば、自然数 i≤j に対して、xu,j_ と xv,i_ が両方 1 であってはならない」と言い換えられるので、
yu→vxu,jxv,i=0for(u,v)∈E, 0≤i≤j<N
と表すことができます。ただし、E は G の辺集合を表し、yu→v は辺 u→v と対応する y の要素です。
条件 2 を表す式は 3 次式であり、イジングマシンで解くためには補助変数を使って 2 次式に変換する必要があります。 Fixstars Amplify SDK はこの変換を自動で行う機能を提供しています。以下では、Fixstars Amplify SDK の次数下げ機能を使う方法と、条件 2 をさらに言い換えて手動で 2 次式に落とす方法の 2 種類を解説します。
実装 (その1: Amplify の次数下げ機能を使用する)¶
上で作成した問題と定式化を使って、実際に問題を解いてみましょう。まずは、Fixstars Amplify の次数下げ機能を使用する方法です。
最初に、VariableGenerator を用いてバイナリ変数 y と x を作成します。
from amplify import VariableGenerator
gen = VariableGenerator()
M = len(G.edges) # 辺の数
y = gen.array("Binary", shape=(M,))
x = gen.array("Binary", shape=(N, N)) # N は頂点の数
次に、目的関数 −∑eye を作成します。
cost = -y.sum()
条件 1 に対応する制約条件を作成します。条件 1 は、x の各行に対する one-hot 制約です。axis パラメータに 1 を指定すると二次元配列の各行に対する
one-hot 制約を一度に生成できます。
from amplify import one_hot
constraint1 = one_hot(x, axis=1)
条件 2 に対応する制約条件を作成します。条件 2 は yu→vxu,jxv,i=0 ((u,v)∈E, 0≤i≤j<N) という制約です。
from amplify import equal_to, sum as amplify_sum
constraint2 = amplify_sum(
equal_to(y[e] * x[u, j] * x[v, i], 0)
for e, (u, v) in enumerate(G.edges)
for i in range(N)
for j in range(i, N)
)
最後に、目的関数と制約条件をまとめて組合せ最適化モデルを作成します。
制約条件は目的関数に対するペナルティ関数としてイジングマシンに与えられるため、制約条件への重みとして、目的関数の取り得る値とおよそ同等の値またはそれより少々大きめの値を推定して決定します。今回は、制約の重みを 2 とします。
penalty_multiplier = 2
model = cost + penalty_multiplier * (constraint1 + constraint2)
クライアントを設定し、Fixstars Amplify Annealing Engine (AE) で実行します。
次数下げを行う際に発行する補助変数の数を削減するため、次数下げアルゴリズムを Substitute
に指定します。このアルゴリズムは、共通の 2 次項 qiqj を因数に持つ 3 次以上の項が多くある場合に有効です。
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,
quadratization_method="Substitute", # 次数下げに使用するアルゴリズム
)
if len(result) == 0:
print("解が見つかりませんでした")
else:
print("解が見つかりました")
最後に求解結果を可視化します。先ほど確認したように、G の帰還辺集合の要素数の最小値は 2 なので、オレンジ色の辺が 2 個であれば最適解が見つかったことになります。
import numpy as np
# バイナリ変数をデコード
values = result.best.values
y_values = y.evaluate(values)
x_values = x.evaluate(values)
# ノードの番号を表示
numbering = dict(np.argwhere(x_values == 1))
# 各辺が F に含まれているかを表示
edge_colors = ["C1" if e == 0 else "k" for e in y_values]
edge_width = [2.0 if e == 0 else 1.0 for e in y_values]
# 描画
nx.draw_networkx(
G,
node_size=600,
font_color="w",
labels=numbering,
edge_color=edge_colors,
width=edge_width,
pos=pos,
)
以上で、最小帰還辺集合を求めるプログラムが作成できました。
次に、同じ問題を Amplify SDK の次数下げ機能を使わずに定式化して解く方法を紹介します。
定式化 (その2: 定式化が 2 次になるようにする)¶
上で行った定式化では、条件 2 が 3 次式になってしまいました。ここでは、バイナリ変数を新たに追加することで条件 2 を 2 次式で表現することを考えます。
前述の通り、条件 2 は、以下のような条件です。
G の辺 u→v について、u→v が帰還辺集合 F に含まれないならば、u の番号は v の番号よりも小さい。
定式化の方針¶
バイナリ変数 y, x については、上で定義した通りとします。
もし頂点 u の番号が i であると分かっているとすると、「u の番号が v の番号よりも小さい」という制約は、
∑j>ixv,j=1
という 1 次式で表せます。したがって、各辺について、帰還辺集合 F に含まれているかどうかと始点の番号を 1 次式で取得できれば、この式と OR をとることで条件 2 を表せそうです。
そこで、M×N のバイナリ変数テーブル z を用意し、F に含まれているかどうかと始点の番号を表すことにします。ここで、M は G の辺数で N は G の頂点数です。z の u→v に対応する行は、辺 u→v が F に含まれている場合はすべて 0 であり、そうでない場合は、 u の番号を i として、i 列目のみが 1 となります。
たとえば、以下の帰還辺の選び方 / 番号の付け方に対しては、z は下の表のようになります。
| z | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| a→b | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
| a→g | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| b→c | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| c→d | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| d→b | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| d→e | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| d→f | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| e→f | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| f→g | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
| g→h | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
| h→a | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
条件 2 の 2 次の定式化¶
バイナリ変数 z を用いて条件 2 を表すためには、z が以下をみたす必要があります。以下、z の 辺 u→v に対応する行を zu→v と書くことにします。
-
条件 2-1：z の各行は、対応する辺が帰還辺集合 F に含まれているかどうかを表す。つまり、z の各行は、対応する辺が F に含まれるならばすべて 0 であり、そうでないならば 1 つだけが 1 である。
-
条件 2-2：z の各行は、対応する辺が帰還辺集合 F に含まれないならば、その辺の始点の番号を表す。つまり、zu→v,i=1 は、 u→v が F に含まれず、かつ u の番号が i であることを意味する。
-
条件 2-3：辺 u→v が F に含まれず、頂点 u の番号が i であるならば、頂点 v の番号は i より大きい。
条件 2-1 は、
N−1∑i=0ze,i=yefore∈E
で表せます。ye は、辺 e が F に含まれるならば 0 となるバイナリ変数であることを思い出してください。
条件 2-2 については、「zu→v,i=1 ならば u の番号は i である」という条件を課せば十分です。「zu→v,i=1 ならば u→v が F に含まれない」ことは条件 2-1 より明らかなためです。したがって、条件 2-2 は
zu→v,i(1−xu,i)=0for(u→v)∈E, i∈{0,1,…,N−1}
で表せます。
条件 2-2 より、条件 2-3 の仮定は zu→v,i=1 と同値です。また、先述の通り「頂点 v の番号が i より大きい」という条件は ∑j>ixv,j=1 で表されるので、条件 2-3 は
zu→v,i(1−∑j>ixv,j)=0for(u→v)∈E, i∈{0,1,…,N−1}
で表せます。
これで制約条件 2-1 ~ 2-3 が定式化できました。これらがすべてみたされていれば制約条件 2 をみたすことは簡単に分かります。
from amplify import VariableGenerator
from amplify import one_hot
gen = VariableGenerator()
M = len(G.edges) # number of edges
y = gen.array("Binary", shape=(M,))
x = gen.array("Binary", shape=(N, N)) # N is number of nodes
cost = -y.sum()
constraint1 = one_hot(x, axis=1)
条件 2 を実装するために、バイナリ変数 z を定義します。
z = gen.array("Binary", shape=(M, N))
条件 2-1 ： ∑N−1i=0ze,i=ye (e∈E) を作成します。equal_to 関数の
axis パラメータに空のタプルを指定することで、第一引数のそれぞれの要素に対する制約条件を一度に生成できます。
from amplify import equal_to
constraint2_1 = equal_to(z.sum(axis=1) - y, 0, axis=())
条件 2-2 ： zu→v,i(1−xu,i)=0 ((u→v)∈E, i∈{0,1,…,N−1}) を作成します。
from amplify import sum
constraint2_2 = sum(
equal_to(z[e, i] * (1 - x[u, i]), 0)
for e, (u, v) in enumerate(G.edges)
for i in range(N)
)
条件 2-3 ：zu→v,i(1−∑j>ixv,j)=0 ((u→v)∈E, i∈{0,1,…,N−1}) を作成します。
条件 1 が成り立っているという条件の下で、条件 2-3 の左辺は最小値 0 をとるので、条件 2-3 の左辺はペナルティ関数として使用することができます。Constraint
クラスのコンストラクタに penalty を手動で指定することで、制約条件を構築します（詳細はこちら）。
一方で、equal_to 関数を用いると、ペナルティ関数を内部的に生成する際に左辺が 2 乗されて 4
次式になってしまいますので、今回の制約条件式においては equal_to 関数の利用を避けます。
from amplify import Constraint, sum
constraint2_3 = sum(
Constraint(
z[e, i] * (1 - x[v, i + 1 :].sum()),
eq=0,
penalty=z[e, i] * (1 - x[v, i + 1 :].sum()),
)
for e, (u, v) in enumerate(G.edges)
for i in range(N)
)
目的関数と制約条件をまとめて組合せ最適化モデルを作成します。条件 1 がみたされていないと constraint2_3 を充足しようとするポテンシャルが働かないので、条件 1
の重みを大きくしておきます。
model = cost + 2 * constraint1 + constraint2_1 + constraint2_2 + constraint2_3
Fixstars Amplify Annealing Engine (AE) で実行し、結果を可視化します。
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 = dict(np.argwhere(x_values == 1))
edge_colors = ["C1" if e == 0 else "k" for e in y_values]
edge_width = [2.0 if e == 0 else 1.0 for e in y_values]
nx.draw_networkx(
G,
node_size=600,
font_color="w",
labels=numbering,
edge_color=edge_colors,
width=edge_width,
pos=pos,
)