はじめに

こんにちは、ユーゴです。今回は、量子アニーリングの手法で遊んでいたので、そちらの紹介をします。

やること

人間関係のスコアを元に、班の組み合わせが最適になるグルーピングを、量子アニーリングの手法で求める。

リポジトリ

先に、リポジトリを共有しておきます。
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. コードに書き起こす

ライブラリをインポート

先に以下のライブラリを導入しましょう。

requirements.txt

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など幅広い技術の紹介をしております。
お役に立てましたら、いいね, フォロー等よろしくお願いいします！！