0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

【量子アニーリング】班決め最適化アルゴリズムを作ってみた

Posted at

はじめに

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

やること

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

リポジトリ

先に、リポジトリを共有しておきます。
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が理論値(だったはず)です。
score.png

私のコードだと、制約は100%守られていましたが、やや理論値がとりにくかったです。もしかしたら、「たまに制約が破られる」くらいのペナルティの方が、理論値がとりやすいかもしれません。とはいえ、結局古典アルゴリズムのご機嫌取りなので、アニーリングマシンだと上手くいくのかもしれません。
ぜひ、ペナルティ係数を調整してみてください。

ちなみに、OpenJijの方が、体感理論値が出やすかったです。ありがとうございます。

リポジトリ

最初にも添付しましたが、ぜひご興味ある方はリポジトリもご覧ください。
mu5dvlp/group-optimization

ちなみにリポジトリでは、Jupyterで書いてます。

余談

どうやら世の中には「量子アニーリング最高!!」派閥と、「量子ゲート方式以外認めん!!」派閥が存在するようですが、私は中道です。とにかく新しい技術で、価値を生み出したい。
今回のトピックは量子アニーリングですが、Qiskitでゲート方式の操作を通じて画像処理なんかもしてます。映像表現について、流体など現状のGPUに限界を感じ始めているためです。

【資格】
IBM Certified Associate Developer - Quantum Computation using Qiskit v0.2X

【記事】
量子コンピュータで画像処理してみた

終わりに

いかがだったでしょうか。今回は、量子アニーリングを用いて、班決め最適化アルゴリズムを作ってみた紹介をしました。修学旅行の班決めで不服な結果となったことを思い出して、作ったという背景がありました。
このように、量子コンピューティングのテーマから、Unity, AWS, GASなど幅広い技術の紹介をしております。
お役に立てましたら、いいね, フォロー等よろしくお願いいします!!

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?