はじめに
前回記事はこちら
引き続き学習を兼ねて問題を作り、QUBOで解いて行きます。
今回のテーマは、前回最大カット問題を題材としてプロジェクトのアサイン最適化としておりましたが、以下のような変更を加えて取り組みます！
- グループ数：2グループ→3グループ以上
- メンバー：10人→30人
前回は10人を5人づつに分けていましたが、今回は30人を10人づつに分割します。
単純に計算すると、$30C10 \times 20C10$ = 5.5兆くらい（3グループにラベルありの場合）となり、Max-k-Cut（k≥2）は古典的にNP困難です。今回の規模では厳密な全列挙は現実的ではありません。
※問題設定は、chatgptを用いて行っております。
問題設定
メンバー30人を想定として、以下のような設定で10人づつで相性の悪いメンバーの切り離したアサインを目指します。
- メンバー30人
- 非相性スコアの設定：前回記事と同様
この記事での記号と定義
グラフ：無向・重み付き
$$
\boxed{G=(V,E,w)},\quad |V|=30,\ V={A,\ldots,Z,A1,B1,C1,D1},\quad w:E\to\mathbb{R}_{\ge 0}
$$
- 辺集合 (E) は与えられたエッジリスト（u, v, w）のペアのみを採用（＝疎グラフ）。
リストにないペアは辺なし（＝評価対象外）とします。 - 重みは対称 $w_{ij}=w_{ji}$、自己ループは $w_{ii}=0$。
分割：3グループ（10人ずつの等分割）
割当変数（one-hot 表現）：
各頂点 $i\in V$ と各クラスタ $c\in{0,1,2}$ に対して
$$
y_{i,c}\in{0,1},\qquad
\sum_{c=0}^{2} y_{i,c}=1\quad(\forall i\in V)
$$
サイズ（等分割）制約：
$$
\sum_{i\in V} y_{i,c}=10\quad(\forall c\in{0,1,2})
$$
問題を解いてみる
前回は Max-Cut（2分割） で 10 人を 5 人ずつに分け、「異なるグループにまたがる辺重み（カット重み）を最大化」しました。
今回はこれを 3 分割（10/10/10） に拡張します。評価の考え方は同じで、なるべく重い辺が“別グループ”になるように分けるのが目標です。
実装では、各ノードに one-hot 変数 $y_{i, c}\in{{0, 1}}$(ノード$i$がクラスタ$c$に属するなら1)を用いる Potts 型の定式化を採用します。これにより、同じクラスタに入った端点を持つ辺重みの総和を最小化する形で扱えます。直感的には、重い辺ほど同じクラスタに置かれにくくなるため、結果として重い辺は“切られやすく”（別クラスタに振り分けられやすく）なります。
目的項
先ほど少し触れたように以下のようにノード$i$と$j$で繋がるエッジ$E$の重みと、そのノードがどのグループに所属されるかを確認して、同じグループであれば1になり、別のグループであれば0となります。
$$
\min H_{same} = \sum_{(i,j)\in E} w_{ij}\sum_c y_{i, c}y_{j,c}
$$
制約項(ワンホット制約)
こちらの制約では、1つのノードに対して、グループを1つだけ振り分けるための制約です。
$$
H_{onehot} = \lambda \sum_i (\sum_c y_{i,c} - 1)^2
$$
※ $\lambda > 0$は one-hot を強制するペナルティ係数。
制約項(バランス制約)
各クラスタの人数がちょうど10人にするための制約です。
$t_c$は各クラスタに入れたい人数が入っています。
$$
H_{bal} = \mu \sum_c (\sum_i y_{i,c} - t_c)^2
$$
※ 本記事では $t_c=10 $(10/10/10等分) $\mu>0$はバランス制約のペナルティ係数。
最終的なQUBO式
同一クラスタに置かれた辺重みの総和を最小化し，one-hot制約とバランス制約でペナルティで課します。
$$
H = H_{same} + H_{onehot} + H_{bal}
$$
プログラム
今回はPyQUBOとOpenJijのSQAを用いてこちらを解きます。
以下のようにプログラムを作成しました。
import csv
import time
from pyqubo import Binary, Constraint, Placeholder
import openjij as oj
import dimod
import numpy as np
# グラフの読み込み
def load_graph_from_csv(path):
with open(path, newline='', encoding='utf-8') as f:
r = csv.DictReader(f)
edges = [(row['u'], row['v'], float(row['w'])) for row in r]
# ノード集合（u,vに現れるラベルを集める）
nodes = sorted({u for u,_,_ in edges} | {v for _,v,_ in edges})
return nodes, edges
nodes, w_edges = load_graph_from_csv("w_edges_30.csv")
def cut_weight(assign):
total = 0.0
for u, v, w in w_edges:
if assign[u] != assign[v]:
total += w
return total
# quboの構築
def build_bqm_kcut(k=3, lambda_onehot=2.0, mu_balance=0.0, targets=None):
# binary vars
y = { (i, c): Binary(f"{i}_{c}") for i in nodes for c in range(k) }
# 目的項：同じグループに入った端点をもつ辺重みの総和
H_same = 0
for u, v, w in w_edges:
for c in range(k):
H_same += w * y[(u,c)] * y[(v,c)]
# 制約項：ワンホット制約、ノード対して1グループしか適用しない
LAMBDA = Placeholder("LAMBDA")
H_onehot = 0
for i in nodes:
s = sum(y[(i,c)] for c in range(k))
H_onehot += Constraint((s - 1) ** 2, label=f"onehot_{i}")
# 制約項：バランス制約、各グループを10人づつにする
MU = Placeholder("MU")
H_bal = 0
n = len(nodes)
if targets is None:
base = n // k
rem = n % k
targets = [base + (1 if c < rem else 0) for c in range(k)]
for c in range(k):
s = sum(y[(i,c)] for i in nodes)
H_bal += Constraint((s - targets[c]) ** 2, label=f"balance_{c}")
# QUBO式
H = H_same + LAMBDA * H_onehot + MU * H_bal
model = H.compile()
Q, offset = model.to_qubo(feed_dict={"LAMBDA": float(lambda_onehot), "MU": float(mu_balance)})
bqm = dimod.BinaryQuadraticModel.from_qubo(Q, offset)
return model, bqm, targets
def solve_kcut(k=3, lambda_onehot=4.0, mu_balance=3.0, targets=None,
num_reads=1000, num_sweeps=2000, beta=5.233147938622838, gamma=0.4551289384713836,
seed=30, trotter=8):
model, bqm, targets = build_bqm_kcut(k=k, lambda_onehot=lambda_onehot, mu_balance=mu_balance, targets=targets)
sampler = oj.SQASampler()
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=num_reads,
num_sweeps=num_sweeps,
beta=beta,
gamma=gamma,
schedule=schedule,
seed=seed,
trotter=trotter)
best = resp.first
decoded = model.decode_sample(best.sample, vartype='BINARY', feed_dict={"LAMBDA": float(lambda_onehot), "MU": float(mu_balance)})
# 最良のサンプルの中身を確認する
assign = {}
for i in nodes:
scores = []
for c in range(k):
scores.append((best.sample[f"{i}_{c}"], c))
# choose the c with value 1 (or fallback to argmax)
scores.sort(reverse=True)
assign[i] = scores[0][1]
groups = [[] for _ in range(k)]
for i in nodes:
groups[assign[i]].append(i)
groups = [sorted(g) for g in groups]
print(f"Energy: {best.energy:.6f}")
for c, g in enumerate(groups):
print(f"Group {c}: {g} (size: {len(g)})")
cw = cut_weight(assign)
print(f"Cut weight = {cw:.2f}")
if __name__ == "__main__":
start_t = time.time()
solve_kcut(k=3, lambda_onehot=2.0, mu_balance=1.0, targets=[10, 10, 10])
print("Execution time:", time.time() - start_t)
パラメタは以下のような設定をしました。
| $\lambda$ (ワンホット制約) | 4.0 |
| $\mu$ (バランス制約) | 3.0 |
| seed | 14 |
| num_reads | 1000 |
| num_sweeps | 2000 |
| beta | 5.233147938622838 |
| gamma | 0.4551289384713836 |
結果
以下の結果が得られました。
本記事の SQA（約20秒）では Cut weight = 105.45 の解が得られました。
比較として、分枝限定（Branch & Bound）は今回時間制限付きの実行で Cut weight ≈ 102 の最良可行解でした。
同一グラフ・10/10/10制約・ほぼ同時間枠で見る限り、SQA が短時間で良い近似解を返せており、この手法の強みが確認できました。
※ Energy は QUBO の値（制約罰込み）で、Cut weight とは別ものです。評価は制約充足＋Cut weightで比較しています。
Energy: 15.220000
Group 0: ['C', 'D1', 'F', 'J', 'L', 'N', 'O', 'Q', 'R', 'U'] (size: 10)
Group 1: ['A', 'A1', 'B1', 'C1', 'D', 'E', 'H', 'K', 'V', 'Y'] (size: 10)
Group 2: ['B', 'G', 'I', 'M', 'P', 'S', 'T', 'W', 'X', 'Z'] (size: 10)
Cut weight = 105.45
Execution time: 19.272079706192017
本稿の評価値は Cut weight $W_{cut}$（異なるグループにまたがる辺重みの総和）です。
なお Energy は QUBO の値（制約ペナルティ込み）で$W_{cut}$と直接比較しません。
$$
W_{cut} = \sum_{(i,j)\in E} w_{ij}(1-\sum_c y_{i, c}y_{j,c})=\sum_{(i,j)\in E} w_{ij} - H_{same}
$$
上式のように、$W_{cut}$の最大化$\equiv$$H_{same}$の最小化が成り立っています。
おわりに
今回は前回の最大カット問題でのアサイン最適化に引き続き、メンバー数を増やして3分割を行うような問題に取り組みました。メンバーが増えたものの意外とあっさり近似解が得られたのは少しびっくりしました。ペナルティ係数や擬似量子アニーリング特有のパラメータの調整などについてはあまり今回ができておらず、また目的となる厳密解もどれくらいになるのかがわからない場合に、どういったパラメータの設定を行うのか、その探索方法なども次回から取り入れていけると良いなと考えております。
付録
グラフのデータ
今回使ったグラフのデータを以下に貼り付けます。
少し長いので閉じてます。
グラフのcsvデータ
u,v,w
A,D1,0.75
A,F,0.75
A,G,0.65
A,H,0.81
A,Q,0.65
A,S,0.81
A,T,0.55
A,U,0.92
A,V,0.62
A,X,0.88
B,B1,0.85
B,E,0.71
B,J,0.9
B,S,0.8
B,V,0.73
B,Y,0.76
B1,C1,0.6
C,A1,0.71
C,D,0.6
C,G,0.69
C,I,0.58
C,Q,0.57
C,U,0.64
C,W,0.82
C,Y,0.59
D,G,0.89
D,I,0.62
D,J,0.85
D,L,0.69
D,N,0.79
D,P,0.65
D,Q,0.71
D,U,0.61
D,W,0.7
E,C1,0.77
E,D1,0.85
E,I,0.77
E,J,0.61
E,N,0.69
E,O,0.71
E,Q,0.63
E,S,0.57
E,V,0.66
E,X,0.85
E,Y,0.6
F,G,0.9
F,H,0.69
F,I,0.67
F,O,0.71
F,P,0.93
F,T,0.69
F,U,0.87
F,X,0.92
F,Y,0.82
F,Z,0.72
G,A1,0.57
G,B1,0.65
G,C1,0.73
G,D1,0.61
G,H,0.88
G,I,0.69
G,L,0.59
G,O,0.86
G,R,0.75
G,T,0.54
G,U,0.6
G,V,0.57
H,D1,0.83
H,L,0.7
H,O,0.69
H,P,0.66
H,R,0.72
H,S,0.87
H,U,0.77
H,W,0.64
H,Z,0.64
I,A1,0.93
I,B1,0.62
I,C1,0.76
I,L,0.78
I,N,0.63
I,O,0.59
I,R,0.66
I,V,0.83
I,Y,0.86
J,C1,0.81
J,K,0.75
J,L,0.75
J,M,0.86
J,T,0.69
J,W,0.8
J,X,0.68
K,C1,0.87
K,L,0.73
K,M,0.58
K,R,0.68
K,S,0.79
K,T,0.7
K,U,0.81
K,W,0.7
K,Y,0.76
L,C1,0.82
L,O,0.8
L,P,0.63
L,S,0.72
L,T,0.7
L,W,0.6
L,X,0.84
L,Y,0.85
M,A1,0.58
M,C1,0.57
M,D1,0.73
M,N,0.85
M,Q,0.79
M,R,0.82
M,U,0.77
N,B1,0.82
N,S,0.67
N,V,0.58
N,W,0.77
N,X,0.67
N,Y,0.83
O,A1,0.76
O,B1,0.68
O,S,0.81
O,W,0.72
P,A1,0.82
P,D1,0.79
P,U,0.86
P,Y,0.75
Q,A1,0.79
Q,C1,0.92
Q,U,0.81
Q,W,0.52
Q,Z,0.83
R,A1,0.6
R,B1,0.81
R,C1,0.82
R,X,0.78
R,Y,0.71
S,B1,0.69
S,U,0.83
S,W,0.94
S,Y,0.91
T,A1,0.79
T,C1,0.77
T,U,0.77
T,W,0.69
U,A1,0.59
U,B1,0.6
U,C1,0.92
U,V,0.67
U,Y,0.66
V,A1,0.72
V,D1,0.81
V,X,0.62
V,Z,0.78
W,A1,0.55
W,B1,0.59
W,C1,0.85
W,Y,0.57
X,Y,0.86
Y,D1,0.82
Y,Z,0.61
Z,D1,0.61
古典的な手法について
今回の問題の比較対象として組み合わせ最適化の古典的な手法として分岐限定法(branch and bound)を使用しました。
以下chatgptでの簡単な要約です。
分枝限定（Branch & Bound）：
ノードを重み付き次数の降順に並べ、ラベル対称性を破る規則を入れた深さ優先探索。各ステップで「同クラスタに確定した辺重みの総和$W_{same}$」を累積し，未確定辺を全てカット可能と仮定した上界 $UB = W_{total} - W_{same}$により枝刈りする。($W_{toral}: $すべてのエッジ重みの総和)さらに 10/10/10 のサイズ制約は残り人数から到達不能な枝を早期に除去。時間制限付きで走らせると最適性証明はしないが，**最良可行解（incumbent）**を更新し続けるため、短時間でも妥当な厳密系ベースラインとして比較可能。