Quantum Native Dojoの5-1. Variational Quantum Eigensolver（VQE）アルゴリズムはpythonのnumpyで実装されていますが、Qulacsを使って同じVQEを実装してみました。
(実装が間違っていたら教えてください。)
まずコードの解説をします。
次にQuantum Native Dojoの5-1で使われているものと同じ変分回路で計算した結果を報告します。
また変分回路でCNOTゲートを作用させなかった場合やパラメータの数を増やした場合の計算結果も報告します。
最後にCNOTゲートがなぜ必要なのか考察してみました。
コードの解説
コード全体は以下のgithubリポジトリにあります。
以下で詳しく説明していきます。
モジュールのインポート
from qulacs import Observable
from qulacs import QuantumState, QuantumCircuit
from qulacs.gate import RX, RZ, CNOT
from scipy.optimize import minimize
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
変分回路の作成
def make_ansatz_circuit(phi, n_bits):
circuit = QuantumCircuit(n_bits)
circuit.add_gate(RX(0, phi[0]))
circuit.add_gate(RZ(0, phi[1]))
circuit.add_gate(RX(1, phi[2]))
circuit.add_gate(RZ(1, phi[3]))
circuit.add_gate(CNOT(0, 1))
circuit.add_gate(RX(1, phi[4]))
circuit.add_gate(RZ(1, phi[5]))
return circuit
phiに最適化するパラメータのリスト、n_bitsに量子ビット数を入力します。
ハミルトニアンの作成
第２量子化を行いパウリ演算しで表現した水素分子のハミルトニアンを準備します。
hamiltonian = [(-3.8505/2, ""), \
(-0.2288/2, "X 1"), \
(-1.0466/2, "Z 1"),\
(-0.2288/2, "X 0"), \
(0.2613/2, "X 0 X 1"), \
(0.2288/2, "X 0 Z 1"), \
(-1.0466/2, "Z 0"), \
(0.2288/2, "Z 0 X 1"), \
(0.2356/2, "Z 0 Z 1")]
QulacsのObservableに上記の演算子を追加していく。
n_bits = 2
observable = Observable(n_bits)
for (coef, Pauli_string) in hamiltonian :
observable.add_operator(coef, Pauli_string)
コスト関数の作成
初期状態にmake_ansatz_circuitで作成した変分回路を作用させてobservableの期待値を計算します。
def cost(phi):
state = QuantumState(n_bits)
state.set_zero_state()
circuit = make_ansatz_circuit(phi, n_bits)
circuit.update_quantum_state(state)
return observable.get_expectation_value(state)
最適化を行う過程で上で計算した期待値をcost_valに順次保存していくための関数です。
def callback(phi):
global cost_val
cost_val.append(cost(phi))
最適化
シミュレーターなので意味があるかどうかは分かりませんが同じ計算を100回を行いその平均値を最後に出力しています。
グラフは最後のiterationの結果だけを表示させています。
n_params = 6
np.random.seed(1111)
n_iter = 100
energy_list = []
for _ in range(n_iter):
cost_val = []
phi = np.random.random(n_params)
callback(phi)
res = minimize(cost, phi, method='Powell', callback=callback)
energy_list.append(cost(res.x))
plt.plot(cost_val)
plt.xlabel("iteration")
plt.ylabel("energy expectation value")
plt.show()
print("Exact Energy: ", -2.8626207640766816)
print("VQE Energy: {}".format(np.mean(energy_list)))
計算結果
出力結果は以下のようになります。
Exact Energy: -2.8626207640766816
VQE Energy: -2.8624166416549484
厳密対角化を行った結果と小数第３位まで一致しています。
計算結果にノイズが入った場合
costとcallbackを以下のように変更します。
cost_noiseでは最後にノイズが追加されています。
def cost_noise(phi):
state = QuantumState(n_bits)
state.set_zero_state()
circuit = make_ansatz_circuit(phi, n_bits)
circuit.update_quantum_state(state)
return observable.get_expectation_value(state) + np.random.normal(0, 0.01)
def callback_noise(phi):
global cost_val
cost_val.append(cost_noise(phi))
最適化は上記で紹介したコードのcostの部分をcost_noise、callbackの部分をcallback_noiseに変更すれば実行できます。
結果は
Exact Energy: -2.8626207640766816
VQE Energy: -2.850086422671339
小数第２位以降で一致しないようになりました。
パラメータを増やした場合
ノイズなし
make_ansatz_circuitを以下のように変更します。
def make_ansatz_circuit(phi, n_bits):
circuit = QuantumCircuit(n_bits)
circuit.add_gate(RX(0, phi[0]))
circuit.add_gate(RZ(0, phi[1]))
circuit.add_gate(RX(1, phi[2]))
circuit.add_gate(RZ(1, phi[3]))
circuit.add_gate(CNOT(0, 1))
circuit.add_gate(RX(1, phi[4]))
circuit.add_gate(RZ(1, phi[5]))
circuit.add_gate(RX(0, phi[6]))
circuit.add_gate(RZ(0, phi[7]))
circuit.add_gate(CNOT(0, 1))
return circuit
最適化は上記で紹介したコードのn_params = 6の部分をn_params = 8に変更すれば実行できます。
結果は
Exact Energy: -2.8626207640766816
VQE Energy: -2.8626194809589203
小数第4位まで一致するようになりました！
パラメータの数が6個の変分回路では-2.8624166416549484だったので精度が向上しています。
ノイズあり
上と同じ変更を施したコードでノイズが入った時に使ったコードをn_params = 8に変更すれば実行できます。
結果は
Exact Energy: -2.8626207640766816
VQE Energy: -2.8457532192348656
パラメータ数が6個でノイズが入った場合の結果が-2.850086422671339だったので逆に精度が下がりました。
CNOTゲートなし、ノイズなし
最後にCNOTゲートなし、ノイズなしで計算を行ってみます。
make_ansatz_circuitからCNOTゲートを削除します。
def make_ansatz_circuit(phi, n_bits):
circuit = QuantumCircuit(n_bits)
circuit.add_gate(RX(0, phi[0]))
circuit.add_gate(RZ(0, phi[1]))
circuit.add_gate(RX(1, phi[2]))
circuit.add_gate(RZ(1, phi[3]))
circuit.add_gate(RX(1, phi[4]))
circuit.add_gate(RZ(1, phi[5]))
return circuit
結果は
Exact Energy: -2.8626207640766816
VQE Energy: -2.85404959554127
CNOTゲートがある場合は-2.8624166416549484だったので精度が下がりました。
CNOTゲートを作用させた方が良いことが分かりました。
考察
CNOTゲートを作用させる場合とさせない場合でどういった違いがあるのかを考察したいと思います。
今回は2量子ビットの状態なので具体的に手計算を行いながら分析したいと思います。
$$
\def\bm#1{\boldsymbol{#1}}
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\brakett#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
\def\braket#1#2#3{\bra{#1}#2\ket{#3}}
$$
初期状態を$\ket{\psi_{init}}=\ket{00}$として、make_ansatz_circuitの演算子を順番に作用させていきます。
(左から0ビット、1ビットと数えることにします。)
def make_ansatz_circuit(phi, n_bits):
circuit = QuantumCircuit(n_bits)
circuit.add_gate(RX(0, phi[0]))
circuit.add_gate(RZ(0, phi[1]))
circuit.add_gate(RX(1, phi[2]))
circuit.add_gate(RZ(1, phi[3]))
circuit.add_gate(CNOT(0, 1))
circuit.add_gate(RX(1, phi[4]))
circuit.add_gate(RZ(1, phi[5]))
return circuit
まず回転ゲートは以下のように定義されています。
\begin{align}
R_{X_i}(\theta) &= e^{-\frac{i\theta}{2} X_i} = \cos{(\frac{\theta}{2})} I_i - i \sin{(\frac{\theta}{2})} X_i \\
R_{Z_i}(\theta) &= e^{-\frac{i\theta}{2} Z_i} = \cos{(\frac{\theta}{2})} I_i - i \sin{(\frac{\theta}{2})} Z_i
\end{align}
$a(\theta) = \cos{(\frac{\theta}{2})}, b(\theta) = - i \sin{(\frac{\theta}{2})}$とおくと、
\begin{align}
R_{X_i}(\theta) &= a(\theta) I_i + b(\theta) X_i \\
R_{Z_i}(\theta) &= a(\theta) I_i + b(\theta) Z_i
\end{align}
と表すことが出来ます。
1ビット状態$\ket{0}, \ket{1}$にRX, RZを順番に作用させて得られる状態をあらかじめ計算しておきます。
まずは$\ket{0}$に対して
\begin{align}
R_X(\theta_0)\ket{0} &= a(\theta_0)\ket{0} + b(\theta_0)\ket{1} \\
R_Z(\theta_1)R_X(\theta_0)\ket{0} &= a(\theta_1)\biggl( a(\theta_0)\ket{0} + b(\theta_0)\ket{1} \biggr) + b(\theta_1) \biggl(a(\theta_0)\ket{0} - b(\theta_0)\ket{1} \biggr) \\
&= \biggl(a(\theta_0)a(\theta_1) + a(\theta_0)b(\theta_1) \biggr)\ket{0} + \biggl(b(\theta_0)a(\theta_1) - b(\theta_0)b(\theta_1) \biggr)\ket{1} \\
& \equiv A(\theta_0, \theta_1)\ket{0} + B(\theta_0, \theta_1)\ket{1}
\end{align}
となります。$\ket{1}$に対して、
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\begin{align}
R_X(\theta_0)\ket{1} &= a(\theta_0)\ket{1} + b(\theta_0)\ket{0} \\
R_Z(\theta_1)R_X(\theta_0)\ket{1} &= a(\theta_1)\biggl( a(\theta_0)\ket{1} + b(\theta_0)\ket{0} \biggr) + b(\theta_1) \biggl(-a(\theta_0)\ket{1} + b(\theta_0)\ket{0} \biggr) \\
&= \biggl(b(\theta_0)a(\theta_1) + b(\theta_0)b(\theta_1) \biggr)\ket{0} + \biggl(a(\theta_0)a(\theta_1) - a(\theta_0)b(\theta_1) \biggr)\ket{1} \\
& \equiv C(\theta_0, \theta_1)\ket{0} + D(\theta_0, \theta_1)\ket{1}
\end{align}
となります。
簡単のために係数$A,B,C,D$は$A(\theta_i, \theta_j)=A_{ij}$とそれぞれ略記することにします。
circuit.add_gate(~))で演算子が作用されるたびに状態を$\ket{\psi_{i-1}}\to\ket{\psi_i}$と更新して表記していきます。
最初は$\ket{\psi_0}=R_{X_0}(\theta_0)\ket{\psi_{init}}$と表すことにします。
CNOTゲートの前までは以下のように状態が更新されていきます。
\begin{align}
\ket{\psi_0} &= R_{X_0}(\theta_0)\ket{\psi_{init}} \\
\ket{\psi_1} &= R_{Z_0}(\theta_1)\ket{\psi_{0}} \\
\ket{\psi_2} &= R_{X_1}(\theta_2)\ket{\psi_{1}} \\
\ket{\psi_3} &= R_{Z_1}(\theta_3)\ket{\psi_{2}}
\end{align}
CNOTゲートの手前まではそれぞれにビットに独立して演算子が作用するので
\begin{align}
\ket{\psi_3} &= R_{Z_0}(\theta_1)R_{X_0}(\theta_0)\ket{0} \otimes R_{Z_1}(\theta_3)R_{X_1}(\theta_2)\ket{0} \\
&= \biggl(A_{01}\ket{0} + B_{01}\ket{1} \biggr) \otimes \biggl(C_{23}\ket{0} + D_{23}\ket{1} \biggr) \\
&= A_{01}C_{23}\ket{00} + A_{01}D_{23}\ket{01} + B_{01}C_{23}\ket{10} + B_{01}D_{23}\ket{11}
\end{align}
となります。
CNOTゲートを作用させると、
\begin{align}
\ket{\psi_4} &= CNOT(0, 1)\ket{\psi_{3}} \\
&= A_{01}C_{23}\ket{00} + A_{01}D_{23}\ket{01} + B_{01}C_{23}\ket{11} + B_{01}D_{23}\ket{10}\\
&= \biggl(A_{01}C_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \ket{0} + \biggl(A_{01}D_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \ket{1}
\end{align}
になります。
最後に1ビット目に$R_X, R_Z$を続けて作用させると
\begin{align}
\ket{\psi_5} &= R_{X_1}(\theta_4)\ket{\psi_{4}} \\
\ket{\psi_6} &= R_{Z_1}(\theta_5)\ket{\psi_{5}}
\end{align}
なので最終的な状態は、
\begin{align}
\ket{\psi_6} &= \biggl(A_{01}C_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \biggl(A_{45}\ket{0} + B_{45}\ket{1} \biggr) \\
& \;\;\; + \biggl(A_{01}D_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \biggl(C_{45}\ket{0} + D_{45}\ket{1} \biggr)\\
&= \biggl(A_{01}C_{23}A_{45} +A_{01}D_{23}C_{45}\biggr)\ket{00} + \biggl(A_{01}C_{23}B_{45} +A_{01}D_{23}D_{45}\biggr)\ket{01} \\
& \;\;\; + \biggl(B_{01}D_{23}A_{45} +B_{01}C_{23}C_{45}\biggr)\ket{10} + \biggl(B_{01}D_{23}B_{45} +B_{01}C_{23}D_{45}\biggr)\ket{11} \\
\end{align}
となります。
もしCNOTゲートを作用させなければ、
\begin{align}
\ket{\psi^{\prime}_4} &= \ket{\psi_{3}} \\
&= A_{01}C_{23}\ket{00} + A_{01}D_{23}\ket{01} + B_{01}C_{23}\ket{10} + B_{01}D_{23}\ket{11}\\
&= \biggl(A_{01}C_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \ket{0} + \biggl(A_{01}D_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \ket{1}
\end{align}
より、
\begin{align}
\ket{\psi^{\prime}_6} &= \biggl(A_{01}C_{23}\ket{0} + B_{01}C_{23}\ket{1} \biggr) \otimes \biggl(A_{45}\ket{0} + B_{45}\ket{1} \biggr) \\
& \;\;\; + \biggl(A_{01}D_{23}\ket{0} + B_{01}D_{23}\ket{1} \biggr) \otimes \biggl(C_{45}\ket{0} + D_{45}\ket{1} \biggr)\\
&= \biggl(A_{01}C_{23}A_{45} +A_{01}D_{23}C_{45}\biggr)\ket{00} + \biggl(A_{01}C_{23}B_{45} +A_{01}D_{23}D_{45}\biggr)\ket{01} \\
& \;\;\; + \biggl(B_{01}C_{23}A_{45} +B_{01}D_{23}C_{45}\biggr)\ket{10} + \biggl(B_{01}C_{23}B_{45} +B_{01}D_{23}D_{45}\biggr)\ket{11} \\
\end{align}
$\ket{\psi_6}, \ket{\psi^{\prime}_6}$を比較すると$\ket{10}, \ket{11}$の真ん中に挟まれている係数($\theta_2, \theta_3$に関する係数)が異なっています。
CNOTゲートを作用させることで、$\theta_2, \theta_3$に関する係数が、$\ket{00}, \ket{01}$はそのままで、$\ket{10}, \ket{11}$だけ反転していることが分かります。
これは単純に回転ゲートを作用させていくだけでは得られない状態でCNOTゲートを作用させることでより複雑な(表現力の高い？)状態を生成していることが分かります。
このような状態を利用して期待値を計算することでより精度の高いエネルギー値が得られているのだと思います。
変分回路を作る統一的な指針みたいなものはあるのでしょうか？
ディープラーニングのニューラルネットワークと同じように層の数や活性化関数の選び方がある程度経験的に基づくのと同様に変分回路も試行錯誤しながら作成していくのでしょうか。