前回に引き続き、量子化学計算における量子計算アルゴリズム「VQE(量子固有値変分法ソルバー)」を実装する際に必要となる「Bravey-Kitaev変換」について、pythonライブラリ「OpenFermion」を用いた計算を行ったのでメモを残して置きます。
Bravey-Kitaev変換については、
「The Bravyi-Kitaev transformation for quantum computation of electronic structure」(J. Chem. Phys 2012)
に詳しい理論的背景、および量子計算での具体例(H$_2$分子の系)が載っています。
OpenFermionについては、公式GitHubを参照してインストールおよび計算を行いました。OpenFermionはPython2.7~3.6に対応していますが、今回はPython3.6上でコード実装を行っています。
やること
前回は、電子系のHamiltonian
\begin{eqnarray}
\mathcal{H}(R) &=& \sum_{p,q} h_{pq} c^{\dagger}_p c_q + \frac{1}{2} \sum_{p,q,r,s} h_{pqrs} c^{\dagger}_p c^{\dagger}_q c_r c_s + \sum_{i,j>i} \frac{Z_i Z_j}{|R_i - R_j|}
\end{eqnarray}
の係数$h_{pq}$と$h_{pqrs}$を計算しました。ここで、量子コンピューターでは実装されている量子ビットに作用・観測を行うことで量子計算を行います。つまり、現在Hamiltonianにあるような電子軌道に対する生成消滅演算子$c^{\dagger},c$ (Fermionic operator)は扱うことは出来ないので、これらを量子ビットに作用するqubit operatorへと変換する必要が出てきます。
Googleの開発したPythonライブラリ「OpenFermion」では、このqubit operatorへのmappingを自動で行ってくれます。自分で実装するにはかなり重い作業になりそうで絶望していたので、開発者様には非常に感謝しております。
今回は、H$_2$分子に関して前回計算した結果
\begin{eqnarray}
\mathcal{H}^{(1)}
&=& \sum_{pq} h_{pq}c^{\dagger}_{p}c_{q} \\
&=& h_{00}c^{\dagger}_{00} + h_{11}c^{\dagger}_{11} + h_{22}c^{\dagger}_{22} + h_{33}c^{\dagger}_{33} \\
\mathcal{H}^{(2)}
&=& \frac{1}{2} \sum_{p,q,r,s} h_{pqrs} c^{\dagger}_p c^{\dagger}_q c_r c_s \\
&=& h_{0110} \, c_0^{\dagger} c_1^{\dagger} c_1 c_0 +
h_{2332} \, c_2^{\dagger} c_3^{\dagger} c_3 c_2 +
h_{0330} \, c_0^{\dagger} c_3^{\dagger} c_3 c_0 +
h_{1221} \, c_1^{\dagger} c_2^{\dagger} c_2 c_1 \\
&+&
(h_{0220}-h_{0202}) \, c_0^{\dagger} c_2^{\dagger} c_2 c_0 +
(h_{1331}-h_{1313}) \, c_1^{\dagger} c_3^{\dagger} c_3 c_1 \\
&+&
h_{0132} \, (c_0^{\dagger} c_1^{\dagger} c_3 c_2 + c_2^{\dagger} c_3^{\dagger} c_1 c_0) +
h_{0312} \, (c_0^{\dagger} c_3^{\dagger} c_1 c_2 + c_2^{\dagger} c_1^{\dagger} c_3 c_0) \\
\end{eqnarray}
(論文引用)
から、OpenFermion上でBravy-Kitaev変換を行って論文中の式(79)を再現することを目的とします。
Bravyi-Kitaev変換
fermionic operatorからqubit operatorへの変換には現在2通りの方法が提唱されています。「Jordan-Wigner変換」と「Bravyi-Kitaev変換」です。最近までは、スピン模型の分野で生み出されたJordan-Wigner変換が量子計算の分野にも流用され使われていましたが、量子コンピューターでの計算にはBravyi-Kitaev変換を使ったほうが効率的であることが判明しました。それぞれの変換から実際に量子計算を行う際、ハードウェアに要求されるqubit数やゲート数に関して$O(n)$ vs $O(\mathrm{log},n)$ぐらいのスケーリングで有利らしいです。
理論はぶっ飛ばして具体的な変換方法について少し書いておきます。論文中の式(63)~にあるように
\begin{eqnarray}
c^{\dagger}_0 c_0 &=& \frac{1}{2} ({\bf 1} - Z_0) \\
c^{\dagger}_1 c_1 &=& \frac{1}{2} ({\bf 1} - Z_0 Z_1) \\
&・& \\
&・& \\
&・&
\end{eqnarray}
みたいな感じでそれぞれの項をqubit operatorに変換していけます。手計算でやるには少しややこしいですが、決まった規則通り変換式が求まるので、プログラミング上で実装することが出来ます。
OpenFermion
公式GitHubにあるように、git cloneからpipでインストール出来ます。cloneしたディレクトリ内の OpenFermion/example/openfermion_demo.ipynb にわかりやすいチュートリアルが載っています。
fermion operatorは次のようにコード上で定義して扱うことが出来ます。
下のコードのoutputにある1.0 [3^ 1]は
a^{\dagger}_3 a_1
を表します。
from openfermion.ops import FermionOperator
my_term = FermionOperator(((3, 1), (1, 0)))
print(my_term)
my_term = FermionOperator('3^ 1')
print(my_term)
1.0 [3^ 1]
1.0 [3^ 1]
複数項を扱うときは+演算子で次々に付け加えたりしていけます。
from openfermion.ops import FermionOperator
term_1 = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)
my_operator = term_1 + term_2
print(my_operator)
(1+2j) [4^ 3^ 9 1] +
-1.7 [3^ 1]
Bravey-Kitaev変換は、次のように簡単に行うことが出来ます。
from openfermion.ops import FermionOperator
from openfermion.transforms import jordan_wigner, bravyi_kitaev
from openfermion.utils import eigenspectrum, hermitian_conjugated
fermion_operator = FermionOperator('0^ 0', 1.0)
print(fermion_operator)
# Transform to qubits under the Bravyi-Kitaev transformation and print its spectrum.
bk_operator = bravyi_kitaev(fermion_operator)
print("")
print(bk_operator)
1.0 [0^ 0]
(0.5+0j) [] +
(-0.5+0j) [Z0]
H2分子でのBravyi-Kitaev変換
実際にH$_2$分子についてBravyi-Kitaev変換を行ってみました。
from openfermion.ops import FermionOperator
from openfermion.transforms import jordan_wigner, bravyi_kitaev
from openfermion.utils import eigenspectrum, hermitian_conjugated
h_00 = h_11 = -1.252477
h_22 = h_33 = -0.475934
h_0110 = h_1001 = 0.674493
h_2332 = h_3323 = 0.697397
h_0220 = h_0330 = h_1221 = h_1331 = h_2002 = h_3003 = h_2112 = h_3113 = 0.663472
h_0202 = h_1313 = h_2130 = h_2310 = h_0312 = h_0132 = 0.181287
# Initialize an operator.
fermion_operator = FermionOperator('0^ 0', h_00)
fermion_operator += FermionOperator('1^ 1', h_11)
fermion_operator += FermionOperator('2^ 2', h_22)
fermion_operator += FermionOperator('3^ 3', h_33)
fermion_operator += FermionOperator('0^ 1^ 1 0', h_0110)
fermion_operator += FermionOperator('2^ 3^ 3 2', h_2332)
fermion_operator += FermionOperator('0^ 3^ 3 0', h_0330)
fermion_operator += FermionOperator('1^ 2^ 2 1', h_1221)
fermion_operator += FermionOperator('0^ 2^ 2 0', h_0220-h_0202)
fermion_operator += FermionOperator('1^ 3^ 3 1', h_1331-h_1313)
fermion_operator += FermionOperator('0^ 1^ 3 2', h_0132)
fermion_operator += FermionOperator('2^ 3^ 1 0', h_0132)
fermion_operator += FermionOperator('0^ 3^ 1 2', h_0312)
fermion_operator += FermionOperator('2^ 1^ 3 0', h_0312)
print(fermion_operator)
# Transform to qubits under the Bravyi-Kitaev transformation and print its spectrum.
bk_operator = bravyi_kitaev(fermion_operator)
print('')
print(bk_operator)
-1.252477 [0^ 0] +
-1.252477 [1^ 1] +
-0.475934 [2^ 2] +
-0.475934 [3^ 3] +
0.674493 [0^ 1^ 1 0] +
0.697397 [2^ 3^ 3 2] +
0.663472 [0^ 3^ 3 0] +
0.663472 [1^ 2^ 2 1] +
0.482185 [0^ 2^ 2 0] +
0.482185 [1^ 3^ 3 1] +
0.181287 [0^ 1^ 3 2] +
0.181287 [2^ 3^ 1 0] +
0.181287 [0^ 3^ 1 2] +
0.181287 [2^ 1^ 3 0]
(-0.8126100000000005+0j) [] +
(0.17120100000000002+0j) [Z0] +
(0.17120100000000002+0j) [Z0 Z1] +
(-0.22279649999999998+0j) [Z2] +
(-0.22279649999999998+0j) [Z1 Z2 Z3] +
(0.16862325+0j) [Z1] +
(0.17434925+0j) [Z1 Z3] +
(0.165868+0j) [Z0 Z1 Z2 Z3] +
(0.165868+0j) [Z0 Z1 Z2] +
(0.12054625+0j) [Z0 Z2] +
(0.12054625+0j) [Z0 Z2 Z3] +
(0.04532175+0j) [Y0 Z1 Y2 Z3] +
(0.04532175+0j) [X0 Z1 X2] +
(0.04532175+0j) [X0 Z1 X2 Z3] +
(0.04532175+0j) [Y0 Z1 Y2]
論文中の式(79):
\begin{eqnarray}
\mathcal{H}_{BK}
&=& -0.81261 \, {\bf 1} + 0.171201 \, Z_0 + 0.16862325 \, Z_1 - 0.2227965 \, Z_2 + 0.171201 \, Z_1 Z_0 \\
&+& 0.12054625 \, Z_2 Z_0 + 0.17434925 \, Z_3 Z_1 + 0.04532175 \, X_2 Z_1 X_0 + 0.04532175 \, Y_2 Z_1 Y_0 \\
&+& 0.165868 \, Z_2 Z_1 Z_0 + 0.12054625 \, Z_3 Z_2 Z_0 - 0.2227965 \, Z_3 Z_2 Z_1 \\
&+& 0.04532175 \, Z_3 X_2 Z_2 X_0 + 0.04532175 \, Z_3 Y_2 Z_1 Y_0 + 0.165868 \, Z_3 Z_2 Z_1 Z_0
\end{eqnarray}
と確かに一致しました！
(補足)
PRX(2016)論文では同じくH$_2$分子を扱っているので、この式と同じ形の式(A5)がBK変換後に出現しています。量子計算のスタートにHartree-Fock stateを用いていることによってqubit1と3がflipしないらしく(イマイチ理由がわからないので偉い人教えてください...)、1と3のインデックスを持つ演算子を恒等演算子に置き換えてさらに簡単なHamiltoninで表現しています。
\begin{eqnarray}
\mathcal{H}_{BK}
&=& g_0 \, {\bf 1} + g_1 \, Z_0 + g_2 \, Z_1 + g_3 \, Z_0 Z_1 + g_4 \, X_0 X_1 + g_5 \, Y_0 Y_1 \\
&=& -0.4804 \, {\bf 1} + 0.3435 \, Z_0 + -0.4347 \, Z_1 + 0.5716 \, Z_0 Z_1 + 0.0910 \, X_0 X_1 + 0.0910 \, Y_0 Y_1 \\
\end{eqnarray}
論文中Table 1を参照して、R=0.75の時の係数gを代入しました。今回の結果においても同様の簡略化を行うと、だいたい同じような係数の式となります。