はじめに
量子アニーリングの座学的な学習を進めておりましたが実際のプログラミングや定式化、パラメータの調整については手を動かさないと把握できないと思い、その辺の学習をかねて実務的な内容の問題設定を行い実際に解いてみる試みをしました。
問題設定は、chatgptを用いて行っております。
今回は最大カット問題を題材としてプロジェクトのアサイン最適化を解こうと思います。
問題設定
メンバー10人を想定として、以下のような設定で5人ずつで相性の悪いメンバーの切り離したアサインを目指します。
- メンバー10人
- 非相性スコアの設定：
- レビュー渋滞度：PRのリードタイム、差し戻し回数
- 依存ブロック：タスクの待ち時間、ブロック発生履歴
- コードオーナー競合：同一クリティカル領域の同時変更頻度
- 時間帯/シフトの重なり：コアタイム一致率、タイムゾーン差
- 会議負荷：同席時間の長さ・未収束率
- 障害連鎖：ある組合せで障害が増えやすいか
非相性スコアは効率的に業務ができていなそうな条件を想定して0~1に正規化した数値として設置しています。(ランダムに数値は決めています)
また、本来は10人分の完全グラフに対して問題を解いたほうが良いかと思いますが、簡単のために相性スコアに閾値をかけ疎なグラフとなるようにしております。
この設定によりグラフを疎にしつつ、よりカットしないといけない部分が切れるイメージです。
この記事での記号と定義
グラフ：無向・重み付き
$$
\boxed{G=(V,E,w)},\quad V={A,\ldots,J},\quad w:E\to[0,1]
$$
- 表にあるペアのみを $E$ に採用（＝疎グラフ）。表にないペアは $w_{ij}=0$ とみなす。
- 重みは対称 $w_{ij}=w_{ji}$、自己ループは $w_{ii}=0$。
割当変数：各ノード $i\in V$ に $x_i\in{0,1}$（2グループのいずれか）。
以降、A〜J はノードID、$i,j$ はそれらを指す添字、$x_i$ はノード $i$ の所属（0/1）、$w_{ij}$ は非相性スコアとします。
共通データ（10人・非相性スコア）
メンバー
A~Jの10人
非相性スコア（同じチームにしないように離す）
（表にあるペアのみエッジとして採用。他は 0 とする）
| $i$ | $j$ | $w_{ij}$ |
|---|---|---|
| G | I | 0.91 |
| B | D | 0.87 |
| C | E | 0.87 |
| B | E | 0.86 |
| A | H | 0.86 |
| E | F | 0.77 |
| C | D | 0.72 |
| A | J | 0.70 |
| C | I | 0.70 |
| B | J | 0.68 |
| C | J | 0.64 |
| B | F | 0.61 |
| E | J | 0.60 |
| C | F | 0.60 |
| B | G | 0.59 |
| B | H | 0.59 |
| D | J | 0.58 |
上の表をグラフにすると以下のような形になっています。
問題を解いてみる
今回の事例は最大カット問題となりますので、d-waveのMax Cutの例を参考にしていきます。
目的項
バイナリ変数 $x_i\in{0,1}$ はノード $i$（メンバー）の所属グループを表します。
D-Wave の Max-Cut に従い、別グループなら 1 になる指示関数（XOR）で数えます：
$$
x_i \oplus x_j = x_i + x_j - 2x_i x_j
$$
これに重み $w_{ij}$ を掛けて目的項は
$$
\max\ H_{\mathrm{obj}}=
\sum_{(i,j)\in E} w_{ij}\big(x_i+x_j-2x_i x_j\big)
$$
制約項
グループをちょうど半分（5/5）にしたいので、$\sum_i x_i$ が 5 からずれた分だけ罰するペナルティを加えます：
$$
H_{\mathrm{bal}} = \rho\Big(\sum_i x_i - 5\Big)^2
$$
本稿では再現性のため $\rho=2.0$ とします。
QUBO式
OpenJij（BQM）は最小化問題として解くため、目的を反転して
$$
\min\ H(x)= -H_{\mathrm{obj}} + H_{\mathrm{bal}}=
-\sum_{(i,j)\in E} w_{ij}(x_i+x_j-2x_i x_j)
+\rho\Big(\sum_i x_i - 5\Big)^2
$$
とします。
プログラム
今回はPyQUBOとOpenJIJのSQAを用いてこちらを解きます。
以下のようにプログラムを作成しました。
from pyqubo import Binary, Constraint, Placeholder
import openjij as oj
import dimod
import numpy as np
nodes = ["A","B","C","D","E","F","G","H","I","J"]
s_edges = [
("G","I",0.91),("B","D",0.87),("C","E",0.87),("B","E",0.86),("A","H",0.86),
("E","F",0.77),("C","D",0.72),("A","J",0.70),("C","I",0.70),("B","J",0.68),
("C","J",0.64),("B","F",0.61),("E","J",0.60),("C","F",0.60),("B","G",0.59),
("B","H",0.59),("D","J",0.58)
]
# Max-Cut用の重み
w_edges = [(u,v,w) for (u,v,w) in s_edges]
x = {n: Binary(n) for n in nodes}
# 目的項
H_obj = sum(w*(x[u] + x[v] - 2*x[u]*x[v]) for (u,v,w) in w_edges)
# 制約項(バランス制約 5/ 5)
RHO = Placeholder("RHO")
H_bal = RHO * Constraint((sum(x[n] for n in nodes) - 5)**2, label="balance")
H = -H_obj + H_bal
model = H.compile()
Q, offset = model.to_qubo(feed_dict={"RHO": 2.0})
bqm = dimod.BinaryQuadraticModel.from_qubo(Q, offset)
sampler = oj.SQASampler()
seed = 33
num_sweeps = 2000
beta = 5.233147938622838
gamma = 0.4551289384713836
# s を 0→1 に単調増加、★ タプル (s, step_length) にする
s_vals = np.linspace(0.0, 1.0, num_sweeps, dtype=float)
schedule = [(float(si), 1) for si in s_vals]
resp = sampler.sample(bqm, num_reads=2000, num_sweeps=num_sweeps, beta=beta, gamma=gamma, schedule=schedule, seed=seed, trotter=8)
best = resp.first
decoded = model.decode_sample(best.sample, vartype='BINARY', feed_dict={"RHO": 2.0})
print("Energy:", best.energy)
team1 = sorted([n for n in nodes if best.sample[n]==1])
team0 = sorted([n for n in nodes if best.sample[n]==0])
print("Team1:", team1, " (size:", len(team1), ")")
print("Team0:", team0, " (size:", len(team0), ")")
パラメタは以下のような設定をしました。
| seed | 33 |
| num_reads | 2000 |
| num_sweeps | 2000 |
| beta | 5.233147938622838 |
| gamma | 0.4551289384713836 |
結果
以下の結果が得られました。
Energy: -9.250000000000036
Team1: ['D', 'F', 'H', 'I', 'J'] (size: 5 )
Team0: ['A', 'B', 'C', 'E', 'G'] (size: 5 )
今回はグラフのノードとエッジが少ない例であるため、古典的な方法で全探索で「カット重みの最大値9.25」を確認し、SQA の最小エネルギーが−9.25 と一致（同じ割当）しました。
古典的な方法で全探索のプログラム
プログラム
# 5/5（10C5=252通り）を総当たりして、カット重みの最大値を確認するスクリプト
from itertools import combinations
# ノードと重み付きエッジ（本文と同じ）
nodes = ["A","B","C","D","E","F","G","H","I","J"]
w_edges = [
("G","I",0.91),("B","D",0.87),("C","E",0.87),("B","E",0.86),("A","H",0.86),
("E","F",0.77),("C","D",0.72),("A","J",0.70),("C","I",0.70),("B","J",0.68),
("C","J",0.64),("B","F",0.61),("E","J",0.60),("C","F",0.60),("B","G",0.59),
("B","H",0.59),("D","J",0.58)
]
def cut_weight(team1):
"""team1 は 0/1 の '1' 側（サイズ5のタプル/リスト/集合）"""
S = set(team1)
return sum(w for u, v, w in w_edges if (u in S) ^ (v in S))
best = -1.0
best_partitions = []
# ちょうど5人を '1' 側に選ぶ全組合せ（10C5=252）
for team1 in combinations(nodes, 5):
val = cut_weight(team1)
if val > best + 1e-12:
best = val
best_partitions = [team1]
elif abs(val - best) <= 1e-12:
best_partitions.append(team1)
# 結果出力（最大カット重みと、代表解を1つ）
print(f"Best cut weight = {best:.2f} (found in {len(best_partitions)} partition(s))")
rep = set(best_partitions[0])
team1 = sorted(rep)
team0 = sorted(set(nodes) - rep)
print("Team1:", team1, "(size:", len(team1), ")")
print("Team0:", team0, "(size:", len(team0), ")")
出力
Best cut weight = 9.25 (found in 2 partition(s))
Team1: ['A', 'B', 'C', 'E', 'G'] (size: 5 )
Team0: ['D', 'F', 'H', 'I', 'J'] (size: 5 )
所感
今回は実務的なイメージのもと最大カット問題として非相性スコアを定義してアサイン最適化をQUBOで解いてみました。
実際にプログラムを作成した所感として、seedやbeta、gammaなどでかなり得られる最大値がブレる印象を持ちました。記事作成までは理論的な部分の勉強などをしておりプログラムに落とし込めば簡単に解が出てくるものと思っておりましたが、パラメータの調整次第で結果にかなり差異が出る点に気づくことができ、やはり手を動かすのは大事だなと再認識しました。
ここのパラメータ探索は手探りにシードを変えて無理やり最小値を探しましたが、実務で扱うような大きいデータであると総当たりでは最適化も見つけることができず、どの値を目指せばいいのかも認識できなくなるため、パラメータの調整方法や、今回は仮で設定した制約項の係数なども勉強していく必要があると感じました。
今回は擬似量子アニーリングで解いておりますが、量子アニーリングではまた別のパメラメータ調整や埋め込みも出てくるようなのでこちらも合わせて学習が必要になりそうです。
付録
今回のアサイン最適化での例はLucasのレビュー論文の2.2. Graph Partitioningと近いものとなります。
論文ではスピン$s_i \in {±1}$で目的項が$\frac{1 - s_i s_j}{2}$となっております。
スピン↔︎バイナリは$s = 2x - 1$から、上記の目的項に代入することで、$x_i+x_j-2x_i x_j$と同等の式が得られます。また、制約についても同様です。