Help us understand the problem. What is going on with this article?

QAOAセミナー#2

はじめに

量子コンピュータで組合せ最適化問題をときたいという要望は多々あります。実際に速度向上を果たす問題は多くありません。しかしその仕組みを理解するためにQAOAを進めてみます。

ちょっと書きかけ感もありますが、3周目に向けてさらにんブラッシュアップをしてみたいと思います。また、VQEセミナーともかぶるところが多いので、今後は資料を統合する予定です。

こちらの記事も参考にしてください。

「VQEセミナー2周目まとめ」
https://qiita.com/YuichiroMinato/items/d4c15de3b13f7d3acc88

QAOAセミナー#2-1(理論編)

2-1-1 量子断熱計算

量子断熱計算は、状態ベクトルを断続的に変化させることで基底状態をキープしたまま時間発展を行うことができる理論です。

初期状態のハミルトニアンを$H_{start}$として、最終的に求めたい問題のハミルトニアンを$H_{fin}$としたときに、時間$t$と全体のスケジュール$T$から、

H_{temp} = (1-\frac{t}{T})H_{start} + \frac{t}{T}H_{fin}

としたときに$T\rightarrow\infty$とすれば、時間発展で変化させた状態ベクトルが、その瞬間瞬間のハミルトニアンに追従し、固有状態をとり固有値$\lambda$を持つようにすることができます。

H_{temp}\mid \psi \rangle = E_{0temp}\mid \psi \rangle

時間発展計算は、

\mid \psi_{t+1} \rangle = U \mid \psi_t \rangle

となります。課題は基底状態と第一励起状態が最接近する部分ですが、$E_1(t)-E_0(t)$のエネルギー差に注意して計算することによって、基底状態をキープできます。

2-1-2 QAOA

QAOAは上記の量子断熱計算の時間発展計算を変分法で利用できるようにしたもので、実際は時間発展の演算子を変分パラメータで書き換え、

1、時間発展パラメータを導入した回路を作る
2、回路を計算し、ハミルトニアンの期待値を計算する
3、上記のハミルトニアンの期待値を最小にするように時間発展パラメータを古典最適化する

というもので、ほぼVQEと同じです。ansatzがQAOA専用の時間発展回路となっています。実際回路は、

1、量子状態の準備、一般的にはHゲートを全ての量子ビットにかけて|+>状態に。
2、ハミルトニアンの時間発展を任意回転ゲートで実装する
3、上記|+>状態に対応するXゲートの時間発展を適用する
4、上記の2と3を変分パラメータを変更しながら繰り返す

というものです。

https___qiita-user-contents.imgix.net_https%3A%2F%2Fqiita-image-store.s3.ap-northeast-1.amazonaws.com%2F0%2F218694%2Fe10f1843-cc16-cdfe-e4a6-e2fbaab6df9f.png_ixlib=rb-1.2.png

用意するansatzは上記の通り時間発展回路を変分パラメータ化したもので、初期状態に対応した初期ハミルトニアンを含む少し長めの回路です。

作るハミルトニアンは、Zゲートを基本としたイジングモデルと呼ばれる物理モデルにおいて主に組合せ最適化問題を想定した、コスト関数の実装を行います。

イジングハミルトニアンの例は、

cost = -Z(0) - Z(0)*Z(1)

のようになります。

from blueqat import vqe
from blueqat.pauli import Z

cost = -Z(0) - Z(0)*Z(1)
step = 2

result = vqe.Vqe(vqe.QaoaAnsatz(cost, step)).run()
print(result.most_common(12))

2-1-3 Zは1か-1をとる

h = -Z(0) - Z(0)*Z(1)

Zの後ろの数字は、量子ビットの通し番号を表します。0番目と1番目の量子ビットの2つを使っています。また、問題設定で大事なのは、Zの前の係数です。

Z(0)の前は-1
Z(0)*Z(1)の前は-1

がついています。業界では、それぞれ「局所磁場」と「相互作用」とか会社によって呼び名が違いますが、ここではあえて、「バイアス」と「ウェイト」という名前にしてみます。

Zは-1か+1のどちらかをとります。
また、hはより小さい値を取ると正解になります。
上記VQEではhの一番小さい数字を取るように計算されます。

手計算で考えてみると、
Z(0)のバイアスが-1で、Z(0)*Z(1)のウェイトが-1の時、
Z(0)=+1で、
Z(1)=+1の時に

h = -1*1 - 1*1 = -2

で最小になるはずです。これを自動でやってくれます。

Z(0),Z(1)の値で場合わけすると、

Z(0) Z(1) h
-1 -1 0
-1 1 2
1 -1 0
1 1 -2

このZ(0)とZ(1)を探すのが、今回の目的なので、達成されました。

2-1-4 この問題を大きくするだけで社会問題がたくさん解ける

そして実際にはより大きな問題を解きます。

上記問題は役にたたなさそうですが、Z(0)をAさん、Z(1)をBさんに見立てて、Aさんはグループ1に属し、BさんはAさんと同じグループに属するという条件をつけた分類問題と同じです。

ただ、このままでは毎回問題を解くのが大変なので様々な工夫がされています。-1と+1の代わりに0と1を使うなど覚えるべきテクニックがたくさんあります。

2-1-5 定式化は01で

定式化とは、上記の記事でもちょっと触れていますが、組合せ最適化問題を+1と-1の値の組合せで表される式に落とし込みます。ただ、課題があります。物理学で使われるイジングモデルは-1,+1を利用しますが、通常の産業用では0と1を計算として利用します(計算基底)。幸い-1と+1は自動的に変換ができますので、01を使っても組合せ最適化問題としては問題がありません。それがこの8年ほど発達してきました。

01で定式化をするのをQUBOと呼びます。QUBOは量子ビットの01を利用できます。係数はまだ名前が確定していませんが、「バイアス」と「ウェイト(結合荷重)」でいいでしょう。定式化で作るのは「コスト関数」と呼びたいと思います。

定式化は社会問題をQUBO形式で、コスト関数を作ることで実現でき、コスト関数は01の値をとる量子ビットに設定するバイアスとウェイトを設定することで実現できます。

この手法は量子コンピュータに限ったことではないので普通の計算に慣れている人でも受け入れられやすいでしょう。

QAOAセミナー2-2(イジングQUBO実装編)

次は実際に01のバイナリ値を使った実装をしてみます。

2-2-1 QUBOの実装

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)
step = 2

result = vqe.Vqe(vqe.QaoaAnsatz(cost, step)).run()
print(result.most_common(12))

上記で設定する必要があるのはコスト関数です。

cost = -3*q(0)-3*q(1)-2*q(0)*q(1)

明らかにq(0)=1,q(1)=1の時に最小値-3-3-2=-8を取りますね。上記のコスト関数は以前のように-1と+1ではなく、0と1で考えることができるので簡単です。

設定するのは、q(0)とq(1)に-3のバイアスを設定し、q(0)*q(1)に-2のウェイトを設定します。

これをQAOAで解くことで、最小値-8が得られます。

(((1, 1), 0.9922852962030535), ((0, 0), 0.0047331664261775755), ((1, 0), 0.0014907686853841205), ((0, 1), 0.0014907686853841168))

2-2-2 実例:交通最適化

交通最適化は度々出てくるテーマです。

youtubeでは、
交通最適化
https://www.youtube.com/watch?v=JU2omOzuOQg

がわかりやすいです。記事は、

「D-WaveでVW社の交通最適化アプリケーションの実装を解く」
https://qiita.com/YuichiroMinato/items/19dd89e381f14ad24134

「交通最適化問題を量子アニーリングで解く」
https://qiita.com/sychocola1/items/646d2c9504e34ed1589f

今回は、交通最適をときます。やることは、

1、1台の車に3つのルート選択を許容する(古典)
2、現在の選択ルートから混雑状況を計算(古典)
3、混雑度を最小にするようにルート選択を最小化(QAOA)

となります。問題は使い回しをして、上記のリンクから。

「Quantum Computing at Volkswagen:
Traffic Flow Optimization using the D-Wave Quantum Annealer」
引用:https://www.dwavesys.com/sites/default/files/VW.pdf

一つ目は混雑度をコスト関数として、

(q_0+q_1+q_3+q_4)^2+(q_0+q_1+q_3+q_4)^2+(q_0+q_3)^2+(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2+(q_2+q_5)^2\\ 
=2(q_0+q_1+q_3+q_4)^2+2(q_0+q_3)^2+(q_1+q_4)^2+(q_1+q_2+q_4+q_5)^2+3(q_2+q_5)^2\\ 
=4q_0^2+4q_0q_1+8q_0q_3+4q_0q_4+4q_1^2+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2^2 +2q_2q_4+8q_2q_5+4q_3^2+4q_3q_4+4q_4^2+2q_4q_5+4q_5^2\\
=4q_0+4q_0q_1+8q_0q_3+4q_0q_4+4q_1+2q_1q_2+4q_1q_3+8q_1q_4+2q_1q_5+4q_2 +2q_2q_4+8q_2q_5+4q_3+4q_3q_4+4q_4+2q_4q_5+4q_5

3つのルートのうちに1つをえらぶ制約条件をこちらに、

K(q_0+q_1+q_2-1)^2+K(q_3+q_4+q_5-1)^2 = -Kq_0+2Kq_0q_1+2Kq_0q_2-Kq_1+2Kq_1q_2-Kq_2-Kq_3+2Kq_3q_4+2Kq_3q_5-Kq_4+2Kq_4q_5-Kq_5+2K
from blueqat import vqe
from blueqat.pauli import qubo_bit as q

K=2
cost2 = -K*q(0)+2*K*q(0)*q(1)+2*K*q(0)*q(2)-K*q(1)+2*K*q(1)*q(2)-K*q(2)-K*q(3)+2*K*q(3)*q(4)+2*K*q(3)*q(5)-K*q(4)+2*K*q(4)*q(5)-K*q(5)

step = 4

result = vqe.Vqe(vqe.QaoaAnsatz(cost2, step)).run()
print(result.most_common(12))

この二つを繋げてQUBO式とします。

そして、実行します。

from strangeworks.blueqat import StrangeworksBackend #quantumcomputing.comを利用する時に必要

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

K=10
cost1 = 4*q(0)+4*q(0)*q(1)+8*q(0)*q(3)+4*q(0)*q(4)+4*q(1)+2*q(1)*q(2)+4*q(1)*q(3)+8*q(1)*q(4)+2*q(1)*q(5)+4*q(2)+2*q(2)*q(4)+8*q(2)*q(5)+4*q(3)+4*q(3)*q(4)+4*q(4)+2*q(4)*q(5)+4*q(5)
cost2 = -K*q(0)+2*K*q(0)*q(1)+2*K*q(0)*q(2)-K*q(1)+2*K*q(1)*q(2)-K*q(2)-K*q(3)+2*K*q(3)*q(4)+2*K*q(3)*q(5)-K*q(4)+2*K*q(4)*q(5)-K*q(5)

step = 8

result = vqe.Vqe(vqe.QaoaAnsatz(cost1+cost2, step)).run()
print(result.most_common(12))

なかなか答えが出ませんが、、、、

(((0, 1, 0, 0, 0, 1), 0.12986970172628198), ((0, 0, 1, 0, 1, 0), 0.1298697017262819), ((0, 0, 0, 0, 0, 0), 0.08232027492914641), ((0, 1, 0, 0, 1, 0), 0.05463053577530196), ((0, 0, 0, 1, 0, 0), 0.04609411831952029), ((1, 0, 0, 0, 0, 0), 0.04609411831952023), ((0, 0, 1, 1, 0, 1), 0.037740176136034545), ((1, 0, 1, 0, 0, 1), 0.0377401761360345), ((1, 0, 0, 0, 1, 1), 0.03131986281902896), ((0, 1, 1, 1, 0, 0), 0.031319862819028946), ((0, 0, 1, 0, 0, 1), 0.03071276141604334), ((1, 0, 0, 1, 0, 0), 0.02140326399538612))

うまくいけば制約条件を満たす解を与えます。VQEセミナーの中ですでにQAOAに触れていましたので、今後は統合しようと思います。

QAOAセミナー#2-3(中級問題編)

QAOAにも考えるべきパラメータがいくつかあります。

1、ansatz
2、ハミルトニアン
3、古典最適化

2-3-1 最適化アルゴリズムの選択

主にscipyから選べるようです。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

'Nelder-Mead', 'Powell', 'CG', 'BFGS', 'L-BFGS-B', 'TNC', 'COBYLA', 'SLSQP'

この辺りを見れば良さそうです。自分の好きなソルバーを組み込むことができます。

def hyperopt_minimizer(objective, n_params):
    from hyperopt import fmin, Trials, tpe, hp
    trials = Trials()
    best = fmin(objective, [hp.uniform(f'p{i}', 0., 2 * np.pi) for i in range(n_params)],
            algo=tpe.suggest, max_evals=10000, trials=trials, verbose=1)
    return list(best.values())

2-3-2 numbaバックエンドの採用

単純に高速化できます。元は、

from blueqat import vqe
from blueqat.pauli import qubo_bit as q

hamiltonian = q(0)*q(1)*q(2)*q(3)
step = 10

start = time.time()
result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run()
end = time.time()

print(end - start)

こちらはハミルトニアンを作ってvqeでstep = 10で実行します。

32.705730676651
from blueqat import BlueqatGlobalSetting
BlueqatGlobalSetting.set_default_backend('numba')

上記を設定します。

hamiltonian = q(0)*q(1)*q(2)*q(3)
step = 10

start = time.time()
result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run()
end = time.time()

print(end - start)

こちらは、

8.867603778839111

3、4倍ほど速くなりました。

詳しくはこちらをみてください。
「blueqatでいろいろベンチマークとってみる。最適化とかGPUとか。」
https://qiita.com/YuichiroMinato/items/f08420957896b96e6874

2-3-3 swap時間発展

初期状態とmixerを変更することで、ハミルトニアンの制約条件を減らすことができます。

「[QAOA]SWAPの時間発展で|10>と|01>とを入れ替える」
https://qiita.com/gyu-don/items/c51a9e3d5d16a6d5baf6

XX+YYゲートを導入することで、01と10を入れ替える操作を通じてmixerを作れます。

from blueqat import Circuit, pauli, vqe
from blueqat.pauli import qubo_bit as q
from math import pi
import numpy as np

def an(index):
    return 0.5 * pauli.X[index] + 0.5j * pauli.Y[index]

def cr(index):
    return 0.5 * pauli.X[index] - 0.5j * pauli.Y[index]

op = (cr(1) * an(0) + cr(0) * an(1)).to_expr().simplify()
print(op)
x = []
y00 = []
y01 = []
y10 = []
y11 = []
evos = [term.get_time_evolution() for term in op.terms]
for i in range(100):
    c = Circuit().h[0].cx[0, 1].x[1]
    for evo in evos:
        evo(c, i / 100 * pi)
    c.rz(i / 100 * pi)[1]
    for evo in evos:
        evo(c, i / 100 * pi)
    c.rz(i / 100 * pi)[1]
    state = c.run()
    _a, _b, _c, _d = list(np.abs(state) ** 2)
    x.append(i)
    y00.append(_a)
    y10.append(_b)
    y01.append(_c)
    y11.append(_d)

import matplotlib.pyplot as plt
plt.plot(x, y00, label='|00>')
plt.plot(x, y10, label='|10>')
plt.plot(x, y01, label='|01>')
plt.plot(x, y11, label='|11>')
plt.legend()
plt.show()

時間発展を見ても、入れ替え操作が見れます。

まとめ

QAOAは初期状態、ハミルトニアンの設計、古典最適化と様々なチューニングができるので今後も新しいテクを開発するごとに計算精度が上がりそうです。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした