QiskitでQAOAを実機実行してみる

# お題
DojoのQAOAの章をqiskitにし、実機実行してみます。
https://dojo.qulacs.org/ja/latest/notebooks/5.3_quantum_approximate_optimazation_algorithm.html

# QAOAパラメータ最適化(state vector simulator)

```python
#必要なモジュールのインポート
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
pauli0 = Pauli(label='ZZII')
pauli1 = Pauli(label='IZZI')
pauli2 = Pauli(label='IIZZ')
pauli3 = Pauli(label='ZIIZ')

obs_op = 0.5*(Operator(pauli0) + Operator(pauli1) + Operator(pauli2) + Operator(pauli3))
obs = obs_op.data
print(obs)

from qiskit.quantum_info import Statevector
backend = Aer.get_backend('statevector_simulator')

# 4-qubitの回路
n = 4

# circuit に U_C(gamma) を加える関数
def add_U_C(circuit, gamma): 
    for i in range(n):
      j = (i+1) % n
      circuit.cx(i, j)  
      circuit.rz(+2*gamma,j) ## qiskitでは RZ(theta)=e^{-i*theta/2*Z}      
      circuit.cx(i, j)
    return circuit

# circuit に U_X(beta) を加える関数
def add_U_X(circuit, beta):
    for i in range(n):
      circuit.rx(+2*beta,i)
    return circuit

# p=1 の |beta, gamma> を作って <beta, gamma| C(Z) |beta, gamma> を返す関数
# x = [beta, gamma]
def QAOA_output_onelayer(x):
    beta, gamma = x

    circuit = QuantumCircuit(n)
    ## 重ね合わせを作るため、アダマールゲートをかける
    for i in range(n):
        circuit.h(i)
    
    circuit.barrier()    
    ## U_C, U_Xをかける
    circuit =  add_U_C(circuit, gamma)
    circuit.barrier()
    circuit =  add_U_X(circuit, beta)

    ## |beta, gamma>を作る
    results = execute(circuit, backend, shots=1024).result()
    result_statevec = Statevector(results.get_statevector())
    expectation = result_statevec.expectation_value(obs_op)
    return np.real(expectation)

from scipy.optimize import minimize
## 初期値
x0 = np.array( [0.1, 0.1 ])

## scipy.minimize を用いて最小化
result = minimize(QAOA_output_onelayer, x0, options={'maxiter':500}, method='powell')
print(result.fun) # 最適化後の値
print(result.x) # 最適化後の(beta, gamma)

```

`-0.9999999994991852 (コスト関数)
[1.17809152 0.39269362] (最適パラメータ)`

このパラメータで回路を作って、状態ベクトルの重みから確率分布を書いてみます。

```python

# 最適なbeta, gammaを使って |beta, gamma> をつくる
beta_opt, gamma_opt = result.x

circuit = QuantumCircuit(n)
## 重ね合わせを作るため、アダマールゲートをかける
for i in range(n):
    circuit.h(i)
circuit.barrier()    
## U_C, U_Xをかける
circuit =  add_U_C(circuit, gamma_opt)
circuit.barrier()
circuit =  add_U_X(circuit, beta_opt)

## |beta, gamma>を作る
results = execute(circuit, backend, shots=1024).result()
final_state = results.get_statevector()

## z方向に観測した時の確率分布を求める. (状態ベクトルの各成分の絶対値の二乗=観測確率)
probs = np.abs(final_state**2)
print(probs)

# プロットする
import matplotlib.pyplot as plt
%matplotlib inline

## z方向に射影測定した時に得られる可能性があるビット列
z_basis = [format(i,"b").zfill(n) for i in range(probs.size)]

plt.figure(figsize=(10, 5))
plt.xlabel("states")
plt.ylabel("probability(%)")
plt.bar(z_basis, probs*100)
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/563a31e5-081a-a86f-b539-b84c7a031995.png)


Dojoと一致してます。0101と1010が解ですね。

原理的には、パラメータ最適化も実機でやりたいところですが、
iterationに時間がかかりすぎるので、この方法を取っています。


# 実機実行

```python

## "MY_API_TOKEN" に、上で取得したtokenを入れる
from qiskit import IBMQ
IBMQ.save_account('hogehoge',overwrite = True)

# least busy な バックエンドを知る
from qiskit.providers.ibmq import least_busy

# 自分のアカウント情報をloadする。（あらかじめ IBMQ.save_account を実行しておく必要がある. 複数のアカウントを使い分ける時はここで行う)
provider = IBMQ.load_account()

# 自分のアカウントで使用できるバックエンドを見る
print(provider.backends())

# Run our circuit on the least busy backend. Monitor the execution of the job in the queue
from qiskit.tools.monitor import job_monitor

circuit.measure_all() # add measurement
circuit.draw('mpl')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/bd6ef838-24d2-c90a-858f-28b566735178.png)

まず、qasm simulator でやってみます。（確率分布からサンプリングし、ヒストグラムを書きます）

```python
results = execute(circuit, backend=Aer.get_backend('qasm_simulator'), shots=8192).result() #qasm simulator
counts=results.get_counts()
plot_histogram(counts)

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/df847631-3a51-ef79-bff5-57b9f3412451.png)

状態ベクトルから得られる確率分布と（当然ですが）一致してますね。
続いて、同じ方法で実機実行します。

```python
job= execute(circuit, backend=provider.get_backend('ibmqx2'), shots=8192, optimization_level = 0) #ibmq
job_monitor(job, interval=2)
results = job.result()
counts=results.get_counts()
plot_histogram(counts)

```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/2bc2f213-7ea6-1214-ac06-115ce1093020.png)

はい、全然ダメでした。あまりにダメなので、ちゃんと動いているか心配ですが。

optimization_level=3 でも同様でした。
error mitigationでどうにかなる気もしないです・・・。

0000がたくさん出てきているということは、基底状態に落っこちてきてしまっているのでしょうか？

**追記**
Noise model + simulator の場合は、このような挙動はしませんでした。
変分パラメータを先行研究[Quantum Algorithm Implementations for Beginners](https://arxiv.org/abs/1804.03719)と同じにし、
更に同じマシン(ibmqx2)かつ同じ量子ビットアサインにしたとしても、やはり動作しない（論文が再現しない）ため、何やら異常な挙動のようです。
しかし原因がわかりません。

# 結論
QAOAを実機で動かすには工夫が必要そう。