A. Lucas, Front. Phys. (2014) 掲載例題の実装と解説 ー 集合パッキング問題¶
本サンプルコードでは、論文 A. Lucas, 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節）
集合パッキング問題とは¶
集合 $S$ があり、$S$ の部分集合 $T_0, T_1, \ldots, T_{N-1}$ が与えられているとします。 $T_0, T_1, \dots, T_{N-1}$ の中からそれらが共通部分を持たないようにいくつかを選び、選んだ部分集合の要素数の総和ができるだけ大きくなるようにする問題を 集合パッキング問題 といいます。
たとえば、下図のように $S = \{1, 2, 3, 4, 5, 6, 7, 8, 9\}$ で、部分集合として $T_0 = \{1, 2, 3, 6, 9\}$、 $T_1 = \{1, 2, 5, 8\}$、 $T_2 = \{4, 7\}$、 $T_3 = \{4, 5\}$、 $T_4 = \{6, 9\}$ の場合を考えます。この時、$T_1$、$T_2$、$T_4$ を選ぶと要素数の総和が $8$ となり、最大となります。一方で、たとえば $T_0$ と $T_1$ は共通部分をもつので、両方を同時に選ぶことはできません。
ここでは、Fixstars Amplify を用いて集合パッキング問題を解くプログラムを作成します。定式化は A. Lucas, Front. Phys. (2014) の 4.2 節のものに沿って行います。
問題の作成¶
例題として、上に挙げた問題をコードで表現しておきます。
S = [1, 2, 3, 4, 5, 6, 7, 8, 9] # 集合 S
T = [[1, 2, 3, 6, 9], [1, 2, 5, 8], [4, 7], [4, 5], [6, 9]] # Sの複数の部分集合
定式化¶
決定変数¶
$N$ 個のバイナリ変数 $q$ を $T_0, T_1, \ldots, T_{N-1}$ と対応付けて、対応する部分集合 $T_i$ を選ぶかどうかを表すことにします。$T_i$ を選ぶなら $q_i$ は $1$ で、選ばないなら $0$ です。
たとえば、$T_1$, $T_2$, $T_4$ の 3 つの部分集合を選ぶときは、決定変数 $q$ は以下のようになります。
| 部分集合 | $$T_0$$ | $$T_1$$ | $$T_2$$ | $$T_3$$ | $$T_4$$ |
|---|---|---|---|---|---|
| $$q$$ | 0 | 1 | 1 | 0 | 1 |
目的関数¶
選んだ部分集合の要素数の和をできるだけ大きくしたいので、目的関数は
$$ -\sum_{i = 0}^{N - 1} q_i \cdot (\# T_i) $$
で表すことができます。ここで、$\# T_i$ は $T_i$ の要素数です。マイナスの符号がついているのは、最大化問題を最小化問題に変換するためです。
制約条件¶
「選んだ部分集合は共通部分 (overlap) を持たない」という条件を $q$ に課す必要があります。 これは、
$$ q_i q_j = 0 \quad \text{if} \quad T_i\ \text{and} \ T_j \ \text{overlap} $$
という式で書くことができます。
実装¶
上で作成した問題と定式化を使って、実際に問題を解いてみましょう。最初に、Fixstars Amplify SDK の BinarySymbolGenerator
を使って部分集合の数だけバイナリ変数 $q$ を作成します。
from amplify import VariableGenerator
N = len(T)
gen = VariableGenerator()
q = gen.array("Binary", N)
次に、目的関数を作成します。前述の通り、目的関数は、$-\sum_{i = 0}^{N - 1} q_i \cdot (\# T_i)$ で表されますが、次のように実装できます。
import numpy as np
subset_lengths = np.array([len(t) for t in T]) # 各 T_i の要素数を表す配列
cost = -(q * subset_lengths).sum()
続いて、制約条件を作成します。制約条件は、$q_i q_j = 0 \ \bigl(\text{if} \:\: T_i \:\: \text{and} \:\: T_j \:\: \text{overlap}\bigr)$ です。
from amplify import equal_to, sum as amplify_sum
import itertools
def overlap(t_i, t_j):
return len(set(t_i) & set(t_j)) > 0
constraints = amplify_sum(
equal_to(q[i] * q[j], 0)
for i, j in itertools.combinations(range(N), 2)
if overlap(T[i], T[j])
)
作成した目的関数と制約条件をまとめて組合せ最適化モデルを構築します。
制約条件は目的関数に対するペナルティ関数としてイジングマシンに与えられるため、制約の重みとして、目的関数の取り得る値と同等の値またはそれより少々大きめの値を推定して決定します。今回は、制約の重みを、$\max(\#T_i)$ とします。
model = cost + np.max(subset_lengths) * constraints
クライアントを設定し、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("解が見つかりました")
最後に、結果を可視化します。余裕があれば、集合 $S$ やその部分集合 $T_i$ を変更して、求解してみましょう。
print(
f"要素数の和：{int(-result.best.objective)}"
) # 得られた最適解に対応する目的関数の値を表示
values = q.evaluate(result.best.values)
for i in np.where(values == 1)[0]:
print(f"T{i} : {T[i]}")
