Qulacsで量子位相推定を試す

# はじめに
この記事では、[前回](https://qiita.com/meguminlove/items/6eaaf8a6cb0ea10499cd)と同じく[Qulacs](http://qulacs.org/index.html)を使って、量子位相推定(<i>Quantum Phase Estimation</i>, QPE)を試します。

QPEとは、ユニタリー行列の固有値を量子ビット上に得るアルゴリズムです。より正確には、ユニタリー行列の固有値は絶対値が1なので、その偏角(argument)が得られれば、固有値を得たことになります。QPEはその2進展開を量子ビット上に実現します。

今回はエルミート行列 $A\in\mathbb{R}^{d\times d}$ に対して、$e^{-itA}$ の固有値を求める、というのを試したいと思います。ここで、

```math
u\colon \mathbb{R}\ni t\mapsto e^{-itA}u(0)\in\mathbb{R}^d
```
は、シュレディンガー方程式

```math
i\frac{du}{dt}=Au, \quad t\in\mathbb{R}
```
を満たします。

さて、なぜこれをやるのかというと、$A$ の固有展開

```math
A=\sum_j \lambda_j P_j
```
に対し、

```math
e^{-itA}=\sum_j e^{-it\lambda_j}P_j
```
であり[^1]、$e^{-itA}$ の固有値の偏角を量子ビット上に得るということは、それはすなわち $A$ の固有値そのものを得るということになるからです。このことは、線型方程式を量子的に解くアルゴリズムとして知られている[HHLアルゴリズム](https://en.wikipedia.org/wiki/Quantum_algorithm_for_linear_systems_of_equations)において重要な役割を果たします。それも見据えてこの部分を動かしてみようというのが目的です。

# 今回の設定
今回は簡単にして、2次元のエルミート行列

```math
A=
\begin{pmatrix}
0 & -i\frac{3\pi}{4} \\
i\frac{3\pi}{4} & 0
\end{pmatrix}
```
で考えます。固有値と固有ベクトルは対応順にそれぞれ $\pm\frac{3\pi}{4}$、$\frac{1}{\sqrt{2}}(1,\pm i)^\text{T}=\frac{|0\rangle\pm i|1\rangle}{\sqrt{2}}$ です。ただし、通常よく用いられるように、

```math
|0\rangle := 
\begin{pmatrix}
1 \\
0
\end{pmatrix}, \quad
|1\rangle := 
\begin{pmatrix}
0 \\
1
\end{pmatrix}
```
です。このとき、

```math
U:=\left. e^{-itA} \right|_{t=1}=\frac{1}{\sqrt{2}}
\begin{pmatrix}
-1 & -1 \\
1 & -1
\end{pmatrix}
```
に対してQPEを行います。

$U$ の固有値は $e^{-i3\pi/4}=e^{i2\pi 0.101}$ と $e^{i3\pi/4}=e^{i2\pi 0.011}$ です。ただし、小数は2進小数です。よって、レジスタービット数を3とすれば、QPEでそれぞれ $|101\rangle$ と $|011\rangle$が(確率1で)得られるはずです。

# Qulacsで使える基本的なゲートについて
ここでは、Qulacsの基本的なゲートについて、備忘録的に記載しておきたいと思います。

- アダマールゲート

```math
\begin{align*}
H &:= \frac{1}{\sqrt{2}}
    \begin{pmatrix}
    1 & 1 \\
    1 & -1
    \end{pmatrix}
\end{align*}
```

- パウリゲート

```math
\begin{align*}
X &:= 
    \begin{pmatrix}
    0 & 1 \\
    1 & 0
    \end{pmatrix}, \quad
Y := 
    \begin{pmatrix}
    0 & -i \\
    i & 0
    \end{pmatrix}, \quad
Z := 
    \begin{pmatrix}
    1 & 0 \\
    0 & -1
    \end{pmatrix}
\end{align*}
```

- ブロッホ球上の回転

Qulacsの動作を見ていると、通常[^2]の $R_X,R_Y,R_Z$ ゲートとは異なり、角度 $\theta$ に対して逆に定義されているようです。

```math
\begin{align*}
R_X(\theta)
&:= e^{i\frac{\theta}{2}X} =
    \begin{pmatrix}
    \cos\frac{\theta}{2} & i\sin\frac{\theta}{2} \\
    i\sin\frac{\theta}{2} & \cos\frac{\theta}{2}
    \end{pmatrix}, \\[.5em]
R_Y(\theta) 
&:= e^{i\frac{\theta}{2}Y} = 
    \begin{pmatrix}
    \cos\frac{\theta}{2} & \sin\frac{\theta}{2} \\
    -\sin\frac{\theta}{2} & \cos\frac{\theta}{2}
    \end{pmatrix}, \\[.5em]
R_Z(\theta) 
&:= e^{i\frac{\theta}{2}Z} = 
    \begin{pmatrix}
    e^{i\frac{\theta}{2}} & 0 \\
    0 & e^{-i\frac{\theta}{2}}
    \end{pmatrix}
\end{align*}
```

# QPE回路図
量子回路図は次のようになります。
![回路図1.jpg](https://qiita-image-store.s3.amazonaws.com/0/330467/da917b63-0972-27cd-27ff-9bc004954048.jpeg)
I-QFTは、逆量子フーリエ変換 (Inverse Quantum Fourier Transform) です。
ただし、

```math
\begin{align*}
|\psi\rangle_\pm = \frac{|0\rangle\pm i|1\rangle}{\sqrt{2}}
\end{align*}
```
であり、

```math
\begin{align*}
R^*_1 &= e^{-i\frac{\pi}{4}}R_Z\left( \frac{\pi}{2} \right), \quad R^*_2 = e^{-i\frac{\pi}{8}}R_Z\left( \frac{\pi}{4} \right)
\end{align*}
```
です。

上の回路図では 初期状態として $|\psi\rangle_\pm$ を描きましたが、プログラムではこれを用意してやらねばなりません。$R_X(\pm\pi/2)|0\rangle=|\psi\rangle_\pm$なので、初期状態を $|0000\rangle$ として、$R_X(\pm\pi/2)$ ゲートを入れて $|000\rangle \otimes |\psi\rangle_\pm$ を作ります。

# プログラム
勉強のため、for文等使わず、回路図の通り愚直に書いています。

```python
import numpy as np
from qulacs import QuantumState
from qulacs import QuantumCircuit
from qulacs import Observable
from qulacs.gate import DenseMatrix
from qulacs.gate import H
from qulacs.gate import RX
from qulacs.gate import RZ
from qulacs.gate import merge
from qulacs.gate import to_matrix_gate

# 量子状態の設定
# 量子ビット数を3(固有値用)+1(固有ベクトル用)=4で
n_qbits = 4
state = QuantumState(n_qbits)
state.set_zero_state()

# 量子回路
circuit = QuantumCircuit(n_qbits)

# A = [[0,-i3π/4],[i3π/4,0]]
# からユニタリー U = e^{-iA}とする
sqr2 = np.sqrt(2)
U = [
    [-1.0 / sqr2, -1.0 / sqr2],
    [+1.0 / sqr2, -1.0 / sqr2]
]

# ========== 前半QPE部 ==========
# |000>⊗|ψ>を準備。|ψ>はUの固有状態。
# H⊗H⊗H⊗RX(±π/2)|0000>=(H⊗H⊗H|000>)⊗|ψ>
# *** U=e^{-iA}の固有値と固有ベクトル情報 ***
# 固有状態:RX(+π/2)|0>=(|0>+i|1>)/√2, 固有値:exp(-i3π/4)=exp(i5π/4)→0.101
# 固有状態:RX(-π/2)|0>=(|0>-i|1>)/√2, 固有値:exp(+i3π/4)→0.011
# * RX(Θ)=exp(iΘX/2)
gate_list = [H(1), H(2), H(3)]
angle = np.pi / 2.0  # (☆) ±π/2
gate_list.append(RX(0, angle))

# controlled-U達のゲート
ctrl_U1_gate = DenseMatrix(0, U)
ctrl_U1_gate.add_control_qubit(3, 1)
ctrl_U2_gate = DenseMatrix(0, U)
ctrl_U2_gate.add_control_qubit(2, 1)
ctrl_U3_gate = DenseMatrix(0, U)
ctrl_U3_gate.add_control_qubit(1, 1)

# controlled gatesをリストに登録
# 1回
gate_list.append(ctrl_U1_gate)
# 2回
gate_list.append(ctrl_U2_gate)
gate_list.append(ctrl_U2_gate)
# 4回
gate_list.append(ctrl_U3_gate)
gate_list.append(ctrl_U3_gate)
gate_list.append(ctrl_U3_gate)
gate_list.append(ctrl_U3_gate)


# ========== 後半、I-QFT部 ==========
# 量子ビット1にHゲート追加
gate_list.append(H(1))
# 量子ビット2に制御R_1^*を追加、Hゲート追加
ctrl_adjoint_R1_gate = RZ(2, np.pi / 2.0)
ctrl_adjoint_R1_gate = to_matrix_gate(ctrl_adjoint_R1_gate)
ctrl_adjoint_R1_gate.multiply_scalar(np.exp(-1j * np.pi / 4.0))
ctrl_adjoint_R1_gate.add_control_qubit(1, 1)
gate_list.append(ctrl_adjoint_R1_gate)
gate_list.append(H(2))

# 量子ビット3に制御R_2^*、制御R_3^*、Hゲート追加
ctrl_adjoint_R2_gate = RZ(3, np.pi / 4.0)
ctrl_adjoint_R2_gate = to_matrix_gate(ctrl_adjoint_R2_gate)
ctrl_adjoint_R2_gate.multiply_scalar(np.exp(-1j * np.pi / 8.0))
ctrl_adjoint_R2_gate.add_control_qubit(1, 1)
ctrl_adjoint_R3_gate = RZ(3, np.pi / 2.0)
ctrl_adjoint_R3_gate = to_matrix_gate(ctrl_adjoint_R3_gate)
ctrl_adjoint_R3_gate.multiply_scalar(np.exp(-1j * np.pi / 4.0))
ctrl_adjoint_R3_gate.add_control_qubit(2, 1)
gate_list.append(ctrl_adjoint_R2_gate)
gate_list.append(ctrl_adjoint_R3_gate)
gate_list.append(H(3))

# ゲートをマージして回路に追加
merged_gate = merge(gate_list)
circuit.add_gate(merged_gate)

# 状態の更新
circuit.update_quantum_state(state)
# 状態のベクトル表示
print(state.get_vector())

# 量子ビット1～3と固有ベクトル|ψ>=(|0>±i|1>)/2の測定
coeff = 1.0
observable3 = Observable(n_qbits)
observable2 = Observable(n_qbits)
observable1 = Observable(n_qbits)
observable_eigen = Observable(n_qbits)
observable3.add_operator(coeff, 'Z 3')
observable2.add_operator(coeff, 'Z 2')
observable1.add_operator(coeff, 'Z 1')
observable_eigen.add_operator(coeff, 'Y 0')
value3 = observable3.get_expectation_value(state)
value2 = observable2.get_expectation_value(state)
value1 = observable1.get_expectation_value(state)
value_eigen = observable_eigen.get_expectation_value(state)
print('3:', value3)
print('2:', value2)
print('1:', value1)
print('eigen state:', value_eigen)
```

# 結果
- 初期状態を $|000\rangle\otimes|\psi\rangle_+$ としたとき

固有値は $e^{-i3\pi/4}=e^{i5\pi/4}$ なので、0.101、つまり $|101\rangle$ が得られるはずです。 

```text
[ 2.77555756e-17-1.38777878e-17j  1.38777878e-17+2.77555756e-17j
  0.00000000e+00+1.96261557e-17j -1.96261557e-17+0.00000000e+00j
  2.77555756e-17+1.38777878e-17j -1.38777878e-17+2.77555756e-17j
  1.96261557e-17+5.55111512e-17j -5.55111512e-17+1.96261557e-17j
  2.77555756e-17+6.93889390e-17j -6.93889390e-17+2.77555756e-17j
  7.07106781e-01-1.96261557e-17j  1.96261557e-17+7.07106781e-01j
  2.77555756e-17-6.93889390e-17j  6.93889390e-17+2.77555756e-17j
  1.96261557e-17-5.55111512e-17j  5.55111512e-17+1.96261557e-17j]
3: -0.9999999999999989
2: 0.9999999999999989
1: -0.9999999999999989
eigen state: 0.9999999999999989
```
この測定はオブザーバブルの固有値が得られるので、$1$ は $|0\rangle$、$-1$ は $|1\rangle$ に対応します。なので、ちゃんと $|101\rangle$ が得られてることが分かります。また、初期状態もちゃんと設定されていることが分かります。

- 初期状態を $|000\rangle\otimes|\psi\rangle_-$ としたとき

固有値は $e^{i3\pi/4}$ なので、0.011、つまり $|011\rangle$ が得られるはずです。計算はプログラムの(☆)部分に"-"を付ければできます。

```text
[ 2.77555756e-17+1.38777878e-17j  1.38777878e-17-2.77555756e-17j
  1.96261557e-17+5.55111512e-17j  5.55111512e-17-1.96261557e-17j
  2.77555756e-17+6.93889390e-17j  6.93889390e-17-2.77555756e-17j
  7.07106781e-01-1.96261557e-17j -1.96261557e-17-7.07106781e-01j
  2.77555756e-17-6.93889390e-17j -6.93889390e-17-2.77555756e-17j
  1.96261557e-17-5.55111512e-17j -5.55111512e-17-1.96261557e-17j
  2.77555756e-17-1.38777878e-17j -1.38777878e-17-2.77555756e-17j
  0.00000000e+00+1.96261557e-17j  1.96261557e-17+0.00000000e+00j]
3: 0.9999999999999996
2: -0.9999999999999996
1: -0.9999999999999996
eigen state: -0.9999999999999996
```
ちゃんと、$|011\rangle$ になっています。初期状態もちゃんと $|\psi\rangle_-$ ですね。

# まとめ
今回は簡単な例でQPEを試してみました。結論としてはちゃんと理論通りに固有値の2進小数が得られました。もっと桁が欲しい場合はレジスター量子ビットを増やして同様にすれば、固有値を得られます。

さて、実際の問題において $A$ のエルミート性を仮定しても一般性を失わない状況を考えます。正定値あるいは負定値であることが分かっていないときは、固有値は正負入り乱れている状況です。今回の結果を見ると分かりますが、得られる固有値の情報は $2\pi$ で規格化した偏角 $\theta\in [0,1)$ です。このアルゴリズムだと元の $A$ の固有値の符号までは区別がつきません。今後調べる必要がありそうです。

[^1]: $P_j$ は $\lambda_j$ に関する固有空間への射影。
[^2]: 通常は $e^{-i\frac{\theta}{2}X},e^{-i\frac{\theta}{2}Y},e^{-i\frac{\theta}{2}Z}$ で定義されることが多い？(いわゆる右手系的に？歴浅のため断言できない)