QAOA(Quantum Approximate Optimization Algorithm)

## この記事に関して
ここでは量子アルゴリズムの一つである、QAOAについて、説明します。

## QAOA
エルミート行列(ハミルトニアン)の最小値を求めます。
量子アニーリングの理論をゲート方式で表現したアルゴリズムです。
量子断熱計算を元に行います。

QAOA に入る前に、
[VQE](https://qiita.com/KeiichiroHiga/items/c900d430dc33d4342b24)、[イジングモデル](https://qiita.com/KeiichiroHiga/items/29a8adfce8ca0d8cf146)、[量子アニーリング](https://qiita.com/KeiichiroHiga/items/bbc60ffb958364923b03)を見ておくとわかりやすいです。

量子断熱計算の式を使い、求めるイジングハミルトニアンを以下のように置きます。

``` math
H(t) = \left(1-\frac{t}{\tau}\right)H_1 + \frac{t}{\tau}H_0
```

``` math
H_0 = -\sum_{i<j}^{}J_{ij} Z_i Z_j \ ,\ 
H_1 = \sum_{i}^{}X_i
```
$H_0$ は求めたいハミルトニアンで $H_1$ は $X$ ゲートの和とします。

### 時間発展
量子アニーリングと同じように初期状態を $\left|\phi(0)\right> = \left|+\right>^{\otimes n}$ として、$t$ 秒後の状態を $\left|\phi(t)\right>$ とします。

このときシュレディンガー方程式

``` math
-i \frac{\partial }{\partial t}\left|\phi(t)\right> = H(t)\left|\phi(t)\right>
```
により時間発展は

``` math
\left|\phi(t)\right> = U(t)\left|\phi(0)\right>
```
と前回は記述しました。( $U(0) = I$ )
この $U(t)$ を詳しく見ます。

#### ダイソン級数
$U(t)$ をシュレディンガー方程式に代入して

``` math
-i \frac{\partial }{\partial t}U(t)\left|\phi(0)\right> = H(t)U(t)\left|\phi(0)\right> → 
-i \frac{\partial }{\partial t}U(t) = H(t)U(t)
```
両辺を $[0,t]$ で積分して

``` math
U(t) = -i\int_0^{t}H(s_1)U(s_1)ds_1 + C
```
$C$ は積分定数です。$U(0) = I$ より

``` math
U(t) = I -i\int_0^{t}H(s_1)U(s_1)ds_1
```
$U(s_1)$ に再び上の式を代入します。

``` math
U(t) = I -i\int_0^{t}H(s_1)\left(
I -i\int_0^{s_1} H(s_2)U(s_2)ds_2
\right)ds_1
```

``` math
= I -i\int_0^{t}H(s_1)ds_1 +
(-i)^2\int_0^t\int_0^{s_1} H(s_1)H(s_2)U(s_2)ds_1ds_2
```
以下繰り返すと $U(t)$ は

``` math
U(t) = I + \sum_{k=1}^{\infty}(-i)^k
\int_0^t\int_0^{s_1} \cdots \int_0^{s_{k-1}}
H(s_1)H(s_2) \cdots H(s_k)
ds_1 ds_2 \cdots ds_k
```
これをダイソン級数と言います。

#### 時間順序積 (T積)
上の式をユニタリ行列にするために指数行列で書き直します。
ここで時間順序積 $T$ を考えます。

``` math
T[H(s_1)H(s_2)] = \left \{ \begin{array}{cc}
H(s_1)H(s_2)& (s_1>s_2) \\ H(s_2)H(s_1)& (s_2>s_1)
\end{array} \right.
```
これを用いると $t > s_1 > \cdots > s_{N-1} > 0$ より $s_l$ の組み合わせを考えると

``` math
k!\int_0^t\int_0^{s_1} \cdots \int_0^{s_{k-1}}
H(s_1)H(s_2) \cdots H(s_k)
ds_1 ds_2 \cdots ds_k \\
=\int_0^t \cdots \int_0^t T[H(s_1)H(s_2) \cdots H(s_k)]ds_1 ds_2 \cdots ds_k
```
$s_l$ は $0 < s_l <s_{l-1}$ でしか定義されていないので積分範囲を $t$ まで拡張しても値は変わりません。
よって $U(t)$ は $T$ を用いると、

``` math
U(t) = I + \sum_{k=1}^{\infty}\frac{(-i)^k}{k!}
\int_0^{t} \cdots \int_0^{t}
T[H(s_1)H(s_2) \cdots H(s_k)]
ds_1 ds_2 \cdots ds_k
```

``` math
=T\left[\sum_{k=0}^{\infty}\frac{(-i)^k}{k!}\left( \int_0^{t}H(s)ds \right)^k\right] = T\left[e^{-i\int_0^t H(s)ds}\right]
```
ここで指数に乗っている積分を近似させます。
区分求積法を今回は用います。

区間 $I = [0,t]$ を $N$ 分割した区間 $I_k = [t_k, t_{k+1}]$ を考えます。$(0 = t_0 < t_1 < \cdots < t_N = t)$

区分求積では各区間 $s \in I_k$ で $H(s) = H(t_k)$ と定数と考えるので、

``` math
-i\int_0^t H(s)ds = -i\left(
H(0)\int_0^{t_1}ds_1 + H(t_1)\int_{t_1}^{t_2}ds_2 + \cdots + 
H(t_{N-1})\int_{t_{N-1}}^{t}ds_N
\right) +
O\left(\frac{1}{N}\right)
```
誤差は分割数 $N$ に反比例します。
$t_k-t_{k-1} = \Delta t_k$ と置けば

``` math
-i\int_0^t H(s)ds \approx 
-i\left(H(0)\Delta t_1 + H(t_1)\Delta t_2 + \cdots +
H(t_{N-1})\Delta t_N \right)
```
以上から $U(t)$ は次のように近似できる。

```  math
U(t) \approx T\left[e^{-iH(0)\Delta t_1 -iH(t_1)\Delta t_2 - \cdots -i
H(t_{N-1})\Delta t_N}\right]
```
ここから更に指数行列を分けるために近似します。
$H(t_k)$ と $H(t_l)$ について $H(t_k)H(t_l)-H(t_l)H(t_k) = 0$ と仮定すれば、$T$ 積を考慮して

``` math
T\left[e^{-iH(0)\Delta t_1 -iH(t_1)\Delta t_2 - \cdots -i
H(t_{N-1})\Delta t_N}\right] \approx
e^{-iH(t_{N-1})\Delta t_N}e^{-iH(t_{N-2})\Delta t_{N-1}} \cdots e^{-iH(0)\Delta t_1}\\
(0 = t_0 < t_1 < \cdots < t_N = t)
```
これによる誤差は $t_k-t_l$ に比例して $\tau$ に反比例します。

よって $e^{-i H(t_k)\Delta t_{k+1}}$ に関して

``` math
e^{-i H(t_k)\Delta t_{k+1}} = e^{-i \left(\left(1-\frac{t_k}{\tau}\right)H_1 +\frac{t_k}{\tau}H_0 \right)\Delta t_{k+1}} =
e^{-i\left(1-\frac{t_k}{\tau}\right)H_1 \Delta t_{k+1}-
\frac{it_k}{\tau}H_0 \Delta t_{k+1}}
```
再び指数部分を分解させます。次は Lie – Trotter 積公式を使います。
$e^{-i H(t_k)\Delta t_{k+1}}$ にこの公式を用いると

``` math
e^{-i H(t_k)\Delta t_{k+1}} \approx
\left(e^{-\frac{i}{m}
\left(1-\frac{t_k}{\tau}\right)H_1 \Delta t_{k+1}}
e^{-\frac{it_k}{\tau m}H_0 \Delta t_{k+1}}\right)^m
```
ここで $\beta_k, \gamma_k$ を以下で定義する。

```  math
\beta_k = \frac{1}{m}\left(1-\frac{t_k}{\tau}\right)\Delta t_{k+1} \ , \ \gamma_k = \frac{t_k}{\tau m}\Delta t_{k+1}
```
これを用いると

``` math
e^{-i H(t_k)\Delta t_{k+1}} \approx
e^{-i\beta_k H_1}e^{-i\gamma_k H_0}e^{-i\beta_k H_1}e^{-i\gamma_k H_0}\cdots
e^{-i\beta_k H_1}e^{-i\gamma_k H_0}
```
$H_0$ と $H_1$ の指数行列が交互に出てくることがわかります。

#### β と γ の意味
$\beta_k$ と $\gamma_k$ については和をとればわかりやすいです。

``` math
\beta_k + \gamma_k = \frac{\Delta t_{k+1}}{m}
```
和は各時間差を更に $m$ 等分に分割していることがわかります。
またトロッター数 $m$ は別で定めるので、各 $\Delta t_{k+1}$ の大きさを決めていることになります。

$t_k$ に関しては初期値が $0$ で徐々に $\tau$ に近づくので、 $\gamma_k$ が始め大きく $\tau$ に近づくにつれて $\beta_k$ の方が大きくなるのがわかります。

### 続き
以上の流れから $H_0$ の基底状態を $\left|\psi_0\right>$ とすれば $t = \tau$ のとき、断熱定理より

``` math
\left|\phi(\tau)\right> = U(\tau)\left|\phi(0)\right> \approx \prod_{k=0}^{N-1}e^{-i H(t_k)\Delta t_{k+1}}\left|\phi(0)\right>
```

``` math
\approx \prod_{k=0}^{N-1} e^{-i\beta_k H_1}e^{-i\gamma_k H_0}e^{-i\beta_k H_1}e^{-i\gamma_k H_0}\cdots e^{-i\beta_k H_1}e^{-i\gamma_k H_0}\left|\phi(0)\right> \approx \left|\psi_0\right>
```
基底状態に近づくことがわかります。
(ここでの積は時間順序積に従うので左側からかけていく)

### アルゴリズム
$H_0$ の最小固有値 $E_0$ を求めます。
$E_0$ は上に述べた基底状態 $\left|\psi_0\right>$ を用いて以下のように表せます。

``` math
E_0 = \left<\psi_0\right|H_0\left|\psi_0\right>
```

$\left|\psi_0\right>$ を作るために $\beta_k,\gamma_k$ を考えます。
時間の分割数を $N$ とします。
今回は $\beta = (\beta_0, \cdots ,\beta_{N-1}),\gamma = (\gamma_0, \cdots ,\gamma_{N-1})$ と置き、

``` math
\left|\phi(\tau)\right> \approx \left|\beta,\gamma\right>
```
とします。

#### 量子変分原理
VQE で使われる量子変分原理を考えます。
任意の状態ベクトル $\left|\Psi\right>$ に関して以下が成り立ちます。

``` math
\left<\Psi\right|H_0\left|\Psi\right> ≧ E_0
```
今回は $\left|\beta,\gamma\right>$ に関して考えると

``` math
\left<\beta,\gamma\right|H_0\left|\beta,\gamma\right> ≧ E_0
```
となることがわかります。

#### ベイズ最適化
$\beta, \gamma$ をうまくとりたい所ですが規則性がありません。
したがって推測して確率的にパラメータを定めます。

このように確率的に適当なパラメータを取り、最適化する手法をベイズ最適化と言います。

以上から $\left<\beta,\gamma\right|H_0\left|\beta,\gamma\right>$ をベイズ最適でパラメータを変えて行くと最終的に

``` math
\left<\beta,\gamma\right|H_0\left|\beta,\gamma\right> \xrightarrow{Bayesian Opt}
\left<\psi_0\right|H_0\left|\psi_0\right> = E_0
```
と説明ができました。

## 量子コンピュータでの実装
$H_0\left|\beta,\gamma\right>$ を量子ゲートで実装します。
$H_0$ を詳しく書きます。

``` math
H_0\left|\beta,\gamma\right> = -\sum_{i<j}^{}J_{ij} Z_i Z_j\left|\beta,\gamma\right>
```
上の式から $Z_i Z_j\left|\beta,\gamma\right>$ を考えます。
$\left|\beta,\gamma\right>$ について、上で述べた時間発展より

``` math
\left|\beta,\gamma\right> =
\prod_{k=0}^{N-1} e^{-i\beta_k H_1}e^{-i\gamma_k H_0}e^{-i\beta_k H_1}e^{-i\gamma_k H_0}\cdots e^{-i\beta_k H_1}e^{-i\gamma_k H_0}\left|\phi(0)\right>
```
まず $\left|\phi(0)\right>$ は

``` math
\left|\phi(0)\right> = \left|+\right>^{\otimes n} = \bigotimes_{}^{n}H\left|0\right>
```
から全体に $H$ ゲートを施せばいいことがわかります。

次に $e^{-i\beta_k H_1}$ について

``` math
e^{-i\beta_k H_1} = e^{-i\beta_k\sum_{i}^{}X_i} = \prod_{i}^{}e^{-i\beta_k X_i}
```
ここで $X = HZH, e^{-iZt} = Rz(-2t)$ より

``` math
e^{-i\beta_k X_i} = H_i e^{-i\beta_k Z_i} H_i = H_i Rz_i(-2\beta_k) H_i
```
以上から

``` math
e^{-i\beta_k H_1} = \prod_{i}^{}H_i Rz_i(-2\beta_k) H_i
```
$H$ ゲートと $Rz$ ゲートで表現できました。

最後に $e^{-i\gamma_k H_0}$ について

``` math
e^{-i\gamma_k H_0} = e^{-i\gamma_k\sum_{i<j}^{}J_{ij}Z_i Z_j} = \prod_{i<j}^{}e^{-i\gamma_k J_{ij}Z_i Z_j}
```
ここで $e^{-iZ_i Z_j t} = CX_{i,j} Rz_j(-2t) CX_{i,j}$ より
( $CX_{i,j}$ は $i$ から $j$ に向かう $CX$ ゲート)

``` math
e^{-i\gamma_k J_{ij}Z_i Z_j} = CX_{i,j} Rz_j(-2\gamma_k J_{ij}) CX_{i,j}
```
以上から

``` math
e^{-i\gamma_k H_0} = \prod_{i<j}^{} CX_{i,j} Rz_j(-2\gamma_k J_{ij}) CX_{i,j}
```
$CX$ ゲートと $Rz$ ゲートで表現できました。
これで $\left|\beta,\gamma\right>$ をゲートで作ることができました。

$Z_iZ_j$ はそのまま $Z$ ゲートを施せばいいので、
これらの流れから各 $Z_i Z_j\left|\beta,\gamma\right>$ を計算して古典計算機で和をとれば
$H_0\left|\beta,\gamma\right>$ が実装できることがわかります。

## 補足
[Rigetti の論文](https://arxiv.org/pdf/1712.05771.pdf)を見てみます。
ここでは $\left|\beta,\gamma\right>$ の作り方として、トロッター数を $m = 1$ にして、更に時間の分割数を $N = 1$ にして行なっています。

求めた $\left|\beta,\gamma\right>$ を使って古典計算機で $\left<\beta,\gamma\right|H_0\left|\beta,\gamma\right>$ を求めベイズ最適をします。
ベイズ最適で得た $\beta, \gamma$ をまた量子コンピュータに入れこれを繰り返しています。

全体図は以下の通り

<img src="https://qiita-user-contents.imgix.net/https%3A%2F%2Fqiita-image-store.s3.amazonaws.com%2F0%2F310258%2F2beb8cd9-36f6-4b20-3f5e-7a2bfc65398e.png?ixlib=rb-1.2.2&auto=format&gif-q=60&q=75&s=dcb0b45ce4bcab80277620f7d116b773" width="80%">

$m,N = 1$ より $\beta, \gamma$ が一回しかループしていないことがわかります。
ループは少ないが量子変分原理から操作を繰り返し行うと最小値が求まります。

## まとめ
今回は QAOA について説明しました。
次回は QUBO について説明しようと思います。

## 参考文献
A Quantum Approximate Optimization Algorithm
https://arxiv.org/pdf/1411.4028.pdf

Quantum Approximate Optimazation Algorithm (QAOA): 量子近似最適化アルゴリズム
https://dojo.qulacs.org/ja/latest/notebooks/5.3_quantum_approximate_optimazation_algorithm.html#QAOA%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0%E3%81%AE%E6%89%8B%E9%A0%86

QAOA - ゲート式量子コンピューターで最適化問題を解く近似アルゴリズム
https://qiita.com/snhrhdt/items/ae55a94b25c06142528a

相互作用表示と時間発展演算子
http://park.itc.u-tokyo.ac.jp/kato-yusuke-lab/nagai/note_070106_in.pdf

時間依存の摂動 - Magnus展開と平均ハミルトニアン
https://whyitsso.net/physics/quantum_mechanics/perturbation5.html

VQE algorithm (Variational Quantum Eigensolver)
https://qiita.com/KeiichiroHiga/items/c900d430dc33d4342b24#%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0

イジングモデル (Ising model)
https://qiita.com/KeiichiroHiga/items/29a8adfce8ca0d8cf146

量子アニーリング (Quantum Anneling)
https://qiita.com/KeiichiroHiga/items/bbc60ffb958364923b03