LoginSignup
19
7

More than 3 years have passed since last update.

量子コンピューターAdvent Calendar 2019Day 11
@nori_autumn(nori)

SSVQEを用いて励起状態との遷移行列を計算してみた。

Last updated at Posted at 2019-12-10

概要

量子系のハミルトニアンの基底状態を求める手法としてVariational Quantum Eigensolver(VQE)と呼ばれるアルゴリズムが利用されています。
VQEは主にハミルトニアンの基底状態を求めることに特化したアルゴリズムになっています。
励起状態を求める手法としてSubspace-Search VQE(SSVQE)と呼ばれる手法が次の論文で提案されています。

この記事ではSSVQEを説明した後に、Qunasysが提供している量子コンピューターのシミュレーターのQulacsを利用して実際に分子の励起状態について計算してみたいと思います。
また上記の論文で説明されているように励起状態と基底状態の遷移行列(要素)を計算してみたいと思います。

基本的なコードはこちら(6-3. 励起状態の探索手法 (subspace-search variational quantum eigensolver))のコードをベースにしました。

VQEとは

VQEはVariational Quantum Eigensolverの略称で日本語では変分量子固有法と呼ばれています
その名の通り変分法を量子力学に適用してシュレディンガー方程式の最小固有値を得るための手法です。

まず量子力学で使われる変分法について説明します。
我々が解きたい方程式は以下の時間に依存しないシュレディンガー方程式になります。

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
\mathcal{H}\ket{\Psi} = E \ket{\Psi}
$$

これに対して以下の汎関数$E[\Psi]$を定義します：

\begin{align}
E[\Psi] := \frac{\bra{\Psi}\mathcal{H} \ket{\Psi}}{\braket{\Psi}{\Psi}}
\end{align}

汎関数$E[\Psi]$を最小化することが基底状態のエネルギー(最小固有値)を求めることに対応しています。

基底状態のエネルギーを$E_0$と表し、エネルギーの小さい順に$E_0\le E_1\le E_2\le\cdots$として、$\ket{\Psi}$を$\mathcal{H}$のエネルギー固有状態$\ket{\psi_n}$を次のように表します。

\begin{align}
\mathcal{H} \ket{\psi_n} &= E_n \ket{\psi_n}\\
 \braket{\psi_n}{\psi_m} &= \delta_{nm}
\end{align}

この$\ket{\psi_n}$で$\ket{\Psi}$を

\begin{align}
\ket{\Psi} = \sum_{n=0}c_n\ket{\psi_n}
\end{align}

と展開して、汎関数$E[\Psi]$に代入すると、

\begin{align}
E[\Psi] &= \frac{\bra{\Psi}\mathcal{H} \ket{\Psi}}{\braket{\Psi}{\Psi}}\\
&= \frac{\sum_{n,m=0}\delta_{nm}c^{*}_nc_mE_n}{\sum_{n,m=0}\delta_{nm}c^{*}_nc_m}\\
&= \frac{\sum_{n=0}E_n|c_n|^2}{\sum_{n=0}|c_n|^2} \ge \frac{E_0\sum_{n=0}|c_n|^2}{\sum_{n=0}|c_n|^2} = E_0\\
\therefore E[\Psi] &\ge E_0
\end{align}

また$\braket{\Psi}{\Psi}=1$の拘束条件を取り入れた

\begin{align}
L[\Psi] := \bra{\Psi}\mathcal{H} \ket{\Psi} - E\biggl(\braket{\Psi}{\Psi}-1 \biggr)
\end{align}

定義し、変分操作$\ket{\Psi}\to \ket{\Psi} + \delta\ket{\Psi}$を行うと、

\begin{align}
\delta L[\Psi] = 0 \Leftrightarrow \mathcal{H}\ket{\Psi} = E \ket{\Psi}
\end{align}

が得られます。
これは$\delta E[\Psi]=0$の極値ではシュレディンガー方程式が成立していることを意味します。
(参考：http://www.th.phys.titech.ac.jp/~muto/lectures/QMII09/QMII09_chap16.pdf)

上記で説明しているように$\ket{\Psi}$をエネルギー固有状態$\ket{\psi_n}$で展開していましたが、そもそもそれ自体が難しくて目標でもあります。
$\ket{\Psi}$を$\ket{\Psi}=\ket{\Psi(\boldsymbol{\theta})}$のように$\boldsymbol{\theta}$でパラメトライズして、$\boldsymbol{\theta}$の値を色々変更しながら$E[\Psi(\boldsymbol{\theta})]$が最小になる値を探索します。

$\ket{\Psi(\boldsymbol{\theta})}$のことを試行関数やVQE ansatzと呼びます。
VQEではこの試行関数を量子回路上で作りだして、$E[\Psi(\boldsymbol{\theta})]$を計算します。

具体的には$\theta_i$でパラメトライズされた1または2量子ビット間に作用する量子ゲート$U_i(\theta_i)$を組み合わせた$U(\boldsymbol{\theta}) = U_1(\theta_1)U_2(\theta_2)\cdots U_N(\theta_N)$を量子回路上に実装して、

\begin{align}
\ket{\Psi(\boldsymbol{\theta})} = U(\boldsymbol{\theta})\ket{\psi_{ref}}
\end{align}

という状態を作り出します。
$\ket{\psi_{ref}}$は$U(\boldsymbol{\theta})$を作用させる前の初期状態になります。
$\ket{\Psi(\boldsymbol{\theta})}$を作り出した後はその状態でハミルトニアンの期待値$E[\Psi(\boldsymbol{\theta})]$を計算します。

\begin{align}
E[\Psi(\boldsymbol{\theta})] = \bra{\psi_{ref}} U^{\dagger}(\boldsymbol{\theta}) \mathcal{H} U(\boldsymbol{\theta})\ket{\psi_{ref}}
\end{align}

期待値を計算した後は古典コンピューターで$\boldsymbol{\theta}$の値を最適化します。
量子回路のパラメータを最適化した$\boldsymbol{\theta}$に更新して再び$E[\Psi(\boldsymbol{\theta})]$を計算します。
VQEでは$E[\Psi(\boldsymbol{\theta})]$の値が収束するまでこのプロセスを繰り返していきます。

VQEをまとめると次のようになります。

  • 量子コンピューター：$E[\Psi(\boldsymbol{\theta})]$を計算する。
  • 古典コンピューター：$\boldsymbol{\theta}$の最適化・更新を行う。

VQEを行うまでの手続き

上の記事でも説明されているようにVQEを行うまでに次の手続きが必要になります。

  1. 電子系のハミルトニアンを定義する。
  2. ボルン・オッペンハイマー近似を行う。
  3. 第2量子化を行う。
  4. Jordan-Wigner変換(Bravyi-Kitaev変換)を行う。
  5. VQEを行う。

1~4で行っていることはまず1~3で解きたいハミルトニアンを定義して、4でそのハミルトニアンを量子コンピューターで解ける形式に変換しています。

実装する際にはこれらの一連の流れをスクラッチで書く必要はなく、OpenFermionと呼ばれるライブラリを使えば簡単に量子コンピューターで解くためのデータを準備することができます。

この記事では水素原子のエネルギー状態をSSVQEを用いて計算していきます。
水素原子のハミルトニアンに対して上記の1~4のプロセスを解説した記事があるので詳しく知りたい方は以下を参照してください。

以降では1~4のプロセスの概要について簡単に説明していきます。

1. ハミルトニアンを定義する。

まず解きたい方程式は次の時間に依存しないシュレディンガー方程式になります。

$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
\mathcal{H}\ket{\Psi} = E \ket{\Psi}
$$

$\mathcal{H}$はハミルトニアンといい、考えている系のエネルギーを決定する演算子です。
$\ket{\Psi}$は波動関数といい、ハミルトニアンが作用する関数になります。

原子核の数をM, 電子の数をNとした時の電子系のハミルトニアンは次のように表せます。

\begin{align}
\mathcal{H} &= -\sum_i\frac{\hbar^2}{2m_e}\nabla^2_i +\frac{1}{2}\sum_{i\neq j}\frac{e^2}{4\pi \epsilon_0}\frac{1}{|\boldsymbol{r}_i - \boldsymbol{r}_j|} \\
&-\sum_I\frac{\hbar^2}{2M_I}\nabla^2_I + \frac{1}{2}\sum_{I\neq J}\frac{e^2}{4\pi \epsilon_0}\frac{1}{|\boldsymbol{R}_I - \boldsymbol{R}_J|} - \sum_{i,I}\frac{e^2}{4\pi \epsilon_0} \frac{Z_I}{|\boldsymbol{r}_i - \boldsymbol{R}_I|} \tag{1}
\end{align}

$m_e, \boldsymbol{r}_i$はそれぞれ電子の質量と位置ベクトルになります。
$M_I, \boldsymbol{R}_I, Z_I$はそれぞれ原子核の質量、位置ベクトル、原子番号になります。

(1)の第1項は電子の運動エネルギー、第2項は電子間のポテンシャルエネルギー(クーロン相互作用)を表しています。
第3項、第4項はこれの原子核バージョンになります。
最後の第5項は電子と原子核間のポテンシャルエネルギーを表しています。

2. ボルン・オッペンハイマー近似を行う。

(1)のハミルトニアンは電子と原子核の両方の影響をフルに考慮したモデルになっています。
原子核の質量は電子に比べて約1840倍重いことから、原子核は位置ベクトルの原点中心で静止していると仮定して、その周りの電子の振る舞いについて考えるという近似を取り入れます。
この近似のことをボルン・オッペンハイマー近似と言います。
(1)の$\boldsymbol{R}_I$のみに依存している項を無視すると次のようになります。

\begin{align}
\mathcal{H}_e = -\sum_i \frac{\nabla^2_i}{2} + \frac{1}{2}\sum_{i\neq j}\frac{1}{|\boldsymbol{r}_i - \boldsymbol{r}_j|} -\sum_{i, I}\frac{Z_I}{|\boldsymbol{r}_i - \boldsymbol{R}_I|} \tag{2}
\end{align}

(1)で現れていた定数は1ハートリー=$e^2/4\pi \epsilon_0 a_0$($a_0$は1ボーアと呼ばれる原子の単位長さ。)という単位で測り無次元化しています。

3. 第2量子化を行う。

(1)や(2)はいわゆる第1量子化されたハミルトニアンです。
実際に計算するときにはこの表示では扱いにくいことが多いので第2量子化を行います。
(2)には電子の位置ベクトル$\boldsymbol{r}_i$が含まれていたり、電子の波動関数$\ket{\Psi}_e$は

\begin{align}
\ket{\Psi}_e = \sum_{\boldsymbol{r}_1,\boldsymbol{r}_2,\cdots,\boldsymbol{r}_N}\psi(\boldsymbol{r}_1,\boldsymbol{r}_2,\cdots,\boldsymbol{r}_N)\ket{\boldsymbol{r}_1,\boldsymbol{r}_2,\cdots,\boldsymbol{r}_N}
\end{align}

と表すことができます。
各電子がどこにいるのかではなく各電子がどの状態に何個いるのかという点に表示の方法を変える手法を第2量子化と呼ばれています。
考え方としてはまず電子が1つも存在しない真空状態を定義して、$\ket{0}$と表すことにします。
この真空状態をベースにしてある状態pに電子を1つ生成する演算子を$a^{\dagger}_p$と定義します。
逆にある状態pの電子を消滅させる演算子を$a_p$と定義します。
電子は粒子の中でフェルミオンと呼ばれるグループに属しており、同じ状態を取れる電子の数は1つにまで制限されます。
この統計性を上記の生成消滅演算子に反映したものが以下の交換関係になります。

\begin{align}
\{a_p, a^{\dagger}_q\} &= a_pa^{\dagger}_q + a^{\dagger}_qa_p = \delta_{pq} \\
\{a_p, a_q\} &= \{a^{\dagger}_p, a^{\dagger}_q\} = 0
\end{align}

この生成消滅演算子を利用して(2)を書き換えると次のようになります。

\begin{align}
\mathcal{H}_e = \sum_{p,q}h_{pq}a^{\dagger}_pa_q + \sum_{p,q,r,s}h_{pqrs}a^{\dagger}_p a^{\dagger}_qa_r a_s \tag{3}
\end{align}

(3)の導出に関しては入門的な場の量子論の教科書に書かれているのですが、自分がわかりやすいと思ったのは「フェッター・ワレッカの教科書(https://www.amazon.co.jp/dp/B008TVLMPO )」です。
無料で読めるpdfとしてはいかが分りやすかったり、詳しかったりします。

4. Jordan-Wigner変換(Bravyi-Kitaev変換)を行う。

最後に(3)から量子コンピューターで解くための形式に変換します。
Jordan-Wigner変換はフェルミオンの演算子を1/2スピン系のパウリ演算子に変換します。

フェルミオンの演算子は次の交換関係に従っていました。

\begin{align}
\{a_p, a^{\dagger}_q\} &= a_pa^{\dagger}_q + a^{\dagger}_qa_p = \delta_{pq} \\
\{a_p, a_q\} &= \{a^{\dagger}_p, a^{\dagger}_q\} = 0
\end{align}

一方、量子コンピューターの回路図でよく見かけるXゲートやZゲートはスピンのアップ方向とダウン方向のそれぞれに作用するパウリ演算子(パウリ行列)で表されています。

\begin{align}
\sigma_x = \left(
    \begin{array}{cc}
      0 & 1 \\
      1 & 0 
    \end{array}
  \right), \; \sigma_y = \left(
    \begin{array}{cc}
      0 & -i \\
      i & 0 
    \end{array}
  \right), \; \sigma_z = \left(
    \begin{array}{cc}
      1 & 0 \\
      0 & -1 
    \end{array}
  \right)
\end{align}

上式の定義よりパウリ行列は次の交換関係を満たします。

\begin{align}
[\sigma_i, \sigma_j] &= \sigma_i\sigma_j - \sigma_j\sigma_i = 2i\sum^3_{k=1}\epsilon_{ijk}\sigma_k \\ 
\{\sigma_i, \sigma_j\} &= \sigma_i\sigma_j + \sigma_j\sigma_i = 2\delta_{ij}\boldsymbol{I}
\end{align}

$1=x, 2=y, 3=z$と考えてください。
$\epsilon_{ijk}$はエディントン記号になり、(1,2,3)の順番を基準にしています。

このようにフェルミオンとパウリ行列では交換関係が異なるため(3)のハミルトニアンのままでは直接量子コンピューターで計算する子は出来ません。
フェルミオンからパウリ行列の変換を行うのがJordan-Wigner変換になります。

一般的にJordan-Wigner変換を表すと次の形をしているようです。

\begin{align}
a^{\dagger}_i &= \frac{1}{2}\biggl(\sigma_{ix} -i \sigma_{iy} \biggr)\displaystyle{\bigotimes_{j<i}}\;\sigma_{jz}\\
a_i &= \frac{1}{2}\biggl(\sigma_{ix} +i \sigma_{iy} \biggr)\displaystyle{\bigotimes_{j<i}}\;\sigma_{jz}
\end{align}

Bravyi-Kitaev変換はJordan-Wigner変換の計算量を効率化したバージョンになっているようです。
詳しくは以下を参照してください。

変換後のハミルトニアンは次のような形をしています。

\mathcal{H}_e = \sum_{i,j,k,l} h_{i,j,k,l}\otimes_i \boldsymbol{I}_i \otimes_j \sigma_{jx} \otimes_k \sigma_{ky} \otimes_l \sigma_{lz}

では実際に水素分子だとどういった形になっているのかをOpenFermionを使って見てみます。
まずは第2量子化された水素分子のハミルトニアンを準備する。

from openfermion.hamiltonians import MolecularData
from openfermionpsi4 import run_psi4
from openfermion.transforms import get_fermion_operator, jordan_wigner

bond_length = 0.7414
geometry = [["H", [0, 0, 0]],
           ["H", [0, 0, bond_length]]]
basis = "sto-3g"
multiplicity = 1
charge = 0
description = "H2"

molecule = MolecularData(geometry, basis, multiplicity, charge, description)

molecule = run_psi4(molecule)
print("Second Quantized Hamiltonian")
print(molecule.get_molecular_hamiltonian())
Output
Second Quantized Hamiltonian
() 0.7137539933504182
((0, 1), (0, 0)) -1.2524635733533045
((1, 1), (1, 0)) -1.2524635733533045
((2, 1), (2, 0)) -0.4759487154406823
((3, 1), (3, 0)) -0.4759487154406823
((0, 1), (0, 1), (0, 0), (0, 0)) 0.33724438314499877
((0, 1), (0, 1), (2, 0), (2, 0)) 0.09064440411555524
((0, 1), (1, 1), (1, 0), (0, 0)) 0.33724438314499877
((0, 1), (1, 1), (3, 0), (2, 0)) 0.09064440411555524
((0, 1), (2, 1), (0, 0), (2, 0)) 0.09064440411555533
((0, 1), (2, 1), (2, 0), (0, 0)) 0.331734048181356
((0, 1), (3, 1), (1, 0), (2, 0)) 0.09064440411555533
((0, 1), (3, 1), (3, 0), (0, 0)) 0.331734048181356
((1, 1), (0, 1), (0, 0), (1, 0)) 0.33724438314499877
((1, 1), (0, 1), (2, 0), (3, 0)) 0.09064440411555524
((1, 1), (1, 1), (1, 0), (1, 0)) 0.33724438314499877
((1, 1), (1, 1), (3, 0), (3, 0)) 0.09064440411555524
((1, 1), (2, 1), (0, 0), (3, 0)) 0.09064440411555533
((1, 1), (2, 1), (2, 0), (1, 0)) 0.331734048181356
((1, 1), (3, 1), (1, 0), (3, 0)) 0.09064440411555533
((1, 1), (3, 1), (3, 0), (1, 0)) 0.331734048181356
((2, 1), (0, 1), (0, 0), (2, 0)) 0.331734048181356
((2, 1), (0, 1), (2, 0), (0, 0)) 0.09064440411555529
((2, 1), (1, 1), (1, 0), (2, 0)) 0.331734048181356
((2, 1), (1, 1), (3, 0), (0, 0)) 0.09064440411555529
((2, 1), (2, 1), (0, 0), (0, 0)) 0.0906444041155553
((2, 1), (2, 1), (2, 0), (2, 0)) 0.3486968836792783
((2, 1), (3, 1), (1, 0), (0, 0)) 0.0906444041155553
((2, 1), (3, 1), (3, 0), (2, 0)) 0.3486968836792783
((3, 1), (0, 1), (0, 0), (3, 0)) 0.331734048181356
((3, 1), (0, 1), (2, 0), (1, 0)) 0.09064440411555529
((3, 1), (1, 1), (1, 0), (3, 0)) 0.331734048181356
((3, 1), (1, 1), (3, 0), (1, 0)) 0.09064440411555529
((3, 1), (2, 1), (0, 0), (1, 0)) 0.0906444041155553
((3, 1), (2, 1), (2, 0), (3, 0)) 0.3486968836792783
((3, 1), (3, 1), (1, 0), (1, 0)) 0.0906444041155553
((3, 1), (3, 1), (3, 0), (3, 0)) 0.3486968836792783

(2,1)が$c^{\dagger}_2$、(2,0)が$c_0$を表しています。
今回の設定では電子が0~3の4状態を取ることができ、真空状態を$\ket{0000}$と定義して左から順番に状態0を取る電子がいると1、いないと0を表します。
例えば、$\ket{1000}=c^{\dagger}_0\ket{0000}, \ket{1100}=c^{\dagger}_0c^{\dagger}_1\ket{0000}$など。

2行目の((0,1), (0,0))は$c^{\dagger}_0 c_0$を表しており、もし既に状態0に電子がいる($\ket{1000}$)とエネルギーは-1.2524635733533045になることを意味しています。

\begin{align}
c^{\dagger}_0 c_0 \ket{1000} &= c^{\dagger}_0 c_0 c^{\dagger}_0 \ket{0000} = c^{\dagger}_0(1-c^{\dagger}_0c_0)\ket{0000} \\
&= c^{\dagger}_0\ket{0000} - c^{\dagger}_0c_0\ket{0000} = c^{\dagger}_0\ket{0000} = \ket{1000}
\end{align}

(上の式全体には係数$h_{00}=-1.2524635733533045$がかかっています。)
1行目の()は恒等演算子を表しています。右の数値はその係数です。

この第2量子化したハミルトニアンを量子コンピューターで計算するためにJordan-Wigner変換を行います。

print("Hamiltonian peformed by Jordan-Wigner Transformation")
jw_hamiltonian = jordan_wigner(get_fermion_operator(molecule.get_molecular_hamiltonian()))
print(jw_hamiltonian)
Output
Hamiltonian peformed by Jordan-Wigner Transformation
(-0.09886396978427328+0j) [] +
(-0.04532220205777764+0j) [X0 X1 Y2 Y3] +
(0.04532220205777764+0j) [X0 Y1 Y2 X3] +
(0.04532220205777764+0j) [Y0 X1 X2 Y3] +
(-0.04532220205777764+0j) [Y0 Y1 X2 X3] +
(0.1711977489805745+0j) [Z0] +
(0.16862219157249939+0j) [Z0 Z1] +
(0.12054482203290035+0j) [Z0 Z2] +
(0.165867024090678+0j) [Z0 Z3] +
(0.17119774898057452+0j) [Z1] +
(0.165867024090678+0j) [Z1 Z2] +
(0.12054482203290035+0j) [Z1 Z3] +
(-0.22278593024287632+0j) [Z2] +
(0.17434844183963916+0j) [Z2 Z3] +
(-0.22278593024287635+0j) [Z3]

XやらYやらZがそれぞれパウリ行列を表していることが分ります。

SSVQE

SSVQEはSubspace-Search VQE(SSVQE)の略称で励起状態を探索する手法です。

上記の論文では第k準位までの励起状態を一気に計算できる手法としてWeighted SSVQEが提案されています。
量子コンピューター側では以下の手順で計算を行います。

  1. 試行回路(ansatz circuit)$U(\boldsymbol{\theta})$を構成する。
  2. 互いに直交するk+1個の状態
\begin{align}
\{ \ket{\varphi_j} \}^k_{j=0}\\
\braket{\varphi_i}{\varphi_j}=\delta_{ij}
\end{align}

を準備する。
3. $i< j$に対して$w_i > w_j$として次の期待値を計算する。

\begin{align}
L_{\boldsymbol{w}}(\boldsymbol{\theta}) := \sum^k_{j=0}w_j\bra{\varphi_j} U^{\dagger}(\boldsymbol{\theta}) \mathcal{H} U(\boldsymbol{\theta})\ket{\varphi_j}
\end{align}

古典コンピューターで$\boldsymbol{\theta}$の最適化を行います。
基底状態のインデックスは$k=0$に取られています。
最適化を行い収束した後では$U(\boldsymbol{\theta}^*)\ket{\varphi_j}$が第j番目の励起状態の波動関数を表します。
収束後の変分汎関数は

\begin{align}
L_{\boldsymbol{w}}(\boldsymbol{\theta}^*) = \sum^k_{j=0}w_jE_j
\end{align}

になっているはずです。

初期状態${ \ket{\varphi_j} }^k_{j=0}$は互いに直交しており、基本的にユニタリ変換しかこの状態には作用していないのでSSVQEの計算過程でこの直交性は保存されているはずです。
ハミルトニアンがエルミートの場合、基底状態と励起状態、また励起状態間は直交しているので、計算で得られた波動関数は直交している必要があります。
このためSSVQEでは初期状態に直交する状態を用意しているのだと思います。

$i< j$に対して$w_i > w_j$という条件があるので、

\begin{align}
w_0 > w_1 > w_2 > \cdots > w_{k} \\
E_0 < E_1 < E_2 < \cdots < E_k
\end{align}

全てを同じウェイトにすると変分後の波動関数がどのエネルギー準位にいるのかを比較して調べる必要があるので、エネルギー準位の大小関係とは逆の関係になるようなウェイトを付与して変分関数に設定して、直交性を利用することで効率的に計算しているのだと思います。(感想。)

コード自体は下のページで用意されており気軽に実行できます。
6-3. 励起状態の探索手法 (subspace-search variational quantum eigensolver)

以降では上記のページのコードを実行したことを前提に話を進めます。

遷移行列の計算方法

論文に記載されている遷移行列の計算方法について説明します。
i番目のエネルギー準位を表す状態$\ket{E_i}$とj番目のエネルギー準位を表す状態$\ket{E_j}$と表し、演算子Aの期待値$\bra{E_i}A\ket{E_j}$を求めることを考えます。

量子コンピューターで期待値を計算するには$\bra{E_i}A\ket{E_i}$という風に同じ状態で演算子を挟んだ形でしか計算することが出来ません。
ゆえに$\bra{E_i}A\ket{E_j}$を同じ状態だけで表すように分解する必要があります。

$\bra{E_i}A\ket{E_j}$は一般的に複素数であるので、$\bra{E_i}A\ket{E_j}=c=a+ib$とします。
この複素数$c$の実部$a$と虚部$b$がそもそもどのような形になるかのを計算する。
演算子Aはエルミート演算子であると仮定すると、

\begin{align}
c^* = \biggl(\bra{E_i}A\ket{E_j}\biggr)^* = \bra{E_j}A^{\dagger}\ket{E_i} = \bra{E_j}A\ket{E_i}
\end{align}

$c^*=a-ib$なので、

\begin{align}
a &= \frac{c+c^*}{2} = \frac{\bra{E_i}A\ket{E_j} + \bra{E_j}A\ket{E_i}}{2} \\
b &= \frac{c-c^*}{2i} = \frac{\bra{E_i}A\ket{E_j} - \bra{E_j}A\ket{E_i}}{2i}
\end{align}

まずaに関して、次の状態$\ket{+^x_{ij}}$を定義します：

\begin{align}
\ket{E+^x_{ij}} := \frac{\ket{E_i} + \ket{E_j}}{\sqrt{2}}
\end{align}

この$\ket{+^x_{ij}}$を用いるとaは以下のように表せます。

\begin{align}
a &= \bra{E+^x_{ij}} A \ket{E+^x_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr)
\end{align}

実際に定義のを代入して計算すると、

\begin{align}
a &= \bra{E+^x_{ij}} A \ket{E+^x_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr) \\
&= \frac{1}{2} \biggl(\bra{E_i} + \bra{E_j} \biggr) A \biggl(\ket{E_i} + \ket{E_j} \biggr) - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr) \\
&= \frac{\bra{E_i}A\ket{E_j} + \bra{E_j}A\ket{E_i}}{2}
\end{align}

となり一致します。
bに関しても次の状態$\ket{+^y_{ij}}$を定義します：

\begin{align}
\ket{E+^y_{ij}} := \frac{\ket{E_i} - i \ket{E_j}}{\sqrt{2}}
\end{align}

この$\ket{+^y_{ij}}$を用いるとbは以下のように表せます。

\begin{align}
b &= \bra{E+^y_{ij}} A \ket{E+^y_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr)
\end{align}

実際に計算すると、

\begin{align}
b &= \bra{E+^y_{ij}} A \ket{E+^y_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr) \\
&= \frac{1}{2} \biggl(\bra{E_i} +i \bra{E_j} \biggr) A \biggl(\ket{E_i} -i \ket{E_j} \biggr) - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr) \\
&= (-i)\cdot\frac{\bra{E_i}A\ket{E_j} - \bra{E_j}A\ket{E_i}}{2} \\
&= \frac{\bra{E_i}A\ket{E_j} - \bra{E_j}A\ket{E_i}}{2i}
\end{align}

となり、bと一致します。
ただ論文では次のように定義されていました：

\begin{align}
\ket{E+^y_{ij}} := \frac{\ket{E_i} + i \ket{E_j}}{\sqrt{2}}
\end{align}

この定義だと上で計算したbとは符号が逆になりました。
検算として、

\begin{align}
\ket{E_i} = \left(
    \begin{array}{c}
      1 \\
      0
    \end{array}
  \right),\; \ket{E_j} = \left(
    \begin{array}{c}
      0 \\
      1
    \end{array}
  \right),\; Y = \left(
    \begin{array}{cc}
      0 & -i \\
      i & 0 
    \end{array}
  \right)
\end{align}

に対して$\bra{E_i} Y \ket{E_j}$を計算する。
直接計算すると、

\begin{align}
\bra{E_i} Y \ket{E_j} &= \left(
    \begin{array}{cc}
      1 & 0
    \end{array}
  \right) \left(
    \begin{array}{cc}
      0 & -i \\
      i & 0 
    \end{array}
  \right) \left(
    \begin{array}{c}
      0 \\
      1
    \end{array}
  \right)\\
&= \left(
    \begin{array}{cc}
      1 & 0
    \end{array}
  \right) \left(
    \begin{array}{c}
      -i \\
      0
    \end{array}
  \right) = -i
\end{align}

よって、$\bra{E_i} Y \ket{E_j}=-i$になります。
次に上で説明した公式で計算してみます。

\begin{align}
\ket{E+^x_{ij}} &= \frac{1}{\sqrt{2}}\left(
    \begin{array}{c}
      1 \\
      1
    \end{array}
  \right) \\
\ket{E+^y_{ij}} &= \frac{1}{\sqrt{2}}\left(
    \begin{array}{c}
      1 \\
      -i
    \end{array}
  \right)
\end{align}

また、

\begin{align}
\bra{E_i} Y \ket{E_i} = \bra{E_j} Y \ket{E_j} = 0
\end{align}

より、まず実部のaは、

\begin{align}
a &= \bra{E+^x_{ij}} Y \ket{E+^x_{ij}} - \frac{1}{2}\biggl(\bra{E_i} Y \ket{E_i} + \bra{E_j} Y \ket{E_j} \biggr)\\
&= \frac{1}{2} \left(
    \begin{array}{cc}
      1 & 1
    \end{array}
  \right) \left(
    \begin{array}{cc}
      0 & -i \\
      i & 0 
    \end{array}
  \right) \left(
    \begin{array}{c}
      1 \\
      1
    \end{array}
  \right) \\
&= \frac{1}{2}\left(
    \begin{array}{cc}
      1 & 1
    \end{array}
  \right) \left(
    \begin{array}{c}
      -i \\
      i
    \end{array}
  \right) = 0
\end{align}

となり実部aは0になります。
次の虚部bは、

\begin{align}
b &= \bra{E+^y_{ij}} A \ket{E+^y_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr) \\
&= \frac{1}{2} \left(
    \begin{array}{cc}
      1 & i
    \end{array}
  \right) \left(
    \begin{array}{cc}
      0 & -i \\
      i & 0 
    \end{array}
  \right) \left(
    \begin{array}{c}
      1 \\
      -i
    \end{array}
  \right) \\
&= \frac{1}{2}\left(
    \begin{array}{cc}
      1 & i
    \end{array}
  \right) \left(
    \begin{array}{c}
      -1 \\
      i
    \end{array}
  \right) = -1
\end{align}

よって、$c=a+ib=-i$にになるので両者の結果は一致しました。

上での議論では

\begin{align}
\ket{E+^x_{ij}} := \frac{\ket{E_i} + \ket{E_j}}{\sqrt{2}}
\end{align}

としていましたが、$\ket{E_i}=U(\boldsymbol{\theta})\ket{\varphi_i}$であったことを利用して、

\begin{align}
\ket{E+^x_{ij}} &=  \frac{U(\boldsymbol{\theta}^*)(\ket{\varphi_i} + \ket{\varphi_j})}{\sqrt{2}} \\
\ket{E+^y_{ij}} &=  \frac{U(\boldsymbol{\theta}^*)(\ket{\varphi_i} -i \ket{\varphi_j})}{\sqrt{2}}
\end{align}

であり、$U(\boldsymbol{\theta}^*)$がかかっている状態をまとめて次のように定義します：

\begin{align}
\ket{+^x_{ij}} &:=  \frac{\ket{\varphi_i} + \ket{\varphi_j}}{\sqrt{2}} \\
\ket{+^y_{ij}} &:=  \frac{\ket{\varphi_i} -i \ket{\varphi_j}}{\sqrt{2}}
\end{align}

論文では最初からこちらで定義されています。
この状態を用いて改めて実部aは、

\begin{align}
a &= \bra{E+^x_{ij}} A \ket{E+^x_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr)\\
&= \bra{+^x_{ij}}U^{\dagger}(\boldsymbol{\theta}^*) A U(\boldsymbol{\theta}^*)\ket{+^x_{ij}}\\
& - \frac{1}{2}\bra{\varphi_i} U^{\dagger}(\boldsymbol{\theta}^*) A U(\boldsymbol{\theta}^*) \ket{\varphi_i} \\
& - \frac{1}{2}\bra{\varphi_j} U^{\dagger}(\boldsymbol{\theta}^*) A U(\boldsymbol{\theta}^*) \ket{\varphi_j}
\end{align}

となります。
続いて虚部bは、

\begin{align}
b &= \bra{E+^y_{ij}} A \ket{E+^y_{ij}} - \frac{1}{2}\biggl(\bra{E_i} A \ket{E_i} + \bra{E_j} A \ket{E_j} \biggr)\\
&= \bra{+^y_{ij}}U^{\dagger}(\boldsymbol{\theta}^*) A U(\boldsymbol{\theta}^*)\ket{+^y_{ij}}\\
& - \frac{1}{2}\bra{\varphi_i} U^{\dagger}(\boldsymbol{\theta}^*) A U(\boldsymbol{\theta}^*) \ket{\varphi_i} \\
& - \frac{1}{2}\bra{\varphi_j} U^{\dagger}(\boldsymbol{\theta}^*) A U(\boldsymbol{\theta}^*) \ket{\varphi_j}
\end{align}

となります。
$\ket{E+^x_{ij}}$を作るのは難しいんですが、$\ket{+^x_{ij}}$であれば最初に用意した自明な直交する状態の線型結合で表現することが出来ます。

6-3. 励起状態の探索手法 (subspace-search variational quantum eigensolver)

上のページで第1励起状態までを求める時の初期状態として、

\begin{align}
\ket{\varphi_0} &= \ket{0000}\\
\ket{\varphi_1} &= \ket{0001}
\end{align}

を利用しています。
この初期状態に対して上で定義している$\ket{+^x_{ij}}, \ket{+^y_{ij}}$は

\begin{align}
\ket{+^x_{01}} &= \frac{1}{\sqrt{2}}\biggl(\ket{0000} + \ket{0001} \biggr)\\
\ket{+^y_{01}} &= \frac{1}{\sqrt{2}}\biggl(\ket{0000} -i \ket{0001} \biggr)
\end{align}

と表せます。
まず$\ket{+^x_{01}}$の作り方について説明します。
Qulacsの量子ビットのインデックスの定義の仕方に倣って量子ビットを次のように分解します。

\begin{align}
\ket{\varphi_0} &= \ket{0000} = \ket{0}_3 \otimes \ket{0}_2 \otimes \ket{0}_1 \otimes \ket{0}_0 \\
\ket{\varphi_1} &= \ket{0001} = \ket{0}_3 \otimes \ket{0}_2 \otimes \ket{0}_1 \otimes \ket{1}_0
\end{align}

0量子ビット目だけが異なるので、0量子ビットの$\ket{0}_1$にアダマールゲートを作用させると、

\begin{align}
H_0\ket{0}_0 = \frac{1}{\sqrt{2}}\biggl(\ket{0}_0 + \ket{1}_0 \biggr)
\end{align}

になるので、$\ket{\varphi_0}$の0量子ビットにアダマールゲートを作用させると

\begin{align}
H_0\ket{\varphi_0} &= \ket{0}_3 \otimes \ket{0}_2 \otimes \ket{0}_1 \otimes H_0 \ket{0}_0\\
&= \ket{0}_3 \otimes \ket{0}_2 \otimes \ket{0}_1 \otimes \frac{1}{\sqrt{2}}\biggl(\ket{0}_0 + \ket{1}_0 \biggr)\\
&= \frac{1}{\sqrt{2}}\biggl(\ket{0000} + \ket{0001} \biggr) = \ket{+^x_{01}}
\end{align}

が得られます。
次に$\ket{+^y_{01}}$の作り方について説明します。
まず回転ゲート

\begin{align}
R_{x_j}(\theta) &= e^{\frac{-i\theta}{2}X} = \cos{\frac{\theta}{2}} I_j - i \sin{\frac{\theta}{2}} X_j
\end{align}

より、$\theta=\frac{\pi}{2}$を代入すると、

\begin{align}
R_{x_0}(\frac{\pi}{2}) &= \frac{1}{\sqrt{2}} \biggl(I_0 - iX_0 \biggr)
\end{align}

これを$\ket{0}_0$に作用させると、

\begin{align}
R_{x_0}(\frac{\pi}{2})\ket{0}_0 &= \frac{1}{\sqrt{2}} \biggl(I_0 - iX_0 \biggr)\ket{0}_0\\
&= \frac{1}{\sqrt{2}} \biggl(\ket{0}_0 - i\ket{1}_0 \biggr)
\end{align}

になるので、$\ket{\varphi_0}$の0量子ビットに$R_{x_0}(\frac{\pi}{2})$を作用させると、

\begin{align}
R_{x_0}(\frac{\pi}{2})\ket{\varphi_0} &= \ket{0}_3 \otimes \ket{0}_2 \otimes \ket{0}_1 \otimes R_{x_0}(\frac{\pi}{2}) \ket{0}_0\\
&= \ket{0}_3 \otimes \ket{0}_2 \otimes \ket{0}_1 \otimes \frac{1}{\sqrt{2}}\biggl(\ket{0}_0 -i \ket{1}_0 \biggr)\\
&= \frac{1}{\sqrt{2}}\biggl(\ket{0000} -i \ket{0001} \biggr) = \ket{+^y_{01}}
\end{align}

が得られます。
後はこれら量子回路上で実装して各期待値を計算すれば遷移行列要素が求められます。

コード

上記のロジックを実際にQulacsを使って実装してみました。
本家のcallback関数の外側にtheta_history = []、内側にtheta_history.append(theta_list)を追加しています。
最適化後の$\theta$にはtheta_history[-1]を利用しています。

from qulacs.gate import RX
def get_transition_matrix(theta_opt, hamiltonian):
    # 初期状態の準備。
    state_0 = QuantumState(n_qubit) #|0000> を準備
    state_1 = QuantumState(n_qubit); state_1.set_computational_basis(1) #|0001> を準備

    # 遷移行列を計算するための状態を生成する回路をそれぞれ準備する。
    circuit_x = QuantumCircuit(n_qubit)
    circuit_x.add_H_gate(0)
    circuit_y = QuantumCircuit(n_qubit)
    circuit_y.add_gate(RX(0, np.pi/2.0))

    # state_xを生成する。
    state_x = QuantumState(n_qubit)
    circuit_x.update_quantum_state(state_x)

    # state_yを生成する。
    state_y = QuantumState(n_qubit)
    circuit_y.update_quantum_state(state_y)

    # 上記の全てのstateに対して最適化したパラメータを利用して期待値を計算する。
    circuit = he_ansatz_circuit(n_qubit, depth, theta_opt) #量子回路を構成
    # 状態に演算子を作用させる。
    circuit.update_quantum_state(state_0)
    circuit.update_quantum_state(state_1)
    circuit.update_quantum_state(state_x)
    circuit.update_quantum_state(state_y)

    # 観測量を定義する。
    #transition_hamiltonian = Observable(n_qubit)
    #transition_hamiltonian.add_operator(coef, Pauli_string)
    transition_hamiltonian = create_observable_from_openfermion_text(hamiltonian)

    exp_state_0 = transition_hamiltonian.get_expectation_value(state_0)
    exp_state_1 = transition_hamiltonian.get_expectation_value(state_1)
    exp_state_x = transition_hamiltonian.get_expectation_value(state_x)
    exp_state_y = transition_hamiltonian.get_expectation_value(state_y)
    #print(exp_state_0, exp_state_1)
    #print(exp_state_x, exp_state_y)
    Re_part = exp_state_x - 0.5 * (exp_state_0 + exp_state_1)
    Im_part = exp_state_y - 0.5 * (exp_state_0 + exp_state_1)
    result = complex(Re_part, Im_part)
    return result

hamiltonian = '(1.0+0j) [Y0 Z1 X2 Z3] '
transition_matrix = get_transition_matrix(theta_history[-1],  hamiltonian)
print(transition_matrix)
Output
(-0.00160728648163498+0.0001124135671282008j)

本当は厳密対角化して得られた固有状態間の遷移行列要素を計算してSSVQEの手法と一致しているかを確認する必要があるのですが今回はここで力尽きました...。
今後の宿題として更新させて下さい。
(そもそもSSVQEのコードを実行後に得られた波動関数の成分が厳密対角化で得られた波動関数と異なっていたのでコード上のバグ、あるいは私の認識ミスなどがありそうです。)

検算のために、上記のコードの関数の中のprintのコメントアウトを解除して元のjw_hamiltonianで計算すると、

transition_matrix = get_transition_matrix(theta_history[-1],  str(jw_hamiltonian))
Output
-1.1059333188180025 -0.7329846241163034
-0.9194589258302568 -0.9194589097535413
(4.5636896151002304e-08+6.171361166540379e-08j)

1行目の出力が水素エネルギーの基底状態と励起状態のエネルギーでこれは厳密対角化の結果と一致しています。
2行目が$\ket{+^x_{01}}, \ket{+^y_{01}}$を使って計算したハミルトニアンの期待値です。
3行目が基底状態と励起状態の遷移行列要素の値で、解析的には基底状態と励起状態が直交しているので0になってほしい値です。
計算結果は0には近いんですが数値的には0とみなせるかどうか微妙な大きさかなと思います。

以上、長々となりましたが終わりにさせて頂きます。
最後まで読んで頂きありがとうございました。
(バグなどが見つかり次第更新します。)

19
7
1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
19
7