電子系のハミルトニアンの基底状態を量子コンピュータを用いて計算する手法としてVQE(Variational Quantum Eigensolver)があります。
量子コンピュータで扱えるビット数を超えたハミルトニアンの基底状態を計算する手法であるDeep VQEについて以下の論文を読んだので手法についてまとめてみました。
https://arxiv.org/abs/2007.10917
https://qunasys.com/news/2020/8/10/deep-vqe
また論文上でも計算されているFIG3(a)のN=2のモデルを実際にQulacsを用いて実装してみました。
※私の理解が不正確な部分もあると思うので指摘して下さると助かります。
手法の概要
論文中に記載されているFig1に全体の流れが示されています。
- ハミルトニアンを部分系と相互作用に分解する。(a)
- 部分系をVQEで計算して基底状態を求める。(b, 1st VQE)
- 相互作用と基底状態から局所基底を作成する。(local basis)
- 局所基底を用いてハミルトニアンの行列要素を計算する。
- 相互作用を取り入れた有効ハミルトニアンを作成する。(effective Hamiltonian)
- 有効ハミルトニアンに対してVQEを実行する。(2nd VQE)
局所基底の数を量子ビットに対応させる段階で使用する量子ビット数の削減が行えると理解しました。
3で局所基底を作成するのですがその数が$K$個の場合、これを表現するのに必要な量子ビット数は$\lceil\log_2{K}\rceil$個になります。
これより5で作成する有効ハミルトニアンを表現するために必要な量子ビット数は2つの部分系をくっつけるので$2\lceil\log_2{K}\rceil$個になります。
論文中で計算されている4量子ビットのハイゼンベルグモデルが2つ繋がった場合を考えます。
分解せずに計算しようとすると8量子ビット必要になります。
分解して計算すると最初は4量子ビットを使ってVQEを行い、相互作用の数が$K=7$なので次の計算では部分系1つあたり3量子ビット必要になります。
この部分系が2つ合体することになるので全体で必要な量子ビット数は6になります。
トイモデルなので微量ですが$8\to 6$に削減できていることが分かります。
次は具体的に局所基底や有効ハミルトニアンの作成方法についてみて行きます。
手法の詳細
ハミルトニアンの定義
ハミルトニアンを部分系と相互作用に分割します。
$$
\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}}
\mathcal{H} = \sum_i \mathcal{H}_i + \sum_{ij}V_{ij}
$$
$\mathcal{H}_i$が部分系のハミルトニアン、$V_{ij}$が部分系$i$と$j$に働く相互作用を表します。
ここで$V_{ij}$は次の形を仮定しています。
$$
V_{ij} = \sum_{k} v_{k} W^{(i)}_{k} W^{(j)}_{k}
$$
$W^{(i)}_{k}$は部分系$i$に、$W^{(j)}_{k}$は部分系$j$に働く演算子を表しています。
局所基底の作成
まずは$\mathcal{H}_i$の基底状態をVQEを用いて求めます。
VQEで用いた変分量子回路を$U_i(\theta^{(i)})$、得られた基底状態を$\ket{\psi^{(i)}_0}$と表します。
($i$は部分系$i$のインデックスに対応しています。)
$\theta^{(i), *}$が最適化されたパラメータ
$$
\theta^{(i), *} \equiv \text{argmin}_{\theta} \braket{0^n}{U_i(\theta^{(i)})^{\dagger} \mathcal{H}_i U_i(\theta^{(i)})}{0^n}
$$
とすると、$\ket{\psi^{(i)}_0}$は
$$
\ket{\psi^{(i)}_0} = U_i(\theta^{(i),*})\ket{0^n}
$$
と表すことが出来ます。
この$\ket{\psi^{(i)}_0}$と相互作用に含まれる$W^{(i)}_{k}$から局所基底(local basis)を作成します。
($k$は相互作用演算子の種類に対応しています。)
$$
\ket{\psi^{(i)}_k} \equiv W^{(i)}_{k}\ket{\psi^{(i)}_0}
$$
この局所基底を用いて元のハミルトニアンを書き換えた有効ハミルトニアンを計算していきます。
$W^{(i)}_1$だけは恒等演算子と定義されており、相互作用演算子の数が$K$個ある場合、全部で$K+1$個の局所規定が作成されます。
(以降は局所基底の数を$K$で表します。)
この段階では$\brakett{\psi^{(i)}_k}{\psi^{(i)}_l}=\delta_{kl}$とは限りません。
これをグラムシュミットの正規直交化法を用いて正規直交化します。
正規直交化された局所基底を$\ket{\tilde{\psi}^{(i)}_k}$と表しています。
$$
\ket{\tilde{\psi}^{(i)}_k} = \frac{1}{C_k}\biggl(\ket{\psi^{(i)}_k} - \sum_{l<k} \brakett{\tilde{\psi}^{(i)}_l}{\psi^{(i)}_k} \ket{\tilde{\psi}^{(i)}_l}\biggr)
$$
$C_k$は$\brakett{\tilde{\psi}^{(i)}_k}{\tilde{\psi}^{(i)}_k}=1$を満たすための規格化定数になります。
これを具体的に書き出してみると、$k=1$はそのまま
\begin{align}
\ket{\tilde{\psi}^{(i)}_1} &= \frac{1}{C_1} \ket{\psi^{(i)}_1} \\
C_1 &= \sqrt{\brakett{\psi^{(i)}_1}{\psi^{(i)}_1}}
\end{align}
となります。$k=2$では
\begin{align}
\ket{\tilde{\psi}^{(i)}_2} &= \frac{1}{C_2}\biggl(\ket{\psi^{(i)}_2} - \brakett{\tilde{\psi}^{(i)}_1}{\psi^{(i)}_2} \ket{\tilde{\psi}^{(i)}_1}\biggr) \\
&= \frac{1}{C_2}\biggl(\ket{\psi^{(i)}_2} - \frac{1}{C_1^2}\brakett{\psi^{(i)}_1}{\psi^{(i)}_2} \ket{\psi^{(i)}_1}\biggr)
\end{align}
$k>2$以降も同様に計算できます。
これから$\ket{\tilde{\psi}^{(i)}_k}$を求めるには局所基底の内積$\brakett{\psi^{(i)}_1}{\psi^{(i)}_2}$を計算する必要があることが分かります。
内積は
\begin{align}
\brakett{\psi^{(i)}_k}{\psi^{(i)}_l} = \braket{\psi^{(i)}_0}{W^{(i)\dagger}_{k}W^{(i)}_{l}}{\psi^{(i)}_0}
\end{align}
と展開できることから$\ket{\psi^{(i)}_0}$で$W^{(i)\dagger}_{k}W^{(i)}_{l}$の期待値を測定すればいいことが分かります。
($\ket{\psi^{(i)}_0}$で期待値を取る処理は後の計算でも出てきます。)
注意点は、例えば$W^{(i)\dagger}=X^{(i)}_0, W^{(i)}_{l}=Y^{(i)}_0$の場合、$W^{(i)\dagger}_{k}W^{(i)}_{l}=X^{(i)}_0Y^{(i)}_0=iZ^{(i)}_0$とパウリ演算子の関係使ってオブザーバブルに変換してから$Z^{(i)}_0$の期待値を測定する必要がある点です。
(この測定ちに$i$をかけた値が行列要素の値になります。)
後の計算のしやすさのために$\ket{\tilde{\psi}^{(i)}_k}$を$\ket{\psi^{(i)}_{k^{\prime}}}$で展開した時の係数を$P_{kk^{\prime}}$と定義して、以下のように表しておきます。
\begin{align}
\ket{\tilde{\psi}^{(i)}_k} &= \sum^K_{k^{\prime}=1}P^{(i)}_{kk^{\prime}} \ket{\psi^{(i)}_{k^{\prime}}} \\
\end{align}
$\brakett{\psi^{(i)}_k}{\psi^{(i)}_l}$をそれぞれ測定すれば行列$P$を作成することが出来ます。
この正規直交化された局所基底を以降では単に局所基底と呼びます。
有効ハミルトニアンの作成
$\ket{\tilde{\psi}^{(i)}_k}$を用いて$\mathcal{H}_i, V_{ij}$の行列要素を計算して行きます。
部分系の有効ハミルトニアン
まずは$\mathcal{H}_i$の有効ハミルトニアンの行列要素$(\mathcal{H}^{\text{eff}}_i)_{kl}$は、
\begin{align}
(\mathcal{H}^{\text{eff}}_i)_{kl} &= \braket{\tilde{\psi}^{(i)}_k }{\mathcal{H}_i}{\tilde{\psi}^{(i)}_l} \\
&= \sum_{k^{\prime},l^{\prime}} \braket{\psi^{(i)}_{k^{\prime}} }{P^{(i)*}_{kk^{\prime}} \mathcal{H}_i P^{(i)}_{ll^{\prime}}}{ \psi^{(i)}_{l^{\prime}}} \\
&= \sum_{k^{\prime},l^{\prime}} P^{(i)*}_{kk^{\prime}} \braket{\psi^{(i)}_0}{ W^{(i)}_{k^{\prime}}\;^{\dagger}\mathcal{H}_i W^{(i)}_{l^{\prime}}}{ \psi^{(i)}_0} P^{(i)}_{ll^{\prime}} \\
&= \biggl(P^{(i)*} \bar{\mathcal{H}}^{\text{eff}}_i P^{(i)T} \biggr)_{kl}
\end{align}
となります。ここで$\bar{\mathcal{H}}^{\text{eff}}_i$は
\begin{align}
(\bar{\mathcal{H}}^{\text{eff}}_i)_{kl} &= \braket{\psi^{(i)}_k }{\mathcal{H}_i}{ \psi^{(i)}_l} \\
&= \braket{\psi^{(i)}_0}{ W^{(i)}_{k}\;^{\dagger}\mathcal{H}_i W^{(i)}_{l}}{ \psi^{(i)}_0}
\end{align}
で定義される行列要素になります。
$\mathcal{H}^{\text{eff}}_i$の行列要素を作成するには$\ket{\psi^{(i)}_0}$で
\begin{align}
W^{(i)}_{k}\;^{\dagger}\mathcal{H}_i W^{(i)}_{l}
\end{align}
の期待値を測定して$\bar{\mathcal{H}}^{\text{eff}}_i$の行列要素を求める必要があります。
これと局所基底を作成する時に用いた行列$P$を用いることで$\mathcal{H}^{\text{eff}}_i$を作成することが出来ます。
相互作用の有効ハミルトニアン
局所基底で相互作用
\begin{align}
V_{ij} = \sum_{k} v_{k} W^{(i)}_{k} W^{(j)}_{k}
\end{align}
を両側から挟んで行列要素を計算していきます。
$W^{(i)}_{k}, W^{(j)}_{k}$は部分系$i$, $j$にそれぞれに作用するので
\begin{align}
(V^{\text{eff}}_{ij})_{kk^{\prime}ll^{\prime}} &= \bra{\tilde{\psi}^{(i)}_k} \bra{\tilde{\psi}^{(j)}_{k^{\prime}}} V_{ij} \ket{\tilde{\psi}^{(i)}_l} \ket{\tilde{\psi}^{(j)}_{l^{\prime}}} \\
&= \sum_{\nu} v_{\nu} \braket{\tilde{\psi}^{(i)}_k }{W^{(i)}_{\nu}}{ \tilde{\psi}^{(i)}_l} \braket{\tilde{\psi}^{(j)}_{k^{\prime}} }{W^{(j)}_{\nu}}{ \tilde{\psi}^{(j)}_{l^{\prime}}}
\end{align}
となり、$\braket{\tilde{\psi}^{(i)}_k }{W^{(i)}_{\nu}}{ \tilde{\psi}^{(i)}_l}$を求める必要があります。
これを先ほどと同様に展開していくと
\begin{align}
\braket{\tilde{\psi}^{(i)}_{k}}{W^{(i)}_{\nu}}{\tilde{\psi}^{(i)}_l} &= \sum_{k^{\prime},l^{\prime}} P^{(i)*}_{kk^{\prime}} \braket{\psi^{(i)}_0}{ W^{(i)}_{k^{\prime}}\;^{\dagger}W^{(i)}_{\nu} W^{(i)}_{l^{\prime}}}{ \psi^{(i)}_0} P^{(i)}_{ll^{\prime}} \\
&= \biggl(P^{(i)*} \bar{W}^{(i), \text{eff}}_{\nu} P^{(i)T} \biggr)_{kl}
\end{align}
であり、ここで
\begin{align}
(\bar{W}^{(i), \text{eff}}_{\nu})_{kl} &\equiv \braket{\psi^{(i)}_k}{ W^{(i)}_{\nu}}{\psi^{(i)}_l} = \braket{\psi^{(i)}_0 }{ W^{(i)}_{k}\;^{\dagger}W^{(i)}_{\nu} W^{(i)}_{l}}{ \psi^{(i)}_0}
\end{align}
と定義されます。
以上より$V^{\text{eff}}_{ij}$を計算するには$\ket{\psi^{(i)}_0}$で
\begin{align}
W^{(i)}_{k}\;^{\dagger}W^{(i)}_{\nu} W^{(i)}_{l}
\end{align}
の期待値を観測して$(\bar{W}^{(i), \text{eff}}_{\nu})$の行列要素を計算します。
これと局所基底を作成する時に用いた行列$P$を用いることで$V^{\text{eff}}_{ij}$を作成することが出来ます。
ゲートで表す
これまでで求めたのは有効ハミルトニアンの行列要素でした。
量子コンピュータを使わなけれ局所基底で表した2つの部分系を統合して得られる$K^2\times K^2$次元の行列を対角化すれば固有値、固有状態から基底状態が求まります。
量子コンピュータを用いて計算するには行列要素に対応するゲートで表現する必要があると思います。
$(\mathcal{H}^{\text{eff}}_i)_{kl}$の場合、対応するハミルトニアンの演算子は
\begin{align}
(\mathcal{H}^{\text{eff}}_i)_{kl}\ket{\tilde{\psi}^{(i)}_k} \bra{\tilde{\psi}^{(i)}_l}
\end{align}
になります。
この$\ket{\tilde{\psi}^{(i)}_k}\bra{\tilde{\psi}^{(i)}_l}$を後はゲートで表現する必要があります。
その前にまずは局所基底$\ket{\tilde{\psi}^{(i)}_k}$と量子ビットを対応させる必要があります。
例えば$K=7$の場合は以下の対応関係を考えれば良いはずです。
\begin{align}
\ket{\tilde{\psi}^{(i)}_1} &= \ket{000} \\
\ket{\tilde{\psi}^{(i)}_2} &= \ket{001} \\
\ket{\tilde{\psi}^{(i)}_3} &= \ket{010} \\
\ket{\tilde{\psi}^{(i)}_4} &= \ket{011} \\
\ket{\tilde{\psi}^{(i)}_5} &= \ket{100} \\
\ket{\tilde{\psi}^{(i)}_6} &= \ket{101} \\
\ket{\tilde{\psi}^{(i)}_7} &= \ket{110}
\end{align}
この対応関係から例えば、
\begin{align}
\ket{\tilde{\psi}^{(i)}_1}\bra{\tilde{\psi}^{(i)}_2} &= \ket{000}\bra{001}\\
&= \ket{0}\bra{0}_0 \otimes \ket{0}\bra{0}_1 \otimes \ket{0}\bra{1}_2
\end{align}
となるので、後はそれぞれの量子ビットの射影演算子$\ket{0}\bra{0}_0$をゲートを用いて表せれば良いことになります。
1量子ビットに注目すると
\begin{align}
Z &= \ket{0}\bra{0} - \ket{1}\bra{1} \\
I &= \ket{0}\bra{0} + \ket{1}\bra{1}
\end{align}
から$\ket{0}\bra{0}, \ket{1}\bra{1}$が$Z, I$を用いて表現できます。
同様にして
\begin{align}
\ket{0}\bra{0} &= \frac{1}{2}(I + Z) \\
\ket{1}\bra{1} &= \frac{1}{2}(I - Z) \\
\ket{0}\bra{1} &= \frac{1}{2}(I + iY) \\
\ket{1}\bra{0} &= \frac{1}{2}(I - iY)
\end{align}
を用いれば変換することができると思います。
ハミルトニアンはエルミートなのでこれらのゲートと実数係数を用いて展開できると思います。
有効ハミルトニアン
上記で求めた$\mathcal{H}_{i}, \mathcal{H}_{j}, V_{ij}$の有効ハミルトニアンを新たな部分系としてVQEで基底状態を求めます。
\begin{align}
\mathcal{H}^{\text{eff}}_{ij} = \mathcal{H}^{\text{eff}}_i + \mathcal{H}^{\text{eff}}_j + V^{\text{eff}}_{ij}
\end{align}
部分系が2つだけの場合は上記の処理を1回だけ実行すれば良いですが、3つ以上ある場合は分割すうに応じて処理を繰り返す必要があります。
量子ビット数
最初の部分系の数を$N$とし、1つの部分系で必要な量子ビット数が$n$の場合、系全体を分割せずに一度で解こうとすると合計で$M=nN$量子ビットが必要になります。
一方で分割を大なう場合、局所基底の数が$K$とすると必要な量子ビット数は$m=N\lceil \log_2{K}\rceil$となり、$M\to m$に削減されています。
論文では具体的に4量子ビットの部分系を2次元状に繰り返したモデルを用いて必要なビット数について議論されています。
以下の画像に必要な量子ビット数の計算方法をまとめてみました。
計算結果
論文上で計算されている計算結果を紹介します。
モデルa
まずは下図のように4量子ビットで構成される部分系を1次元状に繋げたモデルです。
(実装パートではこれが2つ繋がったモデルを計算してみました。)
部分系は次の反強磁性ハイゼンベルグモデルで構成されています。
\begin{align}
\mathcal{H}_i &= \sum_{E=(\mu, \nu)}(X^{(i)}_{\mu}X^{(i)}_{\nu} + Y^{(i)}_{\mu}Y^{(i)}_{\nu} + Z^{(i)}_{\mu}Z^{(i)}_{\nu}) \\
&= \sum_{E=(\mu, \nu)}\sum_{A=X,y,Z}A^{(i)}_{\mu} A^{(i)}_{\nu}
\end{align}
ここでエッジは$E={(0,1), (1,2), (2,3), (3,0), (0,2)}$になります。
部分系間には次の相互作用が働いています。
\begin{align}
V_{ij} = \sum_{A=X,Y,Z}(A^{(i)}_{0} A^{(j)}_{2} + A^{(i)}_{2} A^{(j)}_{0})
\end{align}
計算結果が下の表にまとめられています。
Localが相互作用を無視した部分系の合計エネルギー、
Effectiveが有効ハミルトニアンを対角化して得られるエネルギー、
ITEが虚時間発展法で得られたエネルギーをそれぞれ表しています。
$N=8$ではサイズが大きくなるのでEffective, ITEのあたいは記載されていません。
$N=2$はDeep VQE、Effective、ITEの値が一致しています。
$N>2$ではITEに対してEffective, Deep VQEの値が少し大きくなりずれていきますが、良い近似を与えていることが分かります。
ITEで直接計算出来ないような系もDeep VQEを用いると近似的なエネルギー値が求められていることが分かります。
モデルb
次のモデルは下図になります。
ハミルトニアンは先ほどと同様の形をしていますが12量子ビットに増えており、エッジが図のように変わっています。
相互作用は0番目と6番目の量子ビットに作用しています。
この部分系の厳密解は$-21.78$で、Deep VQEで計算した結果は$-21.72$になるようです。
回路の深さを変えて測定したエネルギーの値が以下に記載されています。
これから部分系を繋げていくと$N=2$では$-43.8$、$N=4$では$-87.9$が得られているようです。
実装
モデルaで$N=2$の場合の基底エネルギーをQulacsを用いてDeepVQEで計算してみました。
部分系のハミルトニアンは先ほどと同様に以下のモデルを用います。
\begin{align}
\mathcal{H}_i &= \sum_{E=(\mu, \nu)}(X^{(i)}_{\mu}X^{(i)}_{\nu} + Y^{(i)}_{\mu}Y^{(i)}_{\nu} + Z^{(i)}_{\mu}Z^{(i)}_{\nu}) \\
&= \sum_{E=(\mu, \nu)}\sum_{A=X,y,Z}A^{(i)}_{\mu} A^{(i)}_{\nu}
\end{align}
エッジは$E={(0,1), (1,2), (2,3), (3,0), (0,2)}$で、相互作用は以下を用います。
\begin{align}
V_{01} = \sum_{A=X,Y,Z}(A^{(i)}_{0} A^{(j)}_{2} + A^{(i)}_{2} A^{(j)}_{0})
\end{align}
上記のモデルをDeep VQEで計算した時のコードはnotebook形式で以下にUpしているので参照してください。
notebook上で説明も記載しています。
コスト関数の最適化にはBFGSを用いておりコスト関数の勾配を本当は入力する必要があるので今回はサボっています。
(論文上でも勾配を入れなさいと言われていますがサボりました。すいません。)
また、厳密対角化を用いて求めた結果も以下にupしています。
厳密対角化の結果が"-14.46"になりTable1の値と一致しています。
Deep VQEの結果も同様の値になっており一応計算自体はうまく実行できていると思います。
以上になります。
間違い等があればご指摘お願いします。
時間があればコードは改良して他の系(モデルaの$N>2$やモデルb)でも計算してみます。