###### tags: `量子計算`
# 量子期待値推定 Python

参考文献: [Qunasys, Quantum Native Dojo 5章-1](https://dojo.qulacs.org/ja/latest/notebooks/5.1_variational_quantum_eigensolver.html#変分法)

## 量子期待値推定の概要
古典コンピュータを用いた変分原理に対して，量子コンピュータで準備しやすい量子状態を使って，エネルギー期待値を求めるアルゴリズム．

> 変分原理:
\begin{eqnarray}
\left<\psi\left|\hat H\right|\psi\right> \geq E
\end{eqnarray}
が成り立つとした原理．

## Pythonによる量子期待値推定
以下では，HeHの**基底エネルギー**期待値を計算するアルゴリズムをかく．使っている量子ビットは2ビットである．量子ゲートや量子状態の定義は，Qulacsは使わずに，Numpyだけを使う．最適化はSciPyを使っている．

> Q: 2ビットなのは，HeとHの２つだから？
> A: 違う．ビットの数はフェルミオンが入ることができるサイト(LCM)の個数のことである．今回は基底エネルギーの期待値を求めていて，LCMは，基底の準位に入っているスピンアップとダウンの軌道の2つである．

> ビット数が少ない時は，Qulacsを使わず，Numpyだけでもいいかもしれない．

### モジュールのインポートと行列の定義

import numpy as np
nqubits = 2
#パウリ演算子を準備する。
pI = np.array([[1+0.0j,0+0.0j],[0+0.0j,1+0.0j]])
pX = np.array([[0+0.0j,1+0.0j],[1+0.0j,0+0.0j]])
pZ = np.array([[1+0.0j,0+0.0j],[0+0.0j,-1+0.0j]])
pY = np.array([[0+0.0j,-1.0j],[0.0+1.0j,0.0+0.0j]])
pHad = (pX+pZ)/np.sqrt(2)
pP0 = (pI+pZ)/2
pP1 = (pI-pZ)/2

\begin{eqnarray}
\mathrm{pI} =&& \left(\begin{array}{cc}1&0\\0&1\end{array}\right), ~~
\mathrm{pX} = \left(\begin{array}{cc}0&1\\1&0\end{array}\right),\\
\mathrm{pY} =&& \left(\begin{array}{cc}0&-i\\i&0\end{array}\right), ~~
\mathrm{pZ} = \left(\begin{array}{cc}1&0\\0&-1\end{array}\right),\\
\mathrm{pHad} =&& \frac{1}{\sqrt{2}}(X+Z) = \left(\begin{array}{cc}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\end{array}\right),\\
\mathrm{pP0} =&& \left|0\right>\left<0\right| = \frac{1}{2}(I+Z) = \left(\begin{array}{cc}1&0\\0&0\end{array}\right),\\
\mathrm{pP1} =&& \left|1\right>\left<1\right| = \frac{1}{2}(I-Z) = \left(\begin{array}{cc}0&0\\0&1\end{array}\right).\\
\end{eqnarray}

### １量子ゲートの定義
１量子ゲート(1つの量子ビットに作用するゲート)を定義する．

#パウリ演算子を1量子ゲートに変換する。
X=[1]*(nqubits)
Y=[1]*(nqubits)
Z=[1]*(nqubits)
H=[1]*(nqubits)
P0=[1]*(nqubits)
P1=[1]*(nqubits)

for i in range(nqubits):
    for j in range(nqubits):
	## np.kronはクロネッカー積を表す．
        if(i != j):
            X[i] = np.kron(pI,X[i])
            Y[i] = np.kron(pI,Y[i])
            Z[i] = np.kron(pI,Z[i])
            H[i] = np.kron(pI,H[i])
            P0[i] = np.kron(pI,P0[i])
            P1[i] = np.kron(pI,P1[i])
        else:
            X[i] = np.kron(pX,X[i])
            Y[i] = np.kron(pY,Y[i])
            Z[i] = np.kron(pZ,Z[i])
            H[i] = np.kron(pHad,H[i])
            P0[i] = np.kron(pP0,P0[i])
            P1[i] = np.kron(pP1,P1[i])
Ide = np.eye(2**nqubits) ## 単位行列を指定する．

配列$Z$を例にfor文でやっていることを追ってみる．まず，
$(i,j)=(0,0)$では，
\begin{eqnarray}
Z[0] = \sigma_Z \otimes 1 = \sigma_Z\\
\therefore Z=[\sigma_Z,1]
\end{eqnarray}
$(i,j)=(0,1)$では，
\begin{eqnarray}
Z[0] = I \otimes Z[0] = I \otimes \sigma_Z\\
\therefore Z=[I \otimes \sigma_Z,1]
\end{eqnarray}
$(i,j)=(1,0)$では，
\begin{eqnarray}
Z[1] = I \otimes 1 = I\\
\therefore Z=[I\otimes\sigma_Z,I]
\end{eqnarray}
$(i,j)=(1,1)$では，
\begin{eqnarray}
Z[1] = \sigma_Z \otimes Z[1] = \sigma_Z \otimes I\\
\therefore Z=[I \otimes \sigma_Z,\sigma_Z \otimes I]
\end{eqnarray}

\begin{eqnarray}
Z = [I \otimes \sigma_Z,\sigma_Z \otimes I]
\end{eqnarray}
$Z[i]$は，$i$ビット目に$\sigma_Z$を作用させるゲートに対応する．他のゲートも同様にして生成する．

> 以降$\sigma_Z$を$Z$と表記することがある．

### ２量子ゲートを準備する

#2量子ゲートを準備する。
CZ = [[0 for i in range(nqubits)] for j in range(nqubits)]
CX = [[0 for i in range(nqubits)] for j in range(nqubits)]

for i in range(nqubits):
    for j in range(nqubits):
        CZ[i][j]= (P0[i]+np.dot(P1[i],Z[j]))
        CX[i][j]= (P0[i]+np.dot(P1[i],X[j]))

\begin{eqnarray}
\mathrm{CZ}
=&& [[\left|0\right>\left<0\right| \otimes I + (\left|1\right>\left<1\right|\otimes I)(Z\otimes I),\\
&&\left|0\right>\left<0\right| \otimes I + (\left|1\right>\left<1\right|\otimes I)(I\otimes Z)],\\
&& [I \otimes\left|0\right>\left<0\right| + (I \otimes\left|1\right>\left<1\right|)(Z\otimes I),\\
&& I \otimes\left|0\right>\left<0\right| + (I \otimes\left|1\right>\left<1\right|)(I\otimes Z)]
]\\
=&& [[Z \otimes I, \left|0\right>\left<0\right|\otimes I + \left|1\right>\left<1\right|\otimes Z],\\
&&[I\otimes \left|0\right>\left<0\right| + Z\otimes \left|1\right>\left<1\right|,I\otimes Z]]
\end{eqnarray}
なお，
\begin{eqnarray}
&& I\otimes \left|0\right>\left<0\right| + Z\otimes \left|1\right>\left<1\right| = \left|0\right>\left<0\right|\otimes I + \left|1\right>\left<1\right|\otimes Z
=&& \left(\begin{array}{ccc}I&\\&Z\end{array}\right)
\end{eqnarray}
である．
つまり，$\mathrm{CX,CZ}$は，対角成分が１量子ゲートで，非対角成分が制御ゲートである．

### 変分量子ゲートを準備する

#変分量子ゲート(X,Y,Zに関する回転の角度を指定できるゲート)を準備する。
from scipy.linalg import expm ## expmは行列関数を計算するモジュール．Exp Matrixの略と思われる．
def RX(target,angle):
    return expm(-0.5*angle*1.j*X[target])
def RY(target,angle):
    return expm(-0.5*angle*1.j*Y[target])
def RZ(target,angle):
    return expm(-0.5*angle*1.j*Z[target])

\begin{eqnarray}
\mathrm{RZ}(k,\theta) = \exp\left(-\frac{i}{2}\theta Z[k]\right)
\end{eqnarray}
ここで，$Z$は前述の配列であり，$Z[k]$は$k$番目の量子ビットに$Z$を作用させる１量子ビット演算子である．

### 初期状態の準備

#初期状態|0000・・・0>を準備する。
def StateZeros(nqubits):
    State = np.zeros(2**nqubits)
    State[0]=1
    return State

\begin{eqnarray}
\left|\underbrace{0..0}_ {n}\right>
\end{eqnarray}
に対応した$2^n$次元ベクトルが出力される．

### ハミルトニアンを準備する

M = (-3.8505 * Ide  - 0.2288 * X[1] - 1.0466 * Z[1] - 0.2288 * X[0] +  0.2613 * np.dot(X[0],X[1]) + \
     0.2288 *np.dot(X[0],Z[1]) - 1.0466*Z[0] + 0.2288* np.dot(Z[0],X[1]) +  0.2356 * np.dot(Z[0],Z[1]) )/2

論文のSupplement Documentから読み取ることができる．電子同士や電子と原子核相互作用の大きさから第二量子化を用いて計算したものである．(と思われる)

### 変分量子回路を準備する
論文をもとに変分量子回路を準備する

n_param = 6
def TwoQubitPQC(phi):
    state = StateZeros(2)
    state = np.dot(RX(0,phi[0]),state)
    state = np.dot(RZ(0,phi[1]),state)
    state = np.dot(RX(1,phi[2]),state)
    state = np.dot(RZ(1,phi[3]),state)
    state = np.dot(CX[1][0],state)
    state = np.dot(RZ(1,phi[4]),state)
    state = np.dot(RX(1,phi[5]),state)
    return state

\begin{eqnarray}
&& \mathrm{TwoQubitPQC}(\vec \phi)\\
=&& \mathrm{RX}(1,\phi_5)\cdot\mathrm{RZ}(1,\phi_4)\cdot\mathrm{CNOT} \cdot \mathrm{RZ}(1,\phi_3) \cdot \mathrm{RX}(1,\phi_2) \cdot \mathrm{RZ}(0,\phi_1) \cdot \mathrm{RX}(0,\phi_0)\left|00\right>
\end{eqnarray}

>なぜこれで変分量子回路が作れるか？~~状態$\left|\psi\right>$の回転角$\phi$を決定するのが目的でそのような回転角を再現する回路が変分量子回路ということ？だとすると任意の回転角が再現される必要があり，それを保証する数学はどうなっているんだろう？(オイラー角とか極座標を意識？そもそも状態の回転角って何？) - 少なくとも論文には，任意の2Qbitの状態を下の計算結果のように6この回転角を用いることで表すことができると書いてあった．~~
>任意の2ビット系を回転角が6こにより表現することができる．これについては[n Qubitの状態の回転角](https://hackmd.io/RMPndG47R_mSxaBXcuNtCA?view)を参照のこと．
>
>Qunasysのページに論文で紹介された変分量子回路をコードにしたって書いてあったけど，どうやって読み取った？ - 下の変形や角の定義を変えながら，量子回路を上手に構築した．

この量子回路を愚直に計算すると，
\begin{eqnarray}
\left|\psi\right> = &&\cos\frac{\phi_0}{2}e^{i\frac{\phi_1+\phi_3+\phi_4}{2}}\left|0\right> \left(\cos\frac{\phi_2+\phi_5}{2}\left|0\right> + i\sin\frac{\phi_2+\phi_5}{2} e^{-i(\phi_3+\phi_4)}\left|1\right> \right)\\
&&+ i\sin \frac{\phi_0}{2}e^{i \frac{\phi_1+\phi_3-\phi_4}{2}}\left|1\right>\left(\sin\frac{\phi_2+\phi_5}{2}\left|0\right> + i\cos\frac{\phi_2+\phi_5}{2} e^{-i(\phi_3-\phi_4)}\left|1\right> \right)\\
\dot= && \cos\frac{\phi_0}{2} \left|0\right> \left(\cos\frac{\phi_2+\phi_5}{2}\left|0\right> + i\sin\frac{\phi_2+\phi_5}{2} e^{-i(\phi_3+\phi_4)}\left|1\right> \right)\\
&&+ i\sin \frac{\phi_0}{2}e^{-i\phi_4} \left|1\right>\left(\cos\frac{\pi-(\phi_2+\phi_5)}{2}\left|0\right> + i\sin\frac{\pi-(\phi_2+\phi_5)}{2} e^{-i(\phi_3-\phi_4)}\left|1\right> \right)
\end{eqnarray}
となる．ここで，
\begin{eqnarray}
\theta_0 =&& \phi_0 \\
\theta_1 =&& \frac{\phi_2+\phi_5}{2} \\
\theta_2 =&& \frac{\pi-(\phi_2+\phi_5)}{2}\\
\omega_0 =&& \phi_4 \\
\omega_1 =&& \phi_3+\phi_4 \\
\omega_2 =&& \phi_3-\phi_4 \\
\end{eqnarray}
とすれば，論文の[サプリメントドキュメント](https://static-content.springer.com/esm/art%3A10.1038%2Fncomms5213/MediaObjects/41467_2014_BFncomms5213_MOESM1050_ESM.pdf)の式(7)(8)と一致する．ただし，
\begin{eqnarray}
(0 \leq \theta \leq \pi, 0 \leq \omega <2\pi)
\end{eqnarray}
である．

### 量子状態のエネルギー期待値を計算する

def ExpectVal(Operator,State):
    BraState = np.conjugate(State.T) #列ベクトルを行ベクトルへ変換
    tmp = np.dot(BraState,np.dot(Operator,State)) #行列を列ベクトルと行ベクトルではさむ
    return np.real(tmp) #要素の実部を取り出す

期待値を次式のように計算する.
\begin{eqnarray}
O = \mathrm{ExpectVal}(\hat O,\psi) = \left<\psi\right|\hat O\left|\psi\right>
\end{eqnarray}

### エネルギー期待値の最小化

import scipy.optimize
import matplotlib.pyplot as plt

def cost(phi):
	"""""
	でたらめな量子状態に対する期待値を出力
	"""""
    return ExpectVal(M, TwoQubitPQC(phi))

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

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

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()

`scipy.optimize.minimize(cost, init,method='Powell',callback=callback)`
costという関数を，initという($\phi$の)初期条件からスタートして，Powell法を用いて，callback関数を使って最適化の経過を記録する．$\vec\phi$はPowell法を用いて更新され，最適化される．

> 変数の詳細については[scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html)を参照のこと．
> callbackとは，scipy.optimize.minimize関数の各Iteration後に呼び込まれる関数である．

## 厳密解との比較
$\hat H$を対角化して得られる厳密な固有値との比較を行う．

import scipy.linalg
l, P = scipy.linalg.eigh(M)
print(l[0]) #最小固有値
print(cost(res.x)) #VQEの結果 ## resは最適化したもので，xは変数組を表している．

出力

最小固有値; -2.8626207640766816
VQEの結果; -2.8529171919339973

これらの結果は近いので，VQEは精度の良い計算アルゴリズムだと考えられる．

## 出力のノイズの影響の検証
量子コンピュータによる計算は，ノイズによる影響を受けやすい．そのため，どの程度のノイズまでなら影響が出ないかを調べておくことは重要である．
### ノイズ付きのエネルギー期待値の定義

def cost_with_noise(phi):
    return ExpectVal(M,TwoQubitPQC(phi))+np.random.normal(0,0.1)

def callback(phi):
    global cost_val
    cost_val.append(cost_with_noise(phi))

ノイズはガウス分布となっているとした．

### ノイズ付きのエネルギー期待値の最小化

cost_val=[] # コスト関数の履歴
init = np.random.rand(6)
callback(init)
res = scipy.optimize.minimize(cost_with_noise, init,
                              method='Powell',
                              callback=callback)
plt.plot(cost_val)
plt.xlabel("iteration")
plt.ylabel("energy expectation value")
plt.show()
print(cost_with_noise(res.x))