Qiskit: VQEを実装してみた

# QiskitでVQEを実装してみた．

今回はIBMが提供している量子コンピュータ用のオープンソースフレームワークを用いてVQE(Variational Quantum Eigensolver)を実装していきます．このアルゴリズムは量子化学への応用が期待されているアルゴリズムです．理論的な解説は[Quantum Native Dojo](https://dojo.qulacs.org/ja/latest/notebooks/5.1_variational_quantum_eigensolver.html)を参照お願いします．
VQEの実装記事は複数あるとは思いますが，日本語記事の殆どはQulacsやBlueqatを用いて書かれているので今回Qiskitで実装を試みてみました．過去にQAOAやQCLの記事も書いているので是非みてみてください．
[Qiskit Aquaを使わないQAOAの実装](https://qiita.com/hotaru_46/items/c4b51e3038d260ef1fc8)
[Quantum Circuit Learningの実装](https://qiita.com/hotaru_46/items/4a7f7ddae1b0e4c7a96b)

# 実装

ライブラリのimport

```python:
from qiskit import Aer, execute
from qiskit import QuantumCircuit
from qiskit.aqua.utils import tensorproduct
from qiskit.quantum_info.analysis import average_data

from scipy.optimize import minimize
import numpy as np
```

初期化

```python:
class VQE:

    def __init__(self, n_qubits, observable, layer, backend):
        self.n_qubits = n_qubits
        self.OBS = observable
        self.Layer = layer
        self.backend = backend
```

今回はHamiltonianも入力として扱います．

回路の作成

```python:
    def make_circuit(self, Params):

        def make_U_circuit(circ, params):
            for n in range(self.n_qubits):
                param = params[3*n:3*(n+1)]
                circ.rz(param[0], n)
                circ.rx(param[1], n)
                circ.rz(param[2], n)
            return circ

        def make_Ent(circ):
            for n in range(self.n_qubits-1):
                circ.cx(n, n+1)
            return circ

        # 回路を設定
        circ = QuantumCircuit(self.n_qubits, self.n_qubits)

        for l in range(self.Layer):
            # parameter circuitの作成
            params = Params[3*self.n_qubits*l:3*self.n_qubits*(l+1)]
            # entanglementの作成
            make_U_circuit(circ, params)
            # Layerの数-1にしなければならないので
            if l != self.Layer-1:
                make_Ent(circ)

        # 測定
        circ.measure(circ.qregs[0], circ.cregs[0])

        return circ
```

VQE実行フェーズ

```python:
    def outputlayer(self, params):
        circ = self.make_circuit(params)
        counts = execute(circ, backend=self.backend, shots=8192).result().get_counts()
        return average_data(counts, self.OBS)

    def initial_params(self):
        # 初期parameterの作成
        init = [0.1 for _ in range(3 * self.Layer * self.n_qubits)]
        return np.array(init)

    def minimize(self):
        initial_params = self.initial_params()
        # 最適化の実行
        opt_params, opt_cost = classica_minimize(self.outputlayer, initial_params, options={'maxiter':500})
        circ = self.make_circuit(opt_params)
        counts = execute(circ, backend=self.backend, shots=8192).result().get_counts()
        ans = sorted(counts.items(), key=lambda x: x[1], reverse=True)
        print(ans)
        return opt_cost
```

なお，他の記事と異なりlayerというのを使用しています．これは参考文献[1]に載っているものです．軽く解説させていただきます．
$\theta=(\theta_1, \theta_2, \cdots, \theta_d)$をパラメータの集合とした時，回路は

```math:
U(\theta) = U_d(\theta_d)U_{ENT} \cdots U_1(\theta_1)U_{ENT}U_0(\theta_0)
```

と書けます．この時の$d$がlayer数となります．$d=2$とすることで他の記事と同じ回路を作成することが可能となります．

# 実行

今回はHamiltonianとして他の記事でも扱っている式を用います．
[Supplementary Information](https://static-content.springer.com/esm/art%3A10.1038%2Fncomms5213/MediaObjects/41467_2014_BFncomms5213_MOESM1050_ESM.pdf)

```python:
def sample_hamiltonian():
    '''
    https://dojo.qulacs.org/ja/latest/notebooks/5.1_variational_quantum_eigensolver.html
    のHamiltonianを使用．
    '''
    I_mat = np.array([[1, 0], [0, 1]])
    X_mat = np.array([[0, 1], [1, 0]])
    Z_mat = np.array([[1, 0], [0, -1]])
    obs = np.zeros((4, 4))
    obs += -3.8505 * tensorproduct(I_mat, I_mat)
    obs += -0.2288 * tensorproduct(I_mat, X_mat)
    obs += -1.0466 * tensorproduct(I_mat, Z_mat)
    obs += 0.2613 * tensorproduct(X_mat, X_mat)
    obs += 0.2288 * tensorproduct(X_mat, Z_mat)
    obs += -1.0466 * tensorproduct(Z_mat, I_mat)
    obs += 0.2288 * tensorproduct(Z_mat, X_mat)
    obs += 0.2356 * tensorproduct(Z_mat, Z_mat)
    return obs / 2
```

また，古典最適化として私が書いている他の記事でも用いている関数を使います．

```python:
def classica_minimize(cost_func, initial_params, options, method='powell'):
    print('classical minimize is starting now... ')
    result = minimize(cost_func, initial_params, options=options, method=method)
    print('opt_cost: {}'.format(result.fun))
    print('opt_params: {}'.format(result.x))
    return result.x, result.fun
```

実行です.

```python:
if __name__ == '__main__':
    backend = Aer.get_backend('qasm_simulator')
    obs = sample_hamiltonian()
    vqe = VQE(2, obs, 2, backend)
    vqe.minimize()
```

# 実行結果

```
-2.85405 # vqeの結果
-2.8626207640766816　# 最小固有値
```

他の記事と比較すると精度は落ちますけど一応成功ですかね...?

# まとめ

Qiskitを用いてVQEの実装を行ってみました．
正直私自身はVQEよりも古典的問題に適しているQAOAの方が好きなので完全に実装を目的として記事を書かせていただきました．
理論等に関してはご自身でお願いします...

# 参考文献

[1] [Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets](https://arxiv.org/abs/1704.05018)
[2] [Quantum Native Dojo](https://dojo.qulacs.org/ja/latest/notebooks/5.1_variational_quantum_eigensolver.html)