この記事は Jij Inc. Advent Calendar 2024 の 10 日目の記事です．
Jij Inc. の坂井です．
Jij が開発している量子最適化計算向けの SDK，Qamomile をザックリ使ってみます．問題を記述することとそれを解くための量子回路を組むことの間にはギャップがあるものですが，それを埋めてくれるに違いありません．
主な特徴
- 問題は JijModeling でなるべく直感的に記述できる．
- 組合せ最適化問題を解く手法である QAOA，古典ビットの量子ビットへの効率的な埋め込み手法である QRAO，といった有名なアルゴリズムをサポートしている．
- 代表的な回路への自動変換だけでなく，独自の中間表現を用いた柔軟なプログラミングもできる．
- Qiskit 回路や Quri Parts 回路への変換もできる[1]．
インストール
求解には Qiskit を経由することにして，
$ pip install "qamomile[qiskit]"
とインストールしましょう．
問題の記述には JijModeling，問題のコンパイル[2]には JijModeling-Transpiler，を使用するので，それらもインストールしてください．
$ pip install jijmodeling
$ pip install jijmodeling-transpiler
数値的なあれこれのために NumPy と SciPy も使うので，これらもお使いの環境にインストールされていることを仮定します．
Multi-car paint shop problem
この記事では Multi-car paint shop (MCPS) 問題と呼ばれる問題を解いてみましょう[3]．ザックリ言えばこの問題は，複数の車種の塗装時における色の切り替え回数を最小化する問題です．人力で塗る場合は違う色のペンキの缶をいちいち持ち替えるのが面倒くさいということですね．この問題は NP 困難というクラスに属することが知られており，常に簡単に解けるわけでは全然ないと言えます．
問題
塗装工程では，発注により決まっている表面の塗装 (ベースコート) と，その下地の塗装 (フィラー) が行なわれます．フィラーの色はベースコートの色によって決定され，
- ベースコートが明るい色ならばフィラーは白色
- ベースコートが暗い色ならばフィラーは黒色
という決まりがあります．MCPS はこのフィラー塗装の最適化問題です．待機列に並ぶ車を，連続するものはできるだけ同じフィラー色で塗装しながら，各車種に対して決まった台数分正しいフィラー色を適用する必要があります．
例として，2 つの車種
これからこの問題を数理モデルとして記述し，QAOA により解いてみましょう．
数理モデル
ただし，
先程の例のように，MCPS 問題では各車種について黒で塗るべき台数が決まっています．車種
この問題の目的は色の切り替え回数を最小することでした．つまり待機列中の連続する車両はできるだけ同じ色で塗装したいです．連続する車両の間で塗る色が変わったときにコストがかかると思うことにし，MCPS はそれを最小化する問題だと思うことにしましょう:
ここで，
- 連続する車両を同じ色に塗る (
s_i = s_{i+1} - s_i s_{i+1} - 連続する車両を違う色に塗る (
s_i \neq s_{i+1} - s_i s_{i+1}
と読むことができます．
同じような意味を持つ
モデルのコーディング
import jijmodeling as jm
N = jm.Placeholder("N") # 車両の数
M = jm.Placeholder("M") # 車種の数
C = jm.Placeholder("C", ndim=2) # C[m] が車種が m である車両のリストを表す 2 次元リスト
K = jm.Placeholder("K", ndim=1) # K[m] が車種が m である車両のうち黒色に塗りたいものの数を表すリスト
q = jm.BinaryVar("q", shape=(N,))
i = jm.Element("i", belong_to=(0, N-1))
m = jm.Element("m", belong_to=(0, M))
j = jm.Element("j", belong_to=C[m])
problem = jm.Problem("MCPS")
problem += jm.sum([i], -(2*x[i] - 1)*(2*x[i+1] - 1)) # 目的関数
problem += jm.Constraint("K-hot", jm.sum([l], q[l]) == K[m], forall=m) # 制約条件
JijModeling では，式の構成要素である定数，変数，添字を予め定義し，それを用いて問題を記述します．Placeholder や Element という語法には慣れが必要かもしれませんが，問題の記述は数式に対応していていくらか直感的といえるのではないでしょうか．
面白いのは配列の添字として使う予定の文字も範囲と共に予め宣言しておく点です．これにより jm.sum([l], q[l]) と，数式に対応した形で記述できるのです．
もう一つ重要なのは，この時点ではインスタンス[7]に依存した具体的な数値が現れていない点です．あくまで文字を用いて式を記述し，その具体的な中身は後で渡す，というのが JijModeling を用いた数理モデリングのやり方です．
添字の範囲が整数の区間としてだけではなく belong_to=C[m] とも書け，
問題が記述できたので，具体的なインスタンスを与えてみましょう．最後に行なっているインスタンスのコンパイルには，ここではそれほど注意する必要はありません．
import random as rand
import jijmodeling_transpiler.core as jtc
carmodel = {0: "V", 1: "T", 2: "S", 3: "P"} # バン，タクシー，SUV，パトカーということにしましょう．
numblack = [1, 1, 1, 1] # 全ての車種に対して，1 車両ずつ黒いフィラーを塗りたいとします．
cars = [0, 0, 1, 1, 2, 2, 3, 3] # 各車種 2 台ずつあるとし，
rand.shuffle(cars) # 待機列にはシャッフルされた順番で入ってくるとします．
C = [[i for i in range(len(cars)) if cars[i] == m] for m in range(len(carmodel))]
# C[m] が車種が m である車両のリストを表す 2 次元リスト
instance_data = {"N": len(cars), "M": len(carmodel), "C": C, "K": numblack}
compiled_model = jmt.compile_model(problem, instance_data)
Qamomile を用いた QAOA 回路とハミルトニアンの生成
Qamomile にはコンパイルされたインスタンスから QAOA を実行するための回路とハミルトニアンを生成する機能が備わっています．もちろんそうした変換は常に自明ではなく，Qamomile に細かなパラメータを渡すこともできます．
QAOA の一般的なワークフローについても簡単に説明しながら進みます．
from qamomile.core.converters.qaoa import QAOAConverter
from qamomile.core.circuit.drawer import plot_quantum_circuit
converter = QAOAConverter(compiled_model)
converter.ising_encode(multipliers={"K-hot": 3})
circuit = converter.get_qaoa_ansatz(p=1)
hamiltonian = converter.get_cost_hamiltonian()
plot_quantum_circuit(qaoa_circuit)
Qamomile により生成された量子回路．
converter.ising_encode() という関数は
種々の制約条件は単純には罰則項として目的関数に埋め込まれます．ここで指定する multipliers が大きいほど対応する制約条件が厳しく守られると思ってください．
ザックリ言うと，QAOA とは変分パラメータを含んだ量子回路を用意し，そのパラメータを調節することで終状態が良い状態，つまり目的関数の期待値が小さい状態，になるようにするアルゴリズムです．converter.get_qaoa_ansatz() に渡している p は変分パラメータの数を決めるものだと思ってください．これが大きいほど回路の表現力は大きくなると言えますが，一般にコストも増大し (実機では) ノイズ源も多くなるという側面もあるため，p はとにかく大きいほど良いというわけではないことに注意が必要です．
また，ここで言うハミルトニアンは，その期待値が今解きたい問題の目的関数の期待値を表すようなものが構築されます．
Qiskit が提供するシミュレータを使用するために，circuit と hamiltonian を Qiskit のクラスに変換しましょう．
import qamomile.qiskit as qm_qk
transpiler = qm_qk.QiskitTranspiler()
qk_circuit = transpiler.transpile_circuit(circuit)
qk_hamiltonian = transpiler.transpile_hamiltonian(hamiltonian)
QAOA の実行と結果
全ての準備が整いました．いよいよ QAOA を動かしてみましょう．
QAOA はハミルトニアンの期待値の評価を量子デバイス (今はシミュレータ) で行ない，上述の変分パラメータの調節を古典デバイスで行なうアルゴリズムです．ここでは古典側の最適化は Scipy の COBYLA で行なうことにしましょう．
import qiskit.primitives as qk_pr
import numpy as np
from scipy.optimize import minimize
history = []
estimator = qk_pr.StatevectorEstimator()
def costfunc(params):
job = estimator.run([(qk_circuit, qk_hamiltonian, params)])
result = job.result()[0]
cost = result.data['evs']
history.append(cost)
return cost
result = minimize(
costfunc,
[np.random.rand(2)*np.pi] * 2,
method="COBYLA",
options={"maxiter": 100},
)
print(result)
costfunc() は QAOA の変分パラメータの具体値を受け取って，それを埋め込んだ回路で評価したハミルトニアンの期待値を返す関数になっています．result の中身は次のようになり，result.x が調節された QAOA の変分パラメータです．
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.4307546700453164
x: [-3.776e-01 1.015e+00]
nfev: 68
maxcv: 0.0
最適化の過程，つまりコスト関数がどのように減少していくかは history に格納されています．
import matplotlib.pyplot as plt
plt.title("Change of Cost", fontsize=15)
plt.xlabel("Iteration", fontsize=15)
plt.ylabel("Cost", fontsize=15)
plt.xscale("log")
plt.xlim(1, result.nfev)
plt.plot(cost_history, label="Cost", color="#2696EB")
plt.show()
最適化の過程．
最適な変分パラメータを使ったときに終状態がどのようなものになるか見てみましょう．
sampler = qk_pr.StatevectorSampler()
qk_circuit.measure_all()
job = sampler.run([(qk_circuit, result.x)], shots=1000)
jobresult = job.result()[0]
counts = jobresult.data["meas"]
ごく簡単に可視化してみます．
color = {
1: "⚫", # 黒塗り
0: " " # 白塗り
}
sampleset = converter.decode(qk_transpiler, jobresult.data["meas"])
maxenergy = -1e9
bestvars = None
for sample in sampleset.feasibles():
if maxenergy < sample.eval.objective:
maxenergy = sample.eval.objective
bestvars = sample.var_values
sol = [0] * 8
for idx in bestvars["q"].values:
sol[idx[0]] = 1
print("cars intake: ", [carmodel[car] for car in cars])
print("color: ", [color[b] for b in sol])
converter.decode() を用いてサンプル結果を取得していることに注意してください．手動でやると間違えやすいので，Qamomile の便利なところの一つです．結果は以下のようになります．
cars intake: ['S', 'P', 'V', 'P', 'T', 'S', 'V', 'T']
color: [' ', ' ', ' ', '⚫', ' ', '⚫', '⚫', '⚫']
パッと見には，少なくとも各車種が 1 台ずつ黒く塗られており，制約条件を満たしていることが分かります．この解では色を切り替える回数が 3 回になっています．このインスタンスはサイズが小さいので，詳細は省略しますが総当たり的な方法で最適解を探してみると，ちゃんと 3 回という結果が得られます．
得たサンプルがどのような度数分布を描くかも見ておきましょう．
from collections import defaultdict
freq = defaultdict(int)
for sample in sampleset:
freq[round(sample.eval.objective, ndigits = 3)] += sample.num_occurrences
plt.bar(frequencies.keys(), frequencies.values(), width=0.25)
plt.xticks(sorted(np.array(list(freq.keys()))))
plt.xlabel('Objective', fontsize=15)
plt.ylabel('Frequency', fontsize=15)
plt.show()
得られたエネルギー状態の度数分布
このようにして見ると最低エネルギー状態は全体のごく僅かな割合でしか観測されておらず，少し心許ない気もします．理想的には上述の converter.get_qaoa_ansatz() 関数に渡していた p というパラメータを増やせば分布の質が改善していくと期待されます．
まとめ
Qamomile を使って，QAOA による MCPS 問題の求解を試してみました．古典ビットを量子ビットに対応付ける (エンコード/デコード) だとか，Pauli 演算子でハミルトニアンを記述するといった部分がかなり単純化されます．JijModeling による数理モデルの快適な記述と併せて，Qamomile が皆さんの量子最適化体験を豊かにすることを期待します．
募集
Jijでは各ポジションを絶賛採用中です！
ご興味がございましたらカジュアル面談からでもぜひ！ ご連絡お待ちしております．
数理最適化エンジニア
量子コンピューティングソフトウェアエンジニア
Rust エンジニア (アルゴリズムエンジニア)
研究チーム PMO
オープンポジション等その他の採用情報はこちら
Discussion