解説記事「量子アニーリングとADMMのハイブリッド方式による不等式制約への対処」では、不等式制約付きの組合せ最適化問題を解くために、量子アニーリング（QA : Quantum Annealing）と ADMM（Alternating Direction Method of Multipliers）を組み合わせた手法を提案した論文を紹介しました。本記事では、そのアルゴリズムを実装し、元論文の再現実験を行います。
文献情報
- タイトル : Solving Inequality-Constrained Binary Optimization Problems on Quantum Annealer
- 著者 : Kouki Yonaga, Masamichi J. Miyama, Masayuki Ohzeki
- 書誌情報 :
https://doi.org/10.48550/arXiv.2012.06119
手法
はじめに、元論文において提案された手法を簡単に振り返ります（詳細は解説記事に記載しています）。
拡張ラグランジュ法
元論文（および本記事）では、以下のような不等式制約付きの最小化問題を解くことを目的としています。
ここで、
式 (1) の不等式制約は、以下のコスト関数として記述することができます（ここで言うコスト関数とは、対象の組合せ最適化問題を解くために最小化する関数のことです）。
ここで、
しかし、式 (2) の
ここで、
式 (3) に対して拡張ラグランジュ法を適用し、新たなコスト関数
ここで、
メインアルゴリズム
それでは、式 (4) を最小化する ADMM アルゴリズムを以下に示します。
- パラメータの初期化：
- 問題サイズ
- 後述する式 (8) を用いて QUBO 行列を計算する
- QUBO 行列をアニーリングしてサンプル
- サンプル
- 収束の確認：以下の条件のうち一つでも満たされたら計算を終了する（条件内のパラメータ
- 手順 (3)~(10) を収束するまで繰り返す
ここで、
すなわち、サンプル
そして、
すなわち、サンプル
実験
準備：各種ライブラリのインストール
まずは、実装に必要な各種ライブラリをインストールします。
!pip -q -q -q install numpy
!pip -q -q -q install matplotlib
!pip -q -q -q install openjij
!pip -q -q -q install gurobipy
!pip -q -q -q install dwave-system
import copy
import math
import random
from statistics import mean, stdev, variance, median
from tqdm import trange
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import openjij as oj
import gurobipy as gp
from gurobipy import GRB
from dwave.system import DWaveCliqueSampler
ADMM アルゴリズムの実装
各種の補助的な関数の定義
次に、ADMM アルゴリズムの中で補助的に用いる関数を定義していきます。
はじめに、式 (4) の
※ QA では変数
さらに第 3 項について展開し、
def compute_qubo_matrix(f, z, lam, G, D, rho, M, N):
Q = np.array([[f[i, j] for j in range(N)] for i in range(N)]) # まず f(x) の代入
# その他の項の代入
for m in range(M):
for i in range(N):
Q[i, i] += (lam[m] - rho * (D[m] + z[m] - 0.5 * G[m, i])) * G[
m, i
] # 対角項
for j in range(i + 1, N):
Q[i, j] += rho * G[m, i] * G[m, j] # 非対角項
return Q
次に、上の compute_qubo_matrix 関数から得られた QUBO 行列を QA によって最適化する関数 annealing_qubo_matrix を定義します。
def annealing_qubo_matrix(Q, N_nu, sampler):
# N_nu はサンプル数
sampleset = sampler.sample_qubo(Q, num_reads=N_nu)
return sampleset.record.sample
続いて、サンプル
def compute_E_aug(x, f, z, lam, G, D, rho, M):
E_aug = np.dot(np.dot(x, f), x) # f(x) の計算
for m in range(M):
Gx = np.dot(G[m, :], x)
E_aug += lam[m] * Gx + rho * 0.5 * (Gx**2) - rho * (D[m] + z[m]) * Gx
return E_aug
式 (5) に基づいて
def compute_x_cost(x_nu, f, z, lam, G, D, rho, M):
# まず最初のインデックスのサンプルを入れておく
E_temp = compute_E_aug(x_nu[0], f, z, lam, G, D, rho, M)
x_cost = x_nu[0]
# 二番目以降のサンプルについて E_aug を最小にするもの (x_cost*) を探索
for nu in range(1, len(x_nu)):
E_aug = compute_E_aug(x_nu[nu], f, z, lam, G, D, rho, M)
if E_aug <= E_temp:
E_temp = E_aug
x_cost = x_nu[nu]
return x_cost
式 (6) に基づいて
def compute_x_feas(x_nu, f, G, D, M):
f_list = []
x_feas_list = []
for x in x_nu: # すべてのサンプルについて
# すべての m について不等式条件を満たしているか確認
if not all(np.dot(G[m, :], x) <= D[m] for m in range(M)):
continue
f_temp = np.dot(x, np.dot(f, x)) # f(x) の計算
f_list.append(f_temp)
x_feas_list.append(x)
if len(x_feas_list) == 0: # 不等式条件を満たすリストが無かった場合
return []
else:
return x_feas_list[
f_list.index(min(f_list))
] # x_feas_list の中で、f(x) を最小にするもの (x_feas*) を返す
サンプル
def compute_E_ineq(f, gamma, G, x, D, M):
E_ineq = np.dot(x, np.dot(f, x)) # f(x) の計算
for m in range(M):
f_step = int(np.dot(G[m, :], x) > D[m]) # ステップ関数 Theta(G*x - D)
E_ineq += gamma * f_step
return E_ineq
アルゴリズムの手順 8. で 3 つの収束条件をチェックする関数 check_convergence を定義します。また、check_convergence 関数内で使用する compute_criteria3 (収束条件 3 をチェックする関数) も併せて定義します。
def compute_criteria3(x_feas, G, D, z, M, epsilon):
sum_c3 = np.sum([(np.dot(G[m, :], x_feas) - D[m] - z[m]) ** 2 for m in range(M)])
return (
math.sqrt(sum_c3) < epsilon
) # 収束条件 3 を満たす場合には True、そうでなければ False を返す
def check_convergence(t, t_max, E_ineq, x_feas, G, D, z, t_conv, epsilon, M, log=True):
if t > t_max: # 収束条件 1 を満たす場合
if log:
print("condition : criteria 1 is satisfied.")
return True
if (
len(E_ineq) >= t_conv and len(set(E_ineq[-t_conv:])) <= 1
): # 収束条件 2 を満たす場合
if log:
print("condition : criteria 2 is satisfied.")
return True
if len(x_feas) != 0 and compute_criteria3(
x_feas, G, D, z, M, epsilon
): # 収束条件 3 を満たす場合
if log:
print("conditon : criteria 3 is satisfied.")
return True
return False
ADMM アルゴリズム本体の定義
これまで作成してきた関数を用いて ADMM アルゴリズムの本体を定義していきます。アルゴリズムの手順 2. に対応する embedding では DWaveCliqueSampler() を用いることで全結合グラフを埋め込むことができます。これにより、全結合グラフ （clique） 用に予め用意された埋め込みが適用されることになります。また、今回は QA として D-Wave Advantage system6.3 を使用します。
def qa_admm(f, G, D, rho, M, N_nu, N, t_max, t_conv, epsilon, gamma):
E_ineq = []
# 手順 1 : パラメータの初期化
z = np.zeros(M)
lam = np.zeros(M)
t = 1
# 手順 2 : 全結合グラフ (サイズ N) の embedding を行う ( = sampler の指定)
dw_sampler = DWaveCliqueSampler(token="YOUR TOKEN", solver="Advantage_system6.3")
# SA の場合 (埋め込みは無し)
# sa_sampler = oj.SASampler()
while True: # 収束条件を満たすまで繰り返す
# print("-" * 50)
# print("number of iterations = ", t)
# 手順 3 : QUBO 行列を計算する
Q = compute_qubo_matrix(f, z, lam, G, D, rho, M, N)
"""
plt.imshow(Q)
plt.colorbar()
"""
# 手順 4 : QUBO 行列をアニーリングしてサンプルを取得
x_nu = annealing_qubo_matrix(Q, N_nu, dw_sampler)
# 手順 5 : x_cost* と x_feas* を計算する
x_cost = compute_x_cost(x_nu, f, z, lam, G, D, rho, M)
# print("x_cost = ", x_cost)
x_feas = compute_x_feas(x_nu, f, G, D, M)
# print("x_feas = ", x_feas)
for m in range(M):
Gx_cost = np.dot(G[m, :], x_cost)
# 手順 6 : z をアップデートする
z[m] = min(0, Gx_cost - D[m])
# print("z[0] = ", z[m])
# 手順 7 : lambda をアップデートする
lam[m] += rho * (Gx_cost - D[m] - z[m])
# print("dot_Gmx_cost = ", Gx_cost)
# print("dot_Gmx_cost-D[m]-z[m] = ", Gx_cost-D[m]-z[m])
# print("lambda[0] = ", lam[m])
# 手順 8 : 収束しているかどうかをチェックする
if len(x_feas) != 0:
E_ineq.append(compute_E_ineq(f, gamma, G, x_feas, D, M))
# print("E_ineq = ", E_ineq)
if check_convergence(
t, t_max, E_ineq, x_feas, G, D, z, t_conv, epsilon, M, log=False
):
break
# 手順 9 : 反復回数 t を更新する
t += 1
return x_feas, E_ineq
解きたい最適化問題（QKP）の定義
実装した ADMM アルゴリズムの精度を評価するために、元論文と同様に二次ナップサック問題（QKP : Quadratic Knapsack Problem）を解いていきます。QKP は以下のように定義されます。
式 (8) において、
実験では、
アルゴリズムの典型的な性能を評価するために、このランダムな QKP を 10 個生成します（すなわちインスタンスの数
なお、
以下に、ランダムな QKP を生成する関数 make_QKP を定義します。
def make_QKP(N, M, Delta):
f = np.zeros((N, N))
G = np.zeros((M, N))
D = np.zeros((M))
for i in range(N):
for j in range(i, N):
if random.uniform(0, 1) <= Delta:
f[j, i] = -random.randint(1, 100) # 最大化問題なので、符号を反転させる
for i in range(N):
G[0, i] = random.randint(1, 50)
D[0] = random.randint(50, np.sum(G[0, :]))
return f, G, D
ADMM アルゴリズムの性能評価
rho = 0.1
t_max = 30
t_conv = 10
epsilon = 10 ** (-3)
gamma = 100
N_nu = 2000 # サンプル数
N_inst = 10 # インスタンスの数
N_list = [8, 16, 32, 64]
len_N_list = len(N_list)
M = 1 # 不等式制約の数（式(8)より１個）
問題サイズ
Delta02 = 0.2
f_inst_D02 = []
G_inst_D02 = []
D_inst_D02 = []
for N in N_list:
for i in range(N_inst):
f, G, D = make_QKP(N, M, Delta02)
f_inst_D02.append(f)
G_inst_D02.append(G)
D_inst_D02.append(D)
ADMM アルゴリズムの実行
生成した QKP（計 40 個）に対して ADMM アルゴリズムを実行し、解を求めていきます。
x_feas_inst_D02 = []
E_ineq_inst_D02 = []
for i in trange(N_inst * len_N_list):
x_feas, E_ineq = qa_admm(
f_inst_D02[i],
G_inst_D02[i],
D_inst_D02[i],
rho,
M,
N_nu,
N_list[i // N_inst],
t_max,
t_conv,
epsilon,
gamma,
)
x_feas_inst_D02.append(x_feas)
E_ineq_inst_D02.append(E_ineq)
Gurobi から厳密解を取得
今度はそれぞれの QKP に対して Gurobi を実行し、厳密解を求めていきます。まずは Gurobi で厳密解を求める関数 grb_exact_qkp を定義します。
def grb_exact_qkp(P, w, c):
N = len(w)
# 問題を設定
QKP = gp.Model(name="QKP")
# 変数を設定（変数単体にかかる制約を含む）
x = {}
for i in range(N):
x[i] = QKP.addVar(vtype="B", name="x[%x]" % i)
# 目的関数を設定
QKP.setObjective(
sum(sum(x[i] * P[i, j] * x[j] for i in range(N)) for j in range(N)),
GRB.MINIMIZE,
)
# 制約を設定
QKP.addConstr(sum(w[i] * x[i] for i in range(N)) <= c, name="constraint")
# ログを切る
QKP.params.LogToConsole = 0
# 解を求める
QKP.optimize()
if QKP.Status != GRB.OPTIMAL:
raise Exception("最適解が得られませんでした")
return np.array([x[i].x for i in range(N)]), QKP.ObjVal # 厳密解とそのコストを返す
この関数を用いて、各 QKP の厳密解を求めます。
x_opt_inst_D02 = []
E_opt_inst_D02 = []
for k in trange(N_inst * len_N_list):
x_exact, val_opt = grb_exact_qkp(f_inst_D02[k], G_inst_D02[k][0], D_inst_D02[k][0])
x_opt_inst_D02.append(x_exact)
E_opt_inst_D02.append(val_opt)
結果の比較
ADMM アルゴリズムで得られた解と Gurobi の厳密解とを比較します。厳密解と比較した ADMM アルゴリズムの精度をここでは MAPE で評価します。MAPE は以下のように定義されます。
ここで、
式(9)に基づき、問題サイズごとに MAPE のインスタンス平均および標準誤差を求めます。
mean_MAPE_D02 = np.zeros(len_N_list)
stder_MAPE_D02 = np.zeros(len_N_list)
for i in trange(len_N_list):
# 最適値の符号を反転させることで、元の最大化問題の最適値に戻す
mean_MAPE_D02[i] = mean(
[
abs(E_ineq_inst_D02[N_inst * i + j][-1] - E_opt_inst_D02[N_inst * i + j])
/ -E_opt_inst_D02[N_inst * i + j]
for j in range(N_inst)
if E_ineq_inst_D02[N_inst * i + j][-1] is not None
]
)
stder_MAPE_D02[i] = stdev(
[
abs(E_ineq_inst_D02[N_inst * i + j][-1] - E_opt_inst_D02[N_inst * i + j])
/ -E_opt_inst_D02[N_inst * i + j]
for j in range(N_inst)
if E_ineq_inst_D02[N_inst * i + j][-1] is not None
]
) / np.sqrt(N_inst)
実験結果として、MAPE の問題サイズ
log2N_list = [3, 4, 5, 6]
plt.rcParams["font.size"] = 13 # フォントの大きさ
plt.errorbar(
log2N_list,
mean_MAPE_D02,
yerr=stder_MAPE_D02,
marker="o",
capthick=1,
capsize=10,
lw=1,
)
plt.gca().xaxis.set_major_locator(
ticker.MaxNLocator(integer=True)
) # 横軸を整数値で表示
plt.ylabel("MAPE")
plt.xlabel("$log_{2}N$")
plt.title("$\Delta=0.2$")
plt.show()
Delta06 = 0.6
f_inst_D06 = []
G_inst_D06 = []
D_inst_D06 = []
for N in N_list:
for i in range(N_inst):
f, G, D = make_QKP(N, M, Delta06)
f_inst_D06.append(f)
G_inst_D06.append(G)
D_inst_D06.append(D)
ADMM アルゴリズムの実行
生成したそれぞれの QKP に対して ADMM アルゴリズムを実行し、解を求めていきます。
x_feas_inst_D06 = []
E_ineq_inst_D06 = []
for i in trange(N_inst * len_N_list):
x_feas, E_ineq = qa_admm(
f_inst_D06[i],
G_inst_D06[i],
D_inst_D06[i],
rho,
M,
N_nu,
N_list[i // N_inst],
t_max,
t_conv,
epsilon,
gamma,
)
x_feas_inst_D06.append(x_feas)
E_ineq_inst_D06.append(E_ineq)
Gurobi から厳密解を取得
続いてそれぞれの QKP に対して Gurobi を実行し、厳密解を求めていきます。
x_opt_inst_D06 = []
E_opt_inst_D06 = []
for k in trange(N_inst * len_N_list):
x_exact, val_opt = grb_exact_qkp(f_inst_D06[k], G_inst_D06[k][0], D_inst_D06[k][0])
x_opt_inst_D06.append(x_exact)
E_opt_inst_D06.append(val_opt)
結果の比較
式 (9) に基づき、問題サイズごとに MAPE のインスタンス平均および標準偏差を求めます。
mean_MAPE_D06 = np.zeros(len_N_list)
stder_MAPE_D06 = np.zeros(len_N_list)
for i in range(len_N_list):
# 最適値の符号を反転させることで、元の最大化問題の最適値に戻す
mean_MAPE_D06[i] = mean(
[
abs(E_ineq_inst_D06[N_inst * i + j][-1] - E_opt_inst_D06[N_inst * i + j])
/ -E_opt_inst_D06[N_inst * i + j]
for j in range(N_inst)
if E_ineq_inst_D06[N_inst * i + j][-1] is not None
]
)
stder_MAPE_D06[i] = stdev(
[
abs(E_ineq_inst_D06[N_inst * i + j][-1] - E_opt_inst_D06[N_inst * i + j])
/ -E_opt_inst_D06[N_inst * i + j]
for j in range(N_inst)
if E_ineq_inst_D06[N_inst * i + j][-1] is not None
]
) / np.sqrt(N_inst)
実験結果として、MAPE の問題サイズ
plt.errorbar(
log2N_list,
mean_MAPE_D06,
yerr=stder_MAPE_D06,
marker="o",
capthick=1,
capsize=10,
lw=1,
)
plt.gca().xaxis.set_major_locator(
ticker.MaxNLocator(integer=True)
) # 横軸を整数値で表示
plt.ylabel("MAPE")
plt.xlabel("$log_{2}N$")
plt.title("$\Delta=0.6$")
plt.show()
問題サイズ
Delta10 = 1.0
f_inst_D10 = []
G_inst_D10 = []
D_inst_D10 = []
for N in N_list:
for i in range(N_inst):
f, G, D = make_QKP(N, M, Delta10)
f_inst_D10.append(f)
G_inst_D10.append(G)
D_inst_D10.append(D)
ADMM アルゴリズムの実行
生成したそれぞれの QKP に対して ADMM アルゴリズムを実行し、解を求めていきます。
x_feas_inst_D10 = []
E_ineq_inst_D10 = []
for i in trange(N_inst * len_N_list):
x_feas, E_ineq = qa_admm(
f_inst_D10[i],
G_inst_D10[i],
D_inst_D10[i],
rho,
M,
N_nu,
N_list[i // N_inst],
t_max,
t_conv,
epsilon,
gamma,
)
x_feas_inst_D10.append(x_feas)
E_ineq_inst_D10.append(E_ineq)
Gurobi から厳密解を取得
続いてそれぞれの QKP に対して Gurobi を実行し、厳密解を求めていきます。この処理には数時間かかります。
x_opt_inst_D10 = []
E_opt_inst_D10 = []
for k in trange(N_inst * len_N_list):
x_exact, val_opt = grb_exact_qkp(f_inst_D10[k], G_inst_D10[k][0], D_inst_D10[k][0])
x_opt_inst_D10.append(x_exact)
E_opt_inst_D10.append(val_opt)
結果の比較
式 (9) に基づき、問題サイズごとにMAPE のインスタンス平均および標準偏差を求めます。
mean_MAPE_D10 = np.zeros(len_N_list)
stder_MAPE_D10 = np.zeros(len_N_list)
for i in range(len_N_list):
# 最適値の符号を反転させることで、元の最大化問題の最適値に戻す
mean_MAPE_D10[i] = mean(
[
abs(E_ineq_inst_D10[N_inst * i + j][-1] - E_opt_inst_D10[N_inst * i + j])
/ -E_opt_inst_D10[N_inst * i + j]
for j in range(N_inst)
if E_ineq_inst_D10[N_inst * i + j][-1] is not None
]
)
stder_MAPE_D10[i] = stdev(
[
abs(E_ineq_inst_D10[N_inst * i + j][-1] - E_opt_inst_D10[N_inst * i + j])
/ -E_opt_inst_D10[N_inst * i + j]
for j in range(N_inst)
if E_ineq_inst_D10[N_inst * i + j][-1] is not None
]
) / np.sqrt(N_inst)
実験結果として、MAPE の問題サイズ
plt.errorbar(
log2N_list,
mean_MAPE_D10,
yerr=stder_MAPE_D10,
marker="o",
capthick=1,
capsize=10,
lw=1,
)
plt.gca().xaxis.set_major_locator(
ticker.MaxNLocator(integer=True)
) # 横軸を整数値で表示
plt.ylabel("MAPE")
plt.xlabel("$log_{2}N$")
plt.title("$\Delta=1.0$")
plt.show()
実験結果のまとめ
plt.errorbar(
log2N_list,
mean_MAPE_D02,
yerr=stder_MAPE_D02,
linestyle="solid",
marker="o",
capthick=1,
capsize=10,
lw=1,
label="$\Delta=0.2$",
)
plt.errorbar(
log2N_list,
mean_MAPE_D06,
yerr=stder_MAPE_D06,
linestyle="dashed",
marker="x",
capthick=1,
capsize=10,
lw=1,
label="$\Delta=0.6$",
)
plt.errorbar(
log2N_list,
mean_MAPE_D10,
yerr=stder_MAPE_D10,
linestyle="dotted",
marker="v",
capthick=1,
capsize=10,
lw=1,
label="$\Delta=1.0$",
)
plt.gca().xaxis.set_major_locator(
ticker.MaxNLocator(integer=True)
) # 横軸を整数値で表示
plt.xlabel("$log_{2}N$")
plt.ylabel("MAPE")
plt.legend()
元論文では
まとめ
本記事では、元論文の ADMM アルゴリズムを実装しました。そして、再現実験として
元論文とは異なり、インスタンスによって MAPE が大きく変化する結果となりました。また、
あとがき
元論文の結果が再現されたとは言い難い結果とはなりましたが、拡張ラグランジュ法や ADMM の勉強にもなり、良い経験になりました。また、今回は QKP という不等式制約が 1 つだけの問題を扱いましたが、より多くの不等式制約が含まれる問題も今回の手法で扱えるかどうかを確かめてみたいと思いました。