A. Lucas, Front. Phys. (2014) 掲載例題の実装と解説ー最小極大マッチング問題¶

本サンプルコードでは、論文 A. Lucas, "Ising formulations of many NP problems", Front. Phys. (2014) で紹介されている『最小極大マッチング問題』に Fixstars Amplify を用いて取り組みます。同論文に紹介されている他の NP 完全・NP 困難な問題も以下で解説しています（カッコ内は論文内で問題に対応する節番号）。

最小極大マッチング問題¶

グラフ $G$ に対して、$G$ の辺の部分集合 $D$ が以下をみたすとき、$D$ を 極大マッチング といいます。

$D$ に含まれる辺同士は、隣接しない。
$D$ に含まれない辺は、必ず $D$ のいずれかの辺と隣接している。

たとえば、以下の図のオレンジ色の辺は、極大マッチングとなっています。オレンジ色の辺同士がつながっていないことと、黒い辺を 1 本でもオレンジ色に塗ったとするとオレンジ色の辺がつながってしまい、極大マッチングではなくなることを確認してください。

最小極大マッチング問題は、与えられたグラフに対して、そのグラフの極大マッチングのうち要素数が最小となるものを求める問題です。

本サンプルプログラムでは、Fixstars Amplify を用いて最小極大マッチングを求めるプログラムを作成します。定式化は A. Lucas, Front. Phys. (2014) の 4.5 節のものに沿って行います。

問題の作成¶

まず、問題として、NetworkX を用いて適当なグラフ $G$ を作成します。

In [ ]:

import networkx as nx
    
N = 6  # グラフの頂点の数
    
G = nx.Graph()
G.add_nodes_from(range(N))
edge_list = [
    (0, 1),
    (0, 5),
    (1, 2),
    (1, 5),
    (2, 3),
    (2, 4),
    (3, 4),
    (4, 5),
]  # 頂点と頂点を結ぶ辺の定義
G.add_edges_from(edge_list)
    
pos = nx.circular_layout(G)
    
nx.draw_networkx(G, node_size=600, font_color="w", pos=pos)

定式化¶

以下、$G$ の頂点の数を $N$、辺の数を $M$ とします。

決定変数¶

$M$ 個のバイナリ変数 $q$ を $G$ の各辺と対応付けて、それぞれの辺が極大マッチング $D$ に含まれるかどうかを表すことにします。 $D$ に含まれるなら $1$, 含まれないなら $0$ です。

たとえば、以下のような極大マッチングに対しては、バイナリ変数 $q$ は下の表のようになります。

辺 $(u, v)$	$$(0, 1)$$	$$(0, 5)$$	$$(1, 2)$$	$$(1, 5)$$	$$(2, 3)$$	$$(2, 4)$$	$$(3, 4)$$	$$(4, 5)$$
$$q$$	0	0	0	1	0	1	0	0

目的関数¶

$D$ の要素数ができるだけ小さくなるようにすればよいので、$ \displaystyle \sum_{i = 0}^{M - 1} q_i$ を最小化します。

制約条件¶

先の説明の通り、$D$ が極大マッチングであるとは、以下が満たされている、ということです。

条件 1 ： $D$ に含まれる辺同士は隣接しない。
条件 2 ： $D$ に含まれない辺は、必ず $D$ のいずれかの辺と隣接している。

これらの条件を言い換えて $q$ で表すことを考えます。

まず、条件 1 は「隣接する $2$ 本の辺がともに $D$ に含まれることはない」と言い換えられます。これは、

$$ q_{v, u} q_{v, w} = 0 \quad \text{for} \quad (v, u), (v, w) \in E $$

と表すことができます。ただし、辺 $(u, v)$ に対応するバイナリ変数配列 $q$ の要素を $q_{u, v}$ と書いています。また、$E$ は $G$ の辺集合を表します。

次に、条件 2 は、「$G$ のすべての辺は、必ず $D$ のいずれかの辺と隣接している」と言い換えられます。これをさらに、「$G$ のどの辺 $(u, v)$ に対しても、$u$ と $v$ のどちらかは $D$ のいずれかの辺の端点となっている」と言い換えます。ある頂点 $v$ が $D$ のいずれかの辺の端点となっているかどうかは、$v$ から出るすべての辺について、対応するバイナリ変数の総和が $1$ であるか $0$ であるかを見れば判定できるので、条件 2 は

$$ (1 - \sum_{(v, x) \in E} q_{v, x}) (1 - \sum_{(u, y) \in E} q_{u, y}) = 0 \quad \text{for} \quad (u, v)\in E $$

で表すことができます。

実装¶

上で作成した問題と定式化を使って、実際に問題を解いてみましょう。最初に、Fixstars Amplify SDK の BinarySymbolGenerator を使って $M$ 個のバイナリ変数 $q$ を作成します。

In [ ]:

from amplify import VariableGenerator
    
M = len(G.edges)
    
gen = VariableGenerator()
q = gen.array("Binary", M)

前述の定式化に沿って目的関数を作成します。目的関数は、極大マッチング $D$ の要素数と等しく、$\displaystyle \sum_{i = 0}^{M - 1} q_i$ で表されます。

In [ ]:

cost = q.sum()

制約条件を作成する準備として、$G$ の各頂点 $v$ に対し、$v$ から出ている辺のインデックスのリストを作っておきます。以下のコードにおいて、edge_indices_list[v] は、ノード v から出ている辺のインデックスのリストとなります。

In [ ]:

edge_indices_list = [[] for _ in range(N)]
for i, (u, v) in enumerate(G.edges):
    edge_indices_list[u].append(i)
    edge_indices_list[v].append(i)

条件 1 に対応する制約条件を作成します。条件 1 は、極大マッチング $D$ に含まれる 2 辺が隣接しない、つまり、隣接する $2$ 本の辺がともに $D$ に含まれないことを意味し、$q_{v, u} q_{v, w} = 0 \ \bigl((v, u), (v, w) \in E\bigr)$ で表されます。

In [ ]:

from itertools import combinations
from amplify import equal_to, sum as amplify_sum
    
constraint1 = amplify_sum(
    equal_to(q[i] * q[j], 0)
    for v in G.nodes
    for i, j in combinations(edge_indices_list[v], 2)
)

条件 2 に対応する制約条件を作成します。条件 2 は、すべての辺が $D$ のいずれかの辺と隣接していることを意味し、 $\displaystyle(1 - \sum_{(v, x) \in E} q_{v, x}) (1 - \sum_{(u, y) \in E} q_{u, y}) = 0 \ \bigl((u, v)\in E\bigr)$ で表されます。

In [ ]:

constraint2 = amplify_sum(
    equal_to(
        (1 - amplify_sum([q[i] for i in edge_indices_list[u]]))
        * (1 - amplify_sum([q[i] for i in edge_indices_list[v]])),
        0,
    )
    for (u, v) in G.edges
)

作成した目的関数と制約条件をまとめて、組合せ最適化モデルを構築します。

今回は必要ありませんが、目的関数と制約条件の両方が存在する場合、通常、制約条件に重みを掛けた方がよい場合があります。これは、制約条件は目的関数に対するペナルティ関数としてイジングマシンに与えられるためです。基本的な考え方として、目的関数の取り得る値と同等の値またはそれより少々大きめの値を推定して決定します。

In [ ]:

model = cost + constraint1 + constraint2

クライアントを設定し、Fixstars Amplify Annealing Engine (AE) で実行します。Amplify SDK は制約条件をみたす解を自動でフィルターするので、result が空でなければ、制約条件をみたす解が見つかったと分かります。

In [ ]:

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("解が見つかりました。")

最後に、結果を可視化します。

In [ ]:

values = q.evaluate(result.best.values)
    
colors = ["k" if i == 0 else "C1" for i in values]
width = [1.0 if i == 0 else 2.0 for i in values]
    
nx.draw_networkx(
    G, node_size=600, font_color="w", edge_color=colors, width=width, pos=pos
)

A. Lucas, Front. Phys. (2014) 掲載例題の実装と解説 ー 最小極大マッチング問題¶

最小極大マッチング問題¶

問題の作成¶

定式化¶

決定変数¶

目的関数¶

制約条件¶

実装¶

A. Lucas, Front. Phys. (2014) 掲載例題の実装と解説ー最小極大マッチング問題¶