QiskitでVQE実機実行を目指した・・・が

# QiskitでVQE実機実行を目指した・・・が

これまで勉強したノウハウを生かして、QiskitでVQE実機実行を目指します。

- [Qiskitで任意オブザーバブルの期待値計算をする方法](https://qiita.com/notori48/items/0bdeba72b8010e800cc8)
- [Qiskitで任意オブザーバブルの期待値を計算する方法 (measurement base)](https://qiita.com/notori48/items/c67c24d2a2fb761f4716)
- [QiskitでVQE初歩](https://qiita.com/notori48/items/4be982db033941437bf0)
- [Qiskitでトロッター分解を用いた量子シミュレーション (statevector__simulator編)](https://qiita.com/notori48/items/396d7148482e005ba507)


実機では以下がポイントです。
- 期待値はZ基底での測定から計算するしかないので、座標変換が必要　→　qiskit.aquaで自動化

なので、VQE実機実行できると思ったのですが、**うまくいきませんでした。**
- そもそもVQE実機実行＝量子回路の繰り返し実行　に時間がかかりすぎる。**1時間待っても終わりませんでした・・・**
- 出来たとしても、測定からの期待値計算はモンテカルロ法であるために、期待値計算の精度があまり高くはない　→　VQE解も精度が良くない

二番目の点が重要なので、見ていきます。

# 期待値計算

```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
from qiskit.aqua.operators.expectations import AerPauliExpectation, MatrixExpectation
from qiskit.ignis.mitigation.measurement import CompleteMeasFitter
from qiskit.providers.aer.noise import NoiseModel
from qiskit.providers.aer import QasmSimulator
import numpy as np


# Hamiltonian
from qiskit.aqua.operators import X, Y, Z, I
from qiskit.aqua.operators.state_fns import CircuitStateFn
from qiskit import Aer
from qiskit.aqua import QuantumInstance
from qiskit.aqua.operators import PauliExpectation, CircuitSampler, StateFn
obs =  (-3.8505/2 * I^I) + (-0.2288/2 * I^X) + (-1.0466/2 * I^Z) + (-0.2288/2 * X^I) + (+0.2613/2 * X^X) + (+0.2288/2 * X^Z) + (-1.0466/2 * Z^I) +  (+0.2288/2 * Z^X) + (0.2356/2 * Z^Z)

# Hamiltonian for state vector simulator
from qiskit.quantum_info.operators import Operator, Pauli
from qiskit.quantum_info import Statevector
pauli0 = Pauli(label='II')
pauli1 = Pauli(label='IX')
pauli2 = Pauli(label='IZ')
pauli3 = Pauli(label='XI')
pauli4 = Pauli(label='XX')
pauli5 = Pauli(label='XZ')
pauli6 = Pauli(label='ZI')
pauli7 = Pauli(label='ZX')
pauli8 = Pauli(label='ZZ')

obs_op = (-3.8505/2*Operator(pauli0)-0.2288/2*Operator(pauli1)-1.0466/2*Operator(pauli2)-0.2288/2*Operator(pauli3)+0.2613/2*Operator(pauli4)+0.2288/2*Operator(pauli5)-1.0466/2*Operator(pauli6)+0.2288/2*Operator(pauli7)+0.2356/2*Operator(pauli8))

# 2-qubitの回路
n = 2
# Variational circuit のパラメータ数
n_param = 6;


Method_expectation_calculation = 4 
# 0: Exact, 1: measurement based sampler, 2: Matrix based sampler, 3:State vector , 4: State vec. with Noisy model


def cost(phi):
  
  #Variational circuit
  circuit = QuantumCircuit(n)
  circuit.rx(-1*phi[0],0)
  circuit.rz(-1*phi[1],0)
  circuit.rx(-1*phi[2],1)
  circuit.rz(-1*phi[3],1)
  circuit.barrier()
  circuit.cx(0,1)
  circuit.barrier()
  circuit.rz(-1*phi[4],1)
  circuit.rx(-1*phi[5],1)
  psi = CircuitStateFn(circuit)

  coupling_map=coupling_map, measurement_error_mitigation_cls=CompleteMeasFitter, cals_matrix_refresh_period=30)

  if Method_expectation_calculation == 0:
    # easy expectation value, use for small systems only!
    Expectation = np.real(psi.adjoint().compose(obs).compose(psi).eval())

  elif Method_expectation_calculation == 1:
    backend = Aer.get_backend('qasm_simulator') 
    q_instance = QuantumInstance(backend, shots=8192)
    # define the state to sample
    measurable_expression = StateFn(obs, is_measurement=True).compose(psi) 
    # convert to expectation value
    expectation = PauliExpectation().convert(measurable_expression)  
    sampler = CircuitSampler(q_instance).convert(expectation)    # get state sampler (you can also pass the backend directly)
    Expectation = sampler.eval().real
  
  elif Method_expectation_calculation == 2:
    backend = Aer.get_backend('qasm_simulator') 
    q_instance = QuantumInstance(backend, shots=8192)
    # define the state to sample
    measurable_expression = StateFn(obs, is_measurement=True).compose(psi) 
    # Matrix based expectation
    expectation = MatrixExpectation().convert(measurable_expression)
    sampler = CircuitSampler(backend).convert(expectation)  
    Expectation = sampler.eval().real

  elif Method_expectation_calculation == 3:
    backend = Aer.get_backend('statevector_simulator')
    results = execute(circuit, backend, shots=1024).result()  
    result_statevec = Statevector(results.get_statevector())
    Expectation = np.real(result_statevec.expectation_value(obs_op))

  elif Method_expectation_calculation == 4:
    backend = provider.get_backend('ibmqx2')
    coupling_map = backend.configuration().coupling_map
    noise_model = NoiseModel.from_backend(backend.properties())
    q_instance = QuantumInstance(backend=Aer.get_backend("qasm_simulator"), shots=8192, noise_model=noise_model,  coupling_map=coupling_map, measurement_error_mitigation_cls=CompleteMeasFitter, cals_matrix_refresh_period=30)
    measurable_expression = StateFn(obs, is_measurement=True).compose(psi) 
    # convert to expectation value
    expectation = PauliExpectation().convert(measurable_expression)  
    sampler = CircuitSampler(q_instance).convert(expectation)    # get state sampler (you can also pass the backend directly)
    Expectation = sampler.eval().real

  return Expectation
```



ここで、期待値計算の手法による違いを見てみます。
期待値を適当な変分パラメータで計算します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/23c5929a-4681-f5b6-f3cc-610081d797ad.png)

 - Exact : 状態ベクトルとオブザーバブルから（線形代数によって）期待値を出しています
 - State vec : 基本的にExcatと同じです。state vector simulator から状態ベクトルを得て、state vectorクラスを作成して expectation_value を使っています。
 - Matrix sampler : Circuit sampler と Matrix を使っています。いまいち挙動がよくわかっていませんが、期待値は決定論的です。
 - Measurement base sampler : 座標変換と測定結果（Z基底）を組み合わせて期待値を計算します。実機でも動きます。
 - Measurement base sampler on Noisy model : ibmqx2の雑音を模擬したモデルを用いて、Measurement base samplerをします
 
Exact, State vec, Matrix sampler は決定論的（なんど実行しても結果が不変）で、値の精度も良さそうです。

他方、Measurement base sampler　は、測定ベース、すなわちモンテカルロ法で期待値を出すため、結果がゆらぎます。
測定回数が8192であっても、$10^{-2}$程度の精度であることがわかります。
これは、量子コンピューターは関係なくて、モンテカルロ法からの当然の帰結です。例えば、
[モンテカルロ法と円周率の近似計算](https://mathtrain.jp/montecarlo)
ここから、精度を10倍上げたければ、測定回数は100倍にする必要があります。

また、Measurement base sampler on Noisy modelの場合、$10^{-1}$の精度しか出なくなっています。
これは、Qiskitのチュートリアルとも大きく違わない結果です。
[Simulating Molecules using VQE](https://qiskit.org/textbook/ch-applications/vqe-molecules.html)

まとめますと、実機実行を考えると
- モンテカルロ法であるために、精度は$10^{-2}$
- 雑音があるので更に劣化して精度は$10^{-1}$になる？
ということで、あまり計算精度はよろしくないことが予想されます。
雑音に関しては、事前キャリブレーションを用いたerror mitigation と呼ばれる方法を使うとある程度抑え込めます。

# VQE

では期待値計算を反復して、変分パラメータを最適化してみます。

```python
import scipy.optimize
import matplotlib.pyplot as plt

cost_val = [] #コスト関数の変化を保存するための関数

#この関数がiteration ごとに呼ばれる。
def callback(phi):
    global cost_val
    cost_val.append(cost(phi))

np.random.seed(0)
init = np.random.rand(n_param)
callback(init)
res = scipy.optimize.minimize(cost, init, method='Powell', callback=callback)

plt.plot(cost_val)
plt.xlabel("iteration")
plt.ylabel("energy expectation value")
plt.show()
print(res.fun)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/9c4ef780-5145-7774-0942-ff37963c6cfd.png)

対角化による厳密解とよく一致したのは、期待値計算を厳密にしたものと、同じことですがState vectorを用いたものです。
Matrix sampler を用いた方法は、期待値計算精度は良かったはずなのに、反復すると解の精度がよろしくありません。理由不明です。。
Measurement base samplerの場合は、期待値計算精度がそもそも$10^{-2}$だったために、VQE解も$10^{-2}$の精度しか出ていません。


# 結論
VQEで精度を出すことは容易ではない。