インターン講座３：QAOA

今回はインターン３回目です。ようやく量子計算にはいります。
なるべくサクッと行きます。ちなみに今回はQAOAを行いますが、以前行ったQAOA説明コンテストの応募内容がとても高度でした。
こちらを参考にしてもらえるとさらに理解が深まります。

量子コンピュータの応用 - QAOAの仕組みとこれを応用した組合せ最適化問題のプログラミング(blueqat編)

量子状態と状態ベクトル

量子コンピュータで利用されるデータの扱いは今のコンピュータと違って、1量子ビットでは球を使って表現されます。量子ビットは|0>状態と|1>状態を持っていてベクトルで表現できます。

\mid 0 \rangle = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \mid 1 \rangle = \begin{bmatrix} 0 \\ 1 \end{bmatrix}

また、任意の量子ビットは重ね合わせて量子状態を表現できます。

\mid \psi \rangle = \alpha \mid 0 \rangle + \beta \mid 1 \rangle = \alpha \begin{bmatrix} 1 \\ 0 \end{bmatrix} + \beta \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} \alpha \\ \beta \end{bmatrix}

また、下記条件を満たします。

|\alpha|^2+|\beta|^2 = 1

$\mid \psi \rangle$ は状態ベクトル(波動関数）と呼ばれ量子ビットの量子状態を表します。 $\alpha$ と $\beta$ は複素数で「確率振幅」と呼ばれ、二乗するとその対応する計算基底の出現確率になり、 $|\alpha|^2$ や $|\beta|^2$ は「測定」をすることで取得できます。

量子状態の極座標表示

ブロッホ球は二つの直交する純粋状態の重ね合わせで表現できる量子状態を単位球面上に表す表記法。

vqe2

\mid \psi \rangle = \begin{bmatrix} cos\theta & -sin\theta e^{-i\phi}\\ sin\theta e^{i\phi} & cos\theta \end{bmatrix} \begin{bmatrix} 1\\0 \end{bmatrix} = \begin{bmatrix} cos\theta \\ sin\theta e^{i\phi} \end{bmatrix}

のように角度を利用したパラメータ表記もできます。量子データは角度で表現されるということを覚えてください。

状態ベクトルの操作と量子ゲート

状態ベクトルを操作するには、量子ゲートを利用します。変分アルゴリズムでは基本のパウリゲート、任意回転ゲート、アダマールゲートを主に利用します。

パウリゲートはXYZの３種類があり、それぞれブロッホ球でのXYZ軸周りでの180度の回転を意味し、対応する行列は、

X= \begin{bmatrix} 0&1 \\ 1&0 \end{bmatrix} , Y= \begin{bmatrix} 0&-i \\ i&0 \end{bmatrix} , Z= \begin{bmatrix} 1&0 \\ 0&-1 \end{bmatrix}

となります。また、XYZ軸周りに任意回転にしたのが、RX,RY,RZゲートで、

Rx(\theta) = \begin{bmatrix} \cos\left(\frac{\theta}{2}\right) & -i\sin\left(\frac{\theta}{2}\right)\\ -i\sin\left(\frac{\theta}{2}\right) & \cos\left(\frac{\theta}{2}\right) \end{bmatrix}, Ry(\theta) = \begin{bmatrix} \cos\left(\frac{\theta}{2}\right) & -\sin\left(\frac{\theta}{2}\right)\\ \sin\left(\frac{\theta}{2}\right) & \cos\left(\frac{\theta}{2}\right) \end{bmatrix}, Rz(\theta) = \begin{bmatrix} e^{-i\frac{\theta}{2}} & 0\\ 0 & e^{i\frac{\theta}{2}} \end{bmatrix}

となります。状態ベクトルにこれらの量子ゲートを作用させることで量子状態を変化させることができます。

X \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0&1 \\ 1&0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}

アダマールゲートはXZ軸45度周りに回転させる便利なゲートです。

H= \frac{1}{\sqrt{2}} \begin{bmatrix} 1&1 \\ 1&-1 \end{bmatrix}

固有値・固有ベクトル問題

ある行列に対する固有値と固有ベクトルを見つけることができれば様々な問題を解くことができます。与えられた行列 $H$ （ハミルトニアンと呼ばれます）に対して、固有値を $E_0$ 、固有ベクトルを $\mid \psi \rangle$ としたときに、

H \mid \psi \rangle = E_0 \mid \psi \rangle

を満たすような、 $E_0$ と $\mid \psi \rangle$ を見つけるのが目的です。 $\mid \psi \rangle$ は状態ベクトル（波動関数）と呼ばれます。

量子コンピュータの問題のほとんどがこの固有値問題に対応しています。行列Hに問題を設定し、固有値Eを探します。

量子変分原理

固有値の探し方はNISQでは主に最適化を使います。インターンでは二回目に学びました。

\langle \psi (\theta) \mid H \mid \psi (\theta) \rangle \geq E_0

固有値を直接求めるのは厳しいので、間接的に固有値の期待値にできるだけ近づける最適化計算に落とし込みます。ハミルトニアン $H$ は与えられた式なので、私たちにできるのは状態ベクトル $\mid \psi (\theta) \rangle$ をできるだけ固有状態に近づける状態ベクトルを角度の最適化計算として探していくのが今回の目的です。状態ベクトルの探索はansatz（アンザッツ）と呼ばれる回路を作ることで実現します。

ハミルトニアンとその期待値

ハミルトニアン $H$ はパウリゲートXYZと単位行列Iから構成されるエルミート行列です。ある量子状態 $\mid \psi\rangle$ が与えられた場合の、ハミルトニアン $H$ の期待値 $E$ は、上記の式の両方に $\langle \psi \mid$ を左からかけてみると、

\langle \psi \mid H \mid \psi \rangle = \langle \psi \mid E \mid \psi \rangle\\ \langle \psi \mid H \mid \psi \rangle = E \langle \psi \mid \psi \rangle\\ \langle \psi \mid H \mid \psi \rangle = E

のように計算できます。1量子ビットの例として、任意の量子状態 $\mid \psi \rangle$ と任意のハミルトニアン $H$ から、

\langle \psi \mid H \mid \psi \rangle = \begin{bmatrix} \alpha^* & \beta^* \end{bmatrix} \begin{bmatrix} a&b\\ c&d \end{bmatrix} \begin{bmatrix} \alpha\\ \beta \end{bmatrix}

のように計算できます。

ハミルトニアンZの期待値を求める

任意の状態ベクトルにおけるハミルトニアンが $H = Z$ の場合の期待値は、

\langle \psi \mid H \mid \psi \rangle = \begin{bmatrix} \alpha^* & \beta^* \end{bmatrix} \begin{bmatrix} 1&0\\ 0&-1 \end{bmatrix} \begin{bmatrix} \alpha\\ \beta \end{bmatrix} = |\alpha|^2 - |\beta|^2

となります。 $|\alpha|^2$ は0がでる確率、 $|\beta|^2$ は1が出る確率に対応します。

ハミルトニアンが $H = X$ の場合のハミルトニアン $H$ の期待値は、

\langle \psi \mid H \mid \psi \rangle = \begin{bmatrix} \alpha^* & \beta^* \end{bmatrix} \begin{bmatrix} 0&1\\ 1&0 \end{bmatrix} \begin{bmatrix} \alpha\\ \beta \end{bmatrix} = \alpha^* \beta + \alpha \beta^*

となりますが、通常の測定ではこの値を求めることができませんので、ハミルトニアンの変形を通じて解を得ます。

ハミルトニアンの変形

ハミルトニアン $H$ が $Z$ の場合には、期待値の計算ができますが、XやYの場合には求めることができません。そこでハミルトニアンの変形を利用します。

X = HZH

より、

\langle \psi \mid X \mid \psi \rangle \\ = \langle \psi \mid HZH \mid \psi \rangle\\ = \langle \psi' \mid Z \mid \psi' \rangle

のように状態ベクトルを変形することで対応できます。Yの場合は、

Y = RX(-\pi/2) Z RX(\pi/2)

と、なり同様に、

\langle \psi \mid Y \mid \psi \rangle \\ = \langle \psi \mid RX(-\pi/2) Z RX(\pi/2) \mid \psi \rangle\\ = \langle \psi'' \mid Z \mid \psi'' \rangle

とすることができます。これらは、量子回路の測定の直前に上記のHやRXゲートを入れて測定を行うことで期待値を計算することができます。

線形結合

実際にはハミルトニアン $H$ はエルミート行列で与えられ、ユニタリ行列の和の形で与えられ、下記のように分解できます。

\langle \psi \mid aH_1 + bH_2 \mid \psi \rangle \\ = \langle \psi \mid aH_1 \mid \psi \rangle + \langle \psi \mid bH_2 \mid \psi \rangle \\ = a\langle \psi \mid H_1 \mid \psi \rangle + b\langle \psi \mid H_2 \mid \psi \rangle

これを使えば、複雑なハミルトニアンの場合でも、

H = 1.2 X_0 Z_2 + 2.5 Z_0 X_1 Y_2 - 3.4 Z_2 X_1

期待値は、

\langle \psi \mid 1.2 X_0 Z_2 + 2.5 Z_0 X_1 Y_2 - 3.4 Z_2 X_1 \mid \psi \rangle\\ = 1.2\langle \psi \mid X_0 Z_2 \mid \psi \rangle + 2.5 \langle \psi \mid Z_0 X_1 Y_2\mid \psi \rangle - 3.4 \langle \psi \mid Z_2 X_1 \mid \psi \rangle

のように個別に計算をして和を求めることで計算できます。

量子断熱計算

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

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

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

としたときに $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^{-iHt} \mid \psi_t \rangle

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

QAOA回路

QAOAは上記の量子断熱計算の時間発展計算をansatzとして変分アルゴリズムに適用したものです。

|0> --H-- --RZZ--RZ-- --RX-- --RZZ--RZ-- --RX--
            |                  |
|0> --H-- --RZZ--RZ-- --RX-- --RZZ--RZ-- --RX--

全体構成を大きく分けて、

1.mixerを選ぶと、初期状態が決まる
2.定式化をQUBOもしくはZで作る

です。2番目の定式化は自動的に時間発展されます。

一番左は初期固有状態|+>を作っています。これに対応するハミルトニアンXは時間発展の形でRxで出てきます。通常はmixerはデフォルトではXを選び、初期状態はそれに対応する|+>を作るために、 $H^ {\otimes N}$ が利用されることが多いです。

また、CX-Rz-CXは問題設定のハミルトニアンのウェイトに対応し、次のRzはハミルトニアンのバイアスに対応しています。

blueqatのQaoaAnsatzはハミルトニアンとステップ分割数を決めれば、自動的にこの時間発展のansatzを構成してくれます。ライブラリ側で適切な計算をしてくれているので、私たちがするのはハミルトニアンの定式化とmixerの選択です。簡単な定式化を解いてみましょう。
##　イジング定式化
イジングモデルと呼ばれる物理モデルを利用して組合せ最適化問題を解きます。イジングモデルでは通常スピンの上下を利用して問題を定式化します。定式化では演算子としてZを利用します。Zは+1と-1を期待値としてもとめられます。

H = -Z_0 -Z_0*Z_1

##　QUBO
つぎに、上記ハミルトニアンZの定式化から派生して、0と1で定式化をするQUBOを確認します。QUBOはZの定式化を書き換えた形になりますが、blueqatでは自動変換してくれるツールがあります。

定式化は、

H = 2q_1-4q_0q_1

としてみたいと思います。

from blueqat import Circuit
from blueqat.utils import qaoa
from blueqat.pauli import qubo_bit as q
from blueqat.pauli import X,Y,Z,I

hamiltonian = 2*q(1)-4*q(0)*q(1)
step = 1

result = qaoa(hamiltonian, step)
result.circuit.run(shots=100)

Counter({'11': 47, '00': 36, '10': 17})