はじめに
通常交通最適化問題では制約条件を設定することで、最適化問題を量子コンピュータで解きますが、概して計算速度や精度に影響します。ここではblueqatでQAOAというアルゴリズムを改造することで、多くの最適化問題において制約条件を定式化することなく最適化する方法を紹介します。
ケーススタディとして交通最適を取り上げます。
量子コンピュータでの交通最適化モデル化
参考資料がありますので、興味ある人は、
「交通最適化問題を量子アニーリングで解く」
今回はスタートAからゴールBへ向かう地図を考え、各道にs0からs11までの通し番号をつけます。
車を二台用意し、それぞれ取りうるルートを設定します。
自動車１
ルート１-1(q0)：s0,s3,s6,s9
ルート1-２(q1)：s0,s3,s8,s11
自動車２
ルート2-1(q2):s0,s3,s8,s11
ルート2-2(q3):s2,s7,s10,s11
とします。カッコ内は今回利用する量子ビットで、そのルートを選ぶ際にはqは1、選ばない時には0となるように計算をします。例えば、自動車１がルート1-1(q0)をとって、自動車２がルート2-2(q3)を取る時には、(q0,q1,q2,q3) = (1,0,0,1)のように配列がなります。
通常の解き方、制約条件項あり
通常は、各道の混雑を足し合わせます。混雑はs0からs11までの道の混雑を量子ビットを使って表現します。混雑度のコスト関数は、
また、制約条件は自動車１、自動車２それぞれルートを１つだけ選ぶという条件で、
通常は、
というふうに、コスト関数に制約条件に調整変数をかけたものを足し合わせます。今回は二番目のConstが不要になるテクニックを導入します。
初期状態を準備
量子アニーリングでは探索はほぼランダムで行われます。初期状態は|+>からスタートし、途中過程の探索は
しかしここでは、ルート選択は１つだけを選ぶわけですから、
ここで、量子回路でまず量子もつれを作ります。量子もつれは多くのパターンから限られた組み合わせだけが出るように調整できるのでした。簡単な例として、
from blueqat import Circuit
Circuit().h[0].cx[0,1].m[:].run(shots=100)
#Counter({'00': 51, '11': 49})
Counter({'00': 52, '11': 48})
こうすると、本来４種類出るところから00と11のみ答えが出ます。
ちょっと形を変えることで、01と10のもつれに変更できます。
from blueqat import Circuit
Circuit().h[0].cx[0,1].x[0].m[:].run(shots=100)
#Counter({'01': 41, '10': 59})
Counter({'01': 53, '10': 47})
これを使います。つまり、自動車１と自動車２にそれぞれもつれを活用して、01と10に制限をかけるように初期状態を設定します。
今回４量子ビットの回路では、
from blueqat import Circuit
Circuit().h[0].cx[0,1].x[0].h[2].cx[2,3].x[2].m[:].run(shots=100)
#Counter({'0101': 27, '0110': 16, '1001': 29, '1010': 28})
Counter({'1010': 20, '0110': 25, '1001': 20, '0101': 35})
この４通りに絞ります。
入れ替え作業
これにより、初期状態で準備した01と10をひたすら入れ替えながら、量子断熱計算で、ハミルトニアンの最小値を求めるQAOAと組み合わせることができます。
具体的コード
具体的な実行コードはこちらです。元々のQAOAのAnsatzを改造しています。
import numpy as np
from blueqat import Circuit
from blueqat.pauli import X, Y, Z, I
from blueqat.pauli import qubo_bit as q
from blueqat.vqe import AnsatzBase, Vqe
def an(index):
return 0.5 * X[index] + 0.5j * Y[index]
def cr(index):
return 0.5 * X[index] - 0.5j * Y[index]
op1 = (cr(1) * an(0) + cr(0) * an(1)).to_expr().simplify()
op2 = (cr(3) * an(2) + cr(2) * an(3)).to_expr().simplify()
class QaoaQaoaAnsatz(AnsatzBase):
def __init__(self, hamiltonian, step=1, init_circuit=None):
super().__init__(hamiltonian, step * 2)
self.hamiltonian = hamiltonian.to_expr().simplify()
if not self.check_hamiltonian():
raise ValueError("Hamiltonian terms are not commutable")
self.step = step
self.n_qubits = self.hamiltonian.max_n() + 1
if init_circuit:
self.init_circuit = init_circuit
if init_circuit.n_qubits > self.n_qubits:
self.n_qubits = init_circuit.n_qubits
else:
self.init_circuit = Circuit(self.n_qubits).h[:]
self.init_circuit.make_cache()
self.time_evolutions = [term.get_time_evolution() for term in self.hamiltonian]
self.time_evolutions1 = [term.get_time_evolution() for term in op1]
self.time_evolutions2 = [term.get_time_evolution() for term in op2]
def check_hamiltonian(self):
"""Check hamiltonian is commutable. This condition is required for QaoaAnsatz"""
return self.hamiltonian.is_all_terms_commutable()
def get_circuit(self, params):
c = self.init_circuit.copy()
betas = params[:self.step]
gammas = params[self.step:]
for beta, gamma in zip(betas, gammas):
beta *= np.pi
gamma *= 2 * np.pi
for evo in self.time_evolutions:
evo(c, gamma)
for evo1 in self.time_evolutions1:
evo1(c, beta)
for evo2 in self.time_evolutions2:
evo2(c, beta)
return c
h = 4*q(0)+4*q(1)+4*q(2)+4*q(3)+4*q(0)*q(1)+4*q(0)*q(2)+8*q(1)*q(2)+2*q(1)*q(3)+2*q(2)*q(3)
h = h.to_expr().simplify()
runner = Vqe(QaoaQaoaAnsatz(h,4,Circuit().h[0].cx[0,1].x[0].h[2].cx[2,3].x[2]))
result = runner.run()
# get probability
print(result.most_common(12))
print('Result by QAOA')
print(runner.ansatz.get_energy(result.circuit, runner.sampler))
# Hamiltonian to matrix
mat = h.to_matrix()
# Calculate by numpy
print('Result by numpy')
print(np.linalg.eigh(mat)[0][0])
(((1, 0, 0, 1), 0.9774267751895629), ((1, 0, 1, 0), 0.012956207185280066), ((0, 1, 0, 1), 0.007812692781706286), ((0, 1, 1, 0), 0.001804324843443824), ((1, 0, 1, 1), 2.5441534565355386e-31), ((1, 0, 0, 0), 8.011868568650901e-32), ((1, 1, 0, 1), 4.3741462176564236e-32), ((0, 0, 0, 1), 3.242873431465313e-32), ((0, 1, 1, 1), 2.8482886397260487e-32), ((1, 1, 1, 0), 1.1917913216226843e-32), ((0, 0, 1, 0), 9.020324969917482e-33), ((0, 1, 0, 0), 6.374376726091035e-33))
Result by QAOA
こうなりました。大きな確率で制約条件を満たしながら混雑の最も少ないルートが選択されています。
(1,0,0,1)
まとめ
量子回路は長くなりますが、初期状態と時間発展演算子を変更するだけでイジング組合せ最適化問題の弱点である制約条件項をなくすことができました。正解にたどり着く確率が大幅に上がりましたので、かなり革新的な結果と言えるでしょう。以上です。