Qiskitでトロッター分解を用いた量子シミュレーション (statevector__simulator編)

# statevector__simulator
Qiskitでトロッター分解を用いた量子シミュレーションの記事を以前に掲載しておりますが、
qasm_simulatorを用いたものでした。
https://qiita.com/notori48/items/dbac7bb0226bcfbfbf40

qasm_simulatorは実機と同様に状態作成後にサンプリングを行います。
対して、statevector__simulatorは状態作成をして、その状態を直接配列として取り出します。
（もちろん、実機ではこのようなことは通常出来ません。状態を推定する別のアルゴリズム、例えば量子状態トモグラフィなどが必要となります。この場合も推定誤差や雑音の影響が問題になります。）

今回はstatevector__simulatorを使います。

# スクリプト

```python
pip install qiskit pylatexenc
#必要なモジュールのインポート
from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute, Aer
from qiskit.qasm import pi
from qiskit.tools.visualization import plot_histogram, circuit_drawer
import numpy as np
from qiskit.quantum_info.operators import Operator, Pauli

# Observable 作成
pauli0 = Pauli(label='ZIIIII')
pauli1 = Pauli(label='IZIIII')
pauli2 = Pauli(label='IIZIII')
pauli3 = Pauli(label='IIIZII')
pauli4 = Pauli(label='IIIIZI')
pauli5 = Pauli(label='IIIIIZ')

obs_op = (1/6)*(Operator(pauli0)+Operator(pauli1)+Operator(pauli2)+Operator(pauli3)+Operator(pauli4)+Operator(pauli5))
obs = obs_op.data
print(obs)
```

```python
# ダイナミクスをシミュレーションする時間
t = 3.0
## 横磁場の強さ
h = 3.
# トロッター分解の分割数
M = 100
# 時間の刻み幅
delta = t/M

## 時間と磁化を記録するリスト
x = [i*delta for i in range(M+1)]
y = []

import matplotlib.pyplot as plt

nqubits = 6

from qiskit.aqua.utils import tensorproduct
from qiskit.quantum_info.analysis import average_data
from qiskit.quantum_info import Statevector
backend = Aer.get_backend('statevector_simulator')

for p in range(M+1):
  circuit_trotter_Ising = QuantumCircuit(nqubits) 
  for k in range(p):
    for i in range(nqubits):
      circuit_trotter_Ising.cx(i,(i+1)%(nqubits))
      circuit_trotter_Ising.rz(2*delta,(i+1)%nqubits) ## RZ(a)=exp(-i*a/2*Z)に注意. rz(rotation_angle, qbit_label)
      circuit_trotter_Ising.cx(i,(i+1)%(nqubits))
      circuit_trotter_Ising.rx(2*delta*h,i) ## RX(a)=exp(-i*a/2*X)に注意
  #circuit_trotter_Ising.measure_all()           
  results = execute(circuit_trotter_Ising, backend, shots=1024).result()
  result_statevec = Statevector(results.get_statevector())
  expectation = result_statevec.expectation_value(obs_op)
  print(expectation)
  y.append(expectation)  
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/4820111a-9c1c-77ac-8fdf-def48cc53b3f.png)

imaginary part が出てきていますが、シミュレーターの計算誤差に起因するものです。
np.real() で排除できます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/d4dce52a-9860-8ed0-1721-eca67e5ae641.png)



ちゃんと出来ています。
https://dojo.qulacs.org/ja/latest/notebooks/4.2_trotter_decomposition.html
の結果と一致してますね

# 解説
statevector_simulatorを使うには、

```python
backend = Aer.get_backend('statevector_simulator')
results = execute(circuit_trotter_Ising, backend, shots=1024).result()
```
でよいです。shotsは特に意味がない（1でもいい）かと思います。
なお、回路の最後にはmeasureを入れないことをおすすめします。
measureがあると、回路実行後に状態が確定（射影）してしまって、それが読み出されてしまいます。

シミュレーション後は、
results.get_statevector()
で状態ベクトルを読み出します。ただ、今回はここからオブザーバブルの期待値を計算したいです。
残念ながら直接はできないので、qiskit.quantum_infoにあるクラス Statevector を作ります。
そして、Statevector のメソッド（機能）としてある .expectation_value を使います。
https://qiskit.org/documentation/locale/ja_JP/stubs/qiskit.quantum_info.Statevector.html

```python
  result_statevec = Statevector(results.get_statevector())
  expectation = result_statevec.expectation_value(obs_op)
```


引数はオブザーバブルと、オブザーバブルが作用する量子ビットのようです。
今回、オブザーバブルは、例えば $ZIIIII$ のように6量子ビットフルの表示をしてますので、量子ビットの指定は必要ないです。もし$Z$と書いた場合は、.expectation_value(obs_Z, 0) のように書くのではないでしょうか。

また、オブザーバブルとしては Operator のクラスをそのまま渡すか、Operator.data で配列に持ってきてから渡すかです。
このあたり、qiskitの関数によってまちまちです。
Operator(U) とすると、行列UをOperatorに変換できます。

# 結論
statevector_simulator でも計算できた。
確かにDojoと一致することがわかった。
なお、他のバックエンドとしてunitary_simulatorというのものあるようです。
これは回路全体を記述するユニタリ行列を取るsimulatorのようです。