はじめに
こんにちは、ユーゴです。今回は、量子アニーリングの手法で遊んでいたので、そちらの紹介をします。
やること
人間関係のスコアを元に、班の組み合わせが最適になるグルーピングを、量子アニーリングの手法で求める。
リポジトリ
先に、リポジトリを共有しておきます。
mu5dvlp/group-optimization
手順
0. 定義
まず、必要な要素を列挙します。
- グループの集合：$G$
- 人間の集合：$V$
- グループ内の人数：$n$
- 人間関係の値：$w$
1. ハミルトニアンの立式
まずは、ハミルトニアンを作成します。
ソフト制約
前提知識として、量子アニーリングはソフト制約による計算となります。
どういうことかというと、「グループ内の人数は5人以下にしたい！」と言っても、この制約を絶対に守らせるのは困難です。
なので、制約を破ると、ペナルティを与えるという方法で、制約を破る選択肢を回避させていきます。
さて、今回守りたい制約は以下2つとなります。
- グループ内の人数はn人
- 1人あたりの所属グループは1つ
制約項の雛形として、以下のような式がよく使われます。
$$
P \left( \sum x - n \right)^2
$$
$P$：ペナルティ係数
$x$：0or1のバイナリ値
$n$：総和の目標値
これを元に立式すると、
グループ内人数：
$$
P_n \sum_{g\in G} \left( \sum_{v\in V} x_{v,g} - n \right) ^2
$$
グループ所属数
$$
P_g \sum_{v\in V} \left( \sum_{g\in G} x_{v,g} - 1 \right) ^2
$$
(余計なお世話かもですが) $\left( \sum_{v\in V} (x_{v,g} - n) \right) ^2$ではありません。最後に$n$引いてます。私はこれの理解が遅れて解釈に悩んでいました。
目的関数
どシンプルです。
AさんとBさんが同じ班にいるときは、積が1になります。それ以外は0です。
以下のように立式されます。
$$
-\sum_{g\in G}\sum_{(i,j)\in V} x_{i,g}x_{j,g}w_{i,j}
$$
マイナスがついているのは、エネルギーの最小化問題と解くためです。人間関係値が高いほどエネルギーが低くなれば良い訳です。
まとめると
以下のようになります。
$$
H = P_g \sum_{v\in V} \left( \sum_{g\in G} x_{v,g}-1 \right)^2 + P_n \sum_{g\in G} \left( \sum_{v\in V} x_{v,g}-n \right)^2 - \sum_{g\in G}\sum_{(i,j)\in V} x_{i,g}x_{j,g}w_{i,j}
$$
2. コードに書き起こす
ライブラリをインポート
先に以下のライブラリを導入しましょう。
openjij
dwave-ocean-sdk
numpy
networkx
pyqubo
matplotlib
matplotlib-inline
特に、今回のアニーリングはD-WaveとOpenJijを使ってシミュレートします。
D-Waveが2025年の1月から、有料化してしまったようです。本物のアニーリングマシンを使った評価はできませんでした😭
ハミルトニアンのコード
さて、ハミルトニアンが立式できました。
次にやるべきは、PyQUBOでハミルトニアンを表記します。
まず、パラメータを用意します：
# 人数
n_verts = 9
# グループ数
n_groups = 3
# 人間関係値の対応表
w = np.array(
[
[0, 10, 1, 1, 1, 1, 10, 1, 1],
[10, 0, 1, 1, 1, 1, 10, 1, 1],
[1, 1, 0, 10, 1, 1, 1, 10, 1],
[1, 1, 10, 0, 1, 1, 1, 10, 1],
[1, 1, 1, 1, 0, 10, 1, 1, 10],
[1, 1, 1, 1, 10, 0, 1, 1, 10],
[10, 10, 1, 1, 1, 1, 0, 1, 1],
[1, 1, 10, 10, 1, 1, 1, 0, 1],
[1, 1, 1, 1, 10, 10, 1, 1, 0],
]
)
自動計算できるパラメータは以下：
# 自動でグループ内の人数を計算
n_group_verts = n_verts // n_groups
# ペナルティ係数。今回は簡易的に自動計算しているが、必要に応じて調整するとよい
P_g = sum([w[i, j] for i, j in itertools.combinations(range(n_verts), 2)]) / 2
print(f"Groups penalty: {P_g}")
P_n = sum([w[i, j] for i, j in itertools.combinations(range(n_verts), 2)]) / 2
print(f"Node in group penalty: {P_n}")
また、バイナリ値の参照＆人間関係値の参照のため、indexの配列を作っておきます。割と可読性のためだけの存在です：
# 0〜人数-1の配列を作る。あとでバイナリ値の参照index、人間関係値の参照indexに使う。
V = [i for i in range(n_verts)]
# 0〜グループ数-1の配列を作る。あとでバイナリ値の参照index、人間関係値の参照indexなどに使う。
G = [i for i in range(n_groups)]
最後に、ハミルトニアンを立式します：
# バイナリ値の集合
x = {(v, g): Binary(f"{v},{g}") for v in V for g in G}
# ハミルトニアン
H = (
P_g
* Constraint(
sum((sum(x[v, g] for g in G) - 1) ** 2 for v in V), label="group_constraint"
)
+ P_n
* Constraint(
sum((sum(x[v, g] for v in V) - n_group_verts) ** 2 for g in G),
label="node_constraint",
)
- sum(
w[i][j] * x[i, g] * x[j, g]
for g in G
for (i, j) in itertools.combinations(V, 2)
)
)
先ほどの式
$$
H = P_g \sum_{v\in V} \left( \sum_{g\in G} x_{v,g}-1 \right)^2 + P_n \sum_{g\in G} \left( \sum_{v\in V} x_{v,g}-n \right)^2 - \sum_{g\in G}\sum_{(i,j)\in V} x_{i,g}x_{j,g}w_{i,j}
$$
を、そのままコードに書き起こしただけなのがお分かりいただけるでしょうか？
全体
ハミルトニアン含め、以下のようになります：
import openjij as oj
from pyqubo import Binary, Constraint
import neal
import numpy as np
import networkx as nx
import itertools
import matplotlib.pyplot as plt
n_verts = 9
n_groups = 3
w = np.array(
[
[0, 10, 1, 1, 1, 1, 10, 1, 1],
[10, 0, 1, 1, 1, 1, 10, 1, 1],
[1, 1, 0, 10, 1, 1, 1, 10, 1],
[1, 1, 10, 0, 1, 1, 1, 10, 1],
[1, 1, 1, 1, 0, 10, 1, 1, 10],
[1, 1, 1, 1, 10, 0, 1, 1, 10],
[10, 10, 1, 1, 1, 1, 0, 1, 1],
[1, 1, 10, 10, 1, 1, 1, 0, 1],
[1, 1, 1, 1, 10, 10, 1, 1, 0],
]
)
n_group_verts = n_verts // n_groups
P_g = sum([w[i, j] for i, j in itertools.combinations(range(n_verts), 2)]) / 2
print(f"Groups penalty: {P_g}")
P_n = sum([w[i, j] for i, j in itertools.combinations(range(n_verts), 2)]) / 2
print(f"Node in group penalty: {P_n}")
V = [i for i in range(n_verts)]
G = [i for i in range(n_groups)]
x = {(v, g): Binary(f"{v},{g}") for v in V for g in G}
H = (
P_g
* Constraint(
sum((sum(x[v, g] for g in G) - 1) ** 2 for v in V), label="group_constraint"
)
+ P_n
* Constraint(
sum((sum(x[v, g] for v in V) - n_group_verts) ** 2 for g in G),
label="node_constraint",
)
- sum(
w[i][j] * x[i, g] * x[j, g]
for g in G
for (i, j) in itertools.combinations(V, 2)
)
)
model = H.compile()
qubo, _ = model.to_qubo()
def get_score(x):
group_penalty = P_g * sum((sum(x[v, g] for g in G) - 1) ** 2 for v in V)
node_penalty = P_n * sum((sum(x[v, g] for v in V) - n_group_verts) ** 2 for g in G)
objective_value = sum(
w[i][j] * x[i, g] * x[j, g]
for g in G
for (i, j) in itertools.combinations(V, 2)
)
return group_penalty + node_penalty - objective_value
labels_row = [f"v{i}" for i in V]
labels_col = [f"g{i}" for i in G]
def show_result(rows):
_, ax = plt.subplots(figsize=(len(G), len(V) / 3))
ax.axis("tight")
ax.axis("off")
table_data = []
for i, row in enumerate(rows):
table_data.append([labels_row[i]] + list(row))
ax.table(
cellText=table_data, colLabels=[""] + labels_col, cellLoc="center", loc="center"
)
plt.show()
def plot_graph(rows):
graph = nx.Graph()
groups = {g: [] for g in G}
for v in V:
row = rows[v]
for g in G:
if row[g] == 1:
groups[g].append(v)
for g in G:
graph.add_nodes_from(groups[g], group=g)
for i, j in itertools.combinations(V, 2):
if rows[i][g] == 1 and rows[j][g] == 1:
graph.add_edge(i, j)
pos = nx.spring_layout(graph, k=1)
nx.draw(graph, pos=pos, with_labels=True, node_color="lightgreen")
plt.show()
sampler_dwave = neal.SimulatedAnnealingSampler()
response_dwave = sampler_dwave.sample_qubo(qubo)
decoded_samples = model.decode_sample(response_dwave.first.sample, vartype="BINARY")
solution_dwave = {key: val for key, val in decoded_samples.sample.items() if val == 1}
rows = np.zeros((len(V), len(G)), dtype=int)
for key, val in solution_dwave.items():
v, g = key.split(",")
rows[int(v), int(g)] = val
score = get_score(rows)
print(f"Score: {score}")
show_result(rows)
plot_graph(rows)
sampler_openjij = oj.SASampler()
response_openjij = sampler_openjij.sample_qubo(qubo)
decoded_samples = model.decode_sample(response_openjij.first.sample, vartype="BINARY")
solution = {key: val for key, val in decoded_samples.sample.items() if val == 1}
rows = np.zeros((len(V), len(G)), dtype=int)
for key, val in solution.items():
v, g = key.split(",")
rows[int(v), int(g)] = val
score = get_score(rows)
print(f"Score: {score}")
show_result(rows)
plot_graph(rows)
3. 結果を確認
スコアが表示されます。
コードの例だと、-90.0が理論値(だったはず)です。
私のコードだと、制約は100%守られていましたが、やや理論値がとりにくかったです。もしかしたら、「たまに制約が破られる」くらいのペナルティの方が、理論値がとりやすいかもしれません。とはいえ、結局古典アルゴリズムのご機嫌取りなので、アニーリングマシンだと上手くいくのかもしれません。
ぜひ、ペナルティ係数を調整してみてください。
ちなみに、OpenJijの方が、体感理論値が出やすかったです。ありがとうございます。
リポジトリ
最初にも添付しましたが、ぜひご興味ある方はリポジトリもご覧ください。
mu5dvlp/group-optimization
ちなみにリポジトリでは、Jupyterで書いてます。
余談
どうやら世の中には「量子アニーリング最高！！」派閥と、「量子ゲート方式以外認めん！！」派閥が存在するようですが、私は中道です。とにかく新しい技術で、価値を生み出したい。
今回のトピックは量子アニーリングですが、Qiskitでゲート方式の操作を通じて画像処理なんかもしてます。映像表現について、流体など現状のGPUに限界を感じ始めているためです。
【資格】
IBM Certified Associate Developer - Quantum Computation using Qiskit v0.2X
【記事】
量子コンピュータで画像処理してみた
終わりに
いかがだったでしょうか。今回は、量子アニーリングを用いて、班決め最適化アルゴリズムを作ってみた紹介をしました。修学旅行の班決めで不服な結果となったことを思い出して、作ったという背景がありました。
このように、量子コンピューティングのテーマから、Unity, AWS, GASなど幅広い技術の紹介をしております。
お役に立てましたら、いいね, フォロー等よろしくお願いいします！！