【Quantum Native Dojo】QulacsでVQEを実装してみた【5-1】

[Quantum Native Dojoの5-1. Variational Quantum Eigensolver（VQE）アルゴリズム](https://dojo.qulacs.org/ja/latest/notebooks/5.1_variational_quantum_eigensolver.html)はpythonのnumpyで実装されていますが、Qulacsを使って同じVQEを実装してみました。
(実装が間違っていたら教えてください。)  

まずコードの解説をします。
次にQuantum Native Dojoの5-1で使われているものと同じ変分回路で計算した結果を報告します。
また変分回路でCNOTゲートを作用させなかった場合やパラメータの数を増やした場合の計算結果も報告します。
最後にCNOTゲートがなぜ必要なのか考察してみました。

# コードの解説
コード全体は以下のgithubリポジトリにあります。 

https://github.com/Noriaki416/Qiita/blob/master/5.1_VQE.ipynb

以下で詳しく説明していきます。

## モジュールのインポート

```python
from qulacs import Observable
from qulacs import QuantumState, QuantumCircuit
from qulacs.gate import RX, RZ, CNOT
from scipy.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
```

## 変分回路の作成

```python
def make_ansatz_circuit(phi, n_bits):
    circuit = QuantumCircuit(n_bits)
    circuit.add_gate(RX(0, phi[0]))
    circuit.add_gate(RZ(0, phi[1]))
    circuit.add_gate(RX(1, phi[2]))
    circuit.add_gate(RZ(1, phi[3]))
    circuit.add_gate(CNOT(0, 1))
    circuit.add_gate(RX(1, phi[4]))
    circuit.add_gate(RZ(1, phi[5]))
    return circuit
```

`phi`に最適化するパラメータのリスト、`n_bits`に量子ビット数を入力します。  

## ハミルトニアンの作成
第２量子化を行いパウリ演算しで表現した水素分子のハミルトニアンを準備します。

```python
hamiltonian = [(-3.8505/2, ""),  \
               (-0.2288/2, "X 1"), \
               (-1.0466/2, "Z 1"),\
               (-0.2288/2, "X 0"), \
               (0.2613/2,  "X 0 X 1"), \
               (0.2288/2,  "X 0 Z 1"), \
               (-1.0466/2,  "Z 0"), \
               (0.2288/2,  "Z 0 X 1"), \
               (0.2356/2,  "Z 0 Z 1")]
```

Qulacsの`Observable`に上記の演算子を追加していく。

```python
n_bits = 2
observable = Observable(n_bits)
for (coef, Pauli_string) in  hamiltonian :
    observable.add_operator(coef, Pauli_string)
```

## コスト関数の作成
初期状態に`make_ansatz_circuit`で作成した変分回路を作用させて`observable`の期待値を計算します。

```python
def cost(phi):
    state = QuantumState(n_bits)
    state.set_zero_state()
    circuit = make_ansatz_circuit(phi, n_bits)
    circuit.update_quantum_state(state)
    return observable.get_expectation_value(state)
```

最適化を行う過程で上で計算した期待値を`cost_val`に順次保存していくための関数です。

```python
def callback(phi):
    global cost_val
    cost_val.append(cost(phi))
```

## 最適化
シミュレーターなので意味があるかどうかは分かりませんが同じ計算を100回を行いその平均値を最後に出力しています。  
グラフは最後のiterationの結果だけを表示させています。

```python
n_params = 6
np.random.seed(1111)
n_iter = 100
energy_list = []
for _ in range(n_iter):  
    cost_val = []
    phi = np.random.random(n_params)
    callback(phi)
    res = minimize(cost, phi, method='Powell', callback=callback)
    energy_list.append(cost(res.x))

plt.plot(cost_val)
plt.xlabel("iteration")
plt.ylabel("energy expectation value")
plt.show()

print("Exact Energy: ", -2.8626207640766816)
print("VQE Energy: {}".format(np.mean(energy_list)))
```

# 計算結果
出力結果は以下のようになります。

```output
Exact Energy:  -2.8626207640766816
VQE Energy:    -2.8624166416549484
```

厳密対角化を行った結果と小数第３位まで一致しています。  

## 計算結果にノイズが入った場合
`cost`と`callback`を以下のように変更します。
`cost_noise`では最後にノイズが追加されています。

```python
def cost_noise(phi):
    state = QuantumState(n_bits)
    state.set_zero_state()
    circuit = make_ansatz_circuit(phi, n_bits)
    circuit.update_quantum_state(state)
    return observable.get_expectation_value(state) + np.random.normal(0, 0.01)

def callback_noise(phi):
    global cost_val
    cost_val.append(cost_noise(phi))
```

最適化は上記で紹介したコードの`cost`の部分を`cost_noise`、`callback`の部分を`callback_noise`に変更すれば実行できます。
結果は

```output
Exact Energy:  -2.8626207640766816
VQE Energy:    -2.850086422671339
```

小数第２位以降で一致しないようになりました。  

## パラメータを増やした場合
### ノイズなし
`make_ansatz_circuit`を以下のように変更します。

```python
def make_ansatz_circuit(phi, n_bits):
    circuit = QuantumCircuit(n_bits)
    circuit.add_gate(RX(0, phi[0]))
    circuit.add_gate(RZ(0, phi[1]))
    circuit.add_gate(RX(1, phi[2]))
    circuit.add_gate(RZ(1, phi[3]))
    circuit.add_gate(CNOT(0, 1))
    circuit.add_gate(RX(1, phi[4]))
    circuit.add_gate(RZ(1, phi[5]))
    circuit.add_gate(RX(0, phi[6]))
    circuit.add_gate(RZ(0, phi[7]))
    circuit.add_gate(CNOT(0, 1))
    return circuit
```

最適化は上記で紹介したコードの`n_params = 6`の部分を`n_params = 8`に変更すれば実行できます。
結果は

```output
Exact Energy:  -2.8626207640766816
VQE Energy:    -2.8626194809589203
```

小数第4位まで一致するようになりました！ 
パラメータの数が6個の変分回路では`-2.8624166416549484`だったので精度が向上しています。

### ノイズあり
上と同じ変更を施したコードでノイズが入った時に使ったコードを`n_params = 8`に変更すれば実行できます。
結果は

```output
Exact Energy:  -2.8626207640766816
VQE Energy:    -2.8457532192348656
```

パラメータ数が6個でノイズが入った場合の結果が`-2.850086422671339`だったので逆に精度が下がりました。

## CNOTゲートなし、ノイズなし
最後にCNOTゲートなし、ノイズなしで計算を行ってみます。
`make_ansatz_circuit`からCNOTゲートを削除します。

```python
def make_ansatz_circuit(phi, n_bits):
    circuit = QuantumCircuit(n_bits)
    circuit.add_gate(RX(0, phi[0]))
    circuit.add_gate(RZ(0, phi[1]))
    circuit.add_gate(RX(1, phi[2]))
    circuit.add_gate(RZ(1, phi[3]))
    circuit.add_gate(RX(1, phi[4]))
    circuit.add_gate(RZ(1, phi[5]))
    return circuit
```

結果は

```output
Exact Energy:  -2.8626207640766816
VQE Energy:    -2.85404959554127
```

CNOTゲートがある場合は`-2.8624166416549484`だったので精度が下がりました。
CNOTゲートを作用させた方が良いことが分かりました。

# 考察
CNOTゲートを作用させる場合とさせない場合でどういった違いがあるのかを考察したいと思います。
今回は2量子ビットの状態なので具体的に手計算を行いながら分析したいと思います。

$$
\def\bm#1{\boldsymbol{#1}}
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\brakett#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
\def\braket#1#2#3{\bra{#1}#2\ket{#3}}
$$

初期状態を$\ket{\psi_{init}}=\ket{00}$として、`make_ansatz_circuit`の演算子を順番に作用させていきます。
(左から0ビット、1ビットと数えることにします。)

```python
def make_ansatz_circuit(phi, n_bits):
    circuit = QuantumCircuit(n_bits)
    circuit.add_gate(RX(0, phi[0]))
    circuit.add_gate(RZ(0, phi[1]))
    circuit.add_gate(RX(1, phi[2]))
    circuit.add_gate(RZ(1, phi[3]))
    circuit.add_gate(CNOT(0, 1))
    circuit.add_gate(RX(1, phi[4]))
    circuit.add_gate(RZ(1, phi[5]))
    return circuit
```

まず回転ゲートは以下のように定義されています。

```math
\begin{align}
R_{X_i}(\theta) &= e^{-\frac{i\theta}{2} X_i} = \cos{(\frac{\theta}{2})} I_i - i \sin{(\frac{\theta}{2})} X_i \\
R_{Z_i}(\theta) &= e^{-\frac{i\theta}{2} Z_i} = \cos{(\frac{\theta}{2})} I_i - i \sin{(\frac{\theta}{2})} Z_i
\end{align}
```

$a(\theta) = \cos{(\frac{\theta}{2})}, b(\theta) = - i \sin{(\frac{\theta}{2})}$とおくと、  

```math
\begin{align}
R_{X_i}(\theta) &= a(\theta) I_i + b(\theta) X_i \\
R_{Z_i}(\theta) &= a(\theta) I_i + b(\theta) Z_i
\end{align}
```

と表すことが出来ます。  
1ビット状態$\ket{0}, \ket{1}$にRX, RZを順番に作用させて得られる状態をあらかじめ計算しておきます。
まずは$\ket{0}$に対して

```math
\begin{align}
R_X(\theta_0)\ket{0} &= a(\theta_0)\ket{0}  + b(\theta_0)\ket{1} \\
R_Z(\theta_1)R_X(\theta_0)\ket{0} &= a(\theta_1)\biggl( a(\theta_0)\ket{0}  + b(\theta_0)\ket{1} \biggr) + b(\theta_1) \biggl(a(\theta_0)\ket{0}  - b(\theta_0)\ket{1} \biggr) \\
&= \biggl(a(\theta_0)a(\theta_1) + a(\theta_0)b(\theta_1) \biggr)\ket{0} + \biggl(b(\theta_0)a(\theta_1) - b(\theta_0)b(\theta_1) \biggr)\ket{1} \\
& \equiv A(\theta_0, \theta_1)\ket{0} + B(\theta_0, \theta_1)\ket{1}
\end{align}
```

となります。$\ket{1}$に対して、

```math
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\begin{align}
R_X(\theta_0)\ket{1} &= a(\theta_0)\ket{1}  + b(\theta_0)\ket{0} \\
R_Z(\theta_1)R_X(\theta_0)\ket{1} &= a(\theta_1)\biggl( a(\theta_0)\ket{1}  + b(\theta_0)\ket{0} \biggr) + b(\theta_1) \biggl(-a(\theta_0)\ket{1}  + b(\theta_0)\ket{0} \biggr) \\
&= \biggl(b(\theta_0)a(\theta_1) + b(\theta_0)b(\theta_1) \biggr)\ket{0} + \biggl(a(\theta_0)a(\theta_1) - a(\theta_0)b(\theta_1) \biggr)\ket{1} \\
& \equiv C(\theta_0, \theta_1)\ket{0} + D(\theta_0, \theta_1)\ket{1}
\end{align}
``` 

となります。
簡単のために係数$A,B,C,D$は$A(\theta_i, \theta_j)=A_{ij}$とそれぞれ略記することにします。

`circuit.add_gate(~))`で演算子が作用されるたびに状態を$\ket{\psi_{i-1}}\to\ket{\psi_i}$と更新して表記していきます。
最初は$\ket{\psi_0}=R_{X_0}(\theta_0)\ket{\psi_{init}}$と表すことにします。
CNOTゲートの前までは以下のように状態が更新されていきます。

```math
\begin{align}
\ket{\psi_0} &= R_{X_0}(\theta_0)\ket{\psi_{init}} \\
\ket{\psi_1} &= R_{Z_0}(\theta_1)\ket{\psi_{0}} \\
\ket{\psi_2} &= R_{X_1}(\theta_2)\ket{\psi_{1}} \\
\ket{\psi_3} &= R_{Z_1}(\theta_3)\ket{\psi_{2}}
\end{align}
```

CNOTゲートの手前まではそれぞれにビットに独立して演算子が作用するので

```math
\begin{align}
\ket{\psi_3} &= R_{Z_0}(\theta_1)R_{X_0}(\theta_0)\ket{0} \otimes R_{Z_1}(\theta_3)R_{X_1}(\theta_2)\ket{0} \\
&= \biggl(A_{01}\ket{0} + B_{01}\ket{1} \biggr) \otimes \biggl(C_{23}\ket{0} + D_{23}\ket{1} \biggr) \\
&= A_{01}C_{23}\ket{00} + A_{01}D_{23}\ket{01} +  B_{01}C_{23}\ket{10} + B_{01}D_{23}\ket{11}
\end{align}
```

となります。
CNOTゲートを作用させると、

```math
\begin{align}
\ket{\psi_4} &= CNOT(0, 1)\ket{\psi_{3}} \\
&= A_{01}C_{23}\ket{00} + A_{01}D_{23}\ket{01} +  B_{01}C_{23}\ket{11} + B_{01}D_{23}\ket{10}\\
&= \biggl(A_{01}C_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \ket{0} + \biggl(A_{01}D_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \ket{1}
\end{align}
```

になります。
最後に1ビット目に$R_X, R_Z$を続けて作用させると

```math
\begin{align}
\ket{\psi_5} &= R_{X_1}(\theta_4)\ket{\psi_{4}} \\
\ket{\psi_6} &= R_{Z_1}(\theta_5)\ket{\psi_{5}}
\end{align}
```

なので最終的な状態は、

```math
\begin{align}
\ket{\psi_6} &= \biggl(A_{01}C_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \biggl(A_{45}\ket{0} + B_{45}\ket{1} \biggr) \\
& \;\;\; + \biggl(A_{01}D_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \biggl(C_{45}\ket{0} + D_{45}\ket{1} \biggr)\\
&= \biggl(A_{01}C_{23}A_{45} +A_{01}D_{23}C_{45}\biggr)\ket{00} + \biggl(A_{01}C_{23}B_{45} +A_{01}D_{23}D_{45}\biggr)\ket{01} \\
& \;\;\; + \biggl(B_{01}D_{23}A_{45} +B_{01}C_{23}C_{45}\biggr)\ket{10} + \biggl(B_{01}D_{23}B_{45} +B_{01}C_{23}D_{45}\biggr)\ket{11} \\
\end{align}
```

となります。
もしCNOTゲートを作用させなければ、


```math
\begin{align}
\ket{\psi^{\prime}_4} &= \ket{\psi_{3}} \\
&= A_{01}C_{23}\ket{00} + A_{01}D_{23}\ket{01} +  B_{01}C_{23}\ket{10} + B_{01}D_{23}\ket{11}\\
&= \biggl(A_{01}C_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \ket{0} + \biggl(A_{01}D_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \ket{1}
\end{align}
```


より、


```math
\begin{align}
\ket{\psi^{\prime}_6} &= \biggl(A_{01}C_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \biggl(A_{45}\ket{0} + B_{45}\ket{1} \biggr) \\
& \;\;\; + \biggl(A_{01}D_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \biggl(C_{45}\ket{0} + D_{45}\ket{1} \biggr)\\
&= \biggl(A_{01}C_{23}A_{45} +A_{01}D_{23}C_{45}\biggr)\ket{00} + \biggl(A_{01}C_{23}B_{45} +A_{01}D_{23}D_{45}\biggr)\ket{01} \\
& \;\;\; + \biggl(B_{01}C_{23}A_{45} +B_{01}D_{23}C_{45}\biggr)\ket{10} + \biggl(B_{01}C_{23}B_{45} +B_{01}D_{23}D_{45}\biggr)\ket{11} \\
\end{align}
```

$\ket{\psi_6}, \ket{\psi^{\prime}_6}$を比較すると$\ket{10}, \ket{11}$の真ん中に挟まれている係数($\theta\_2, \theta\_3$に関する係数)が異なっています。
CNOTゲートを作用させることで、$\theta\_2, \theta\_3$に関する係数が、$\ket{00}, \ket{01}$はそのままで、$\ket{10}, \ket{11}$だけ反転していることが分かります。
これは単純に回転ゲートを作用させていくだけでは得られない状態でCNOTゲートを作用させることでより複雑な(表現力の高い？)状態を生成していることが分かります。
このような状態を利用して期待値を計算することでより精度の高いエネルギー値が得られているのだと思います。

変分回路を作る統一的な指針みたいなものはあるのでしょうか？
ディープラーニングのニューラルネットワークと同じように層の数や活性化関数の選び方がある程度経験的に基づくのと同様に変分回路も試行錯誤しながら作成していくのでしょうか。