反復的量子位相推定（IQPE）で水素分子の基底エネルギーを計算する

論文[1]をもとに反復的量子位相推定（IQPE）で水素分子の基底エネルギーを計算することを目標とします。
Quantum Native Dojoの7-1を参考にコードを実装し、実際の挙動を確かめてみました。
あくまで自分の勉強用のメモとして書いたので間違っている箇所などがあれば教えていただけると幸いです。

##反復的量子位相推定について

反復的量子位相推定では補助ビット１つで任意の桁数までの位相推定を行うことができます。 １回の計算過程で固有値の位相$0.j_1...j_t$の１つの桁($j_k$)を求めることができます。詳しくはQuantum Native Dojoの7-1を参考にしてください。

##求める手順
0.ハミルトニアンのサイズを対照性などを利用して削減する。

1. ハミルトニアンの時間発展演算子$U=e^{-iH\tau}$を精度よく近似する。
2. 制御時間発展演算子を量子コンピュータで計算できるサイズに分解し実装する。
3. 基底状態と十分重なりのある初期状態を準備する。
4. IQPEでエネルギー固有値を求める。

以下ではこの手順に沿って説明していきます。

##0. ハミルトニアンのサイズを対照性などを利用して削減する。

水素分子のハミルトニアンを第２量子化し、STO-6G基底関数のもとでBravyi-Kitaev変換を行うと、ハミルトニアンは４量子ビットを用いて表されます。ただ、この量子ビット数は電子数一定という条件を使って２個まで削減することができ、ハミルトニアンは以下の形で表現できます。 $$H=g_0\boldsymbol{1}+g_1Z_0+g_2Z_1+g_3Z_0Z_1+g_4X_0X_1+g_5Y_0Y_1$$
##1. ハミルトニアンの時間発展演算子を精度よく近似する。

IQPEで用いる制御ユニタリ演算$\Lambda(U^{2^{k}})$を実装するため、まずは時間発展演算子$U=e^{-iH\tau}$を量子回路に実装します。まず$g_0\boldsymbol{1}$と$g_3Z_0Z_1$についてはハミルトニアンのほかの項と可換であるので
$$U=\exp \left[-i\tau\sum_{i}g_iH_i \right]=\exp \left[-i\tau g_0\textbf{1} \right] \exp[-i\tau g_3Z_0Z_1]{\rm{exp}}[-i\tau H_{\rm{eff}}]$$
となります。（可換な演算子AとBについてはexp(A+B)=exp(A)exp(B)となることを用いました。）
ここで
$$H_{\rm{eff}}=g_1Z_0+g_2Z_1+g_4X_0X_1+g_5Y_0Y_1$$
です。$g_0\boldsymbol{1}$と$g_3Z_0Z_1$の寄与は後から足し合わせればよいので以下では$H_{\rm{eff}}$の固有値を$U_{\rm{eff}}:={\rm{exp}}[-i\tau H_{\rm{eff}}]$のIQPEを用いて求めることを考えます。

$U_{\rm{eff}}$をトロッター分解すると
$$U_{\text{eff}} = \exp \left[−i \tau \sum_{i=1,2,4,5} g_i H_i \right] \approx U_{\text{Trot}}^{(N)} (\tau) := \left( \prod_{i=1,2,4,5} \exp[-i g_i H_i \tau/N] \right)^N$$
となります。ここで$ U_{\text{Trot}}^{(N)} (\tau)$に表れる積の各項はパウリ行列の積の指数関数の形$\rm{exp}( i\theta P)$の形をしているので量子ゲートによって容易に実装できます。

今回は扱う行列が4×4と小さいので、実際に$H_{\textrm{eff}}$を対角化してみてその最小固有値$E_{\textrm{min}}$を求め、$ U_{\text{Trot}}^{(N)} (\tau)$の固有値の$e^{i \lambda_{\textrm{Trot}}\tau}$の$\lambda_{\textrm{Trot}}$と比較してみます。

まずは$H_{\textrm{eff}}$の対角化を行います。
以下では量子回路シミュレータであるQulacsを用いて計算を行っていきます。
ここでエネルギーの単位はハートリー(Hartree)が使われています。

pip install qulacs pyscf openfermion openfermionpyscf
from functools import reduce
import numpy as np
from numpy.linalg import matrix_power, eig
from scipy.sparse.linalg import eigsh
from openfermion.ops import QubitOperator
from openfermion.linalg import get_sparse_operator
from qulacs import QuantumState, Observable, QuantumCircuit
import matplotlib.pyplot as plt

def hamiltonian_eff():
    """
    distance = 0.70 A
    removed I and 'Z0 Z1' terms, which add up to -1.31916027
    """
    n_qubits = 2
    g_list = [0.3593, 0.0896, -0.4826, 0.0896]## taken from table 1 of paper [1]
    pauli_strings = ['Z0','Y0 Y1','Z1','X0 X1']
    hamiltonian = QubitOperator()
    for g,h in zip(g_list,pauli_strings):
        hamiltonian += g *QubitOperator(h)
    sparse_matrix = get_sparse_operator(hamiltonian, n_qubits=n_qubits)#get_sparse_operatorは演算子から疎行列を得る
    vals, vecs = eigsh(sparse_matrix, k=1, which="SA") ## 最も小さい固有値と対応する固有ベクトルが得られる
    return sparse_matrix, vals

gomi,eigs = hamiltonian_eff()
exact_eigenvalue = eigs[0]
print('exact_eigenvalue:{} Ha'.format(exact_eigenvalue))

次に$ U_{\text{Trot}}^{(N)} (\tau)$を対角化します。後のステップでは$ U_{\text{Trot}}^{(N)} (\tau)$を量子回路として具体的に実装しますが、ここでは$H_i^2=I$のとき
$$
\left( \prod_{i=1,2,4,5} \exp[-i g_i H_i \tau/N] \right)^N = \left( \prod_{i=1,2,4,5} \left( \cos(g_i\tau/N) I -i \sin(g_i\tau/N) H_i \right) \right)^N
$$
となる性質を用いて計算します。そして$N=1,3,...,9$において$ U_{\text{Trot}}^{(N)} (\tau)$の固有値$e^{i \lambda_{\textrm{Trot}};\tau}$の$\lambda_{\textrm{Trot}}$と$E_{\textrm{min}}$を比較します。

def order_n_trotter_approx(t,n_trotter_steps):
    """
    ordering: 'Z0', 'Y0 Y1', 'Z1', 'X0 X1'
    Returns:
        sparse_matrix: trotterized [exp(iHt/n)]^n
        args: list of phases of each eigenvalue, exp(i*phase)
    """
    n_qubits = 2
    g_list = [0.3593, 0.0896, -0.4826, 0.0896]
    pauli_strings = ['Z0','Y0 Y1','Z1','X0 X1']
    terms = []
    for g,h in zip(g_list,pauli_strings):
        arg = g * t / n_trotter_steps
        qop = complex(np.cos(arg),0)* QubitOperator("") - complex(0,np.sin(arg))*QubitOperator(h)
        terms += [get_sparse_operator(qop, n_qubits=n_qubits)]
    
    sparse_matrix = terms[0]
    for i in range(3):
        sparse_matrix = np.dot(sparse_matrix,terms[i+1])
    """
    sparse_matrix = reduce(np.dot, terms)でも同じ
    """
    matrix = matrix_power(sparse_matrix.toarray(),n_trotter_steps)##U_{trot}^{(N)}
    vals, vecs = eig(matrix)## e^{i lambda_{trot} }
    args = np.angle(vals)## [-pi, pi] の範囲で固有値(今回はλ×τに当たる)が得られる->今回、求める固有値が-0.86付近と分かっているので修正不要。
    return sparse_matrix, sorted(args)

tau = 0.640 ## taken from table 1 of paper [1]
print('N, E_trot, |exact_eig - E_trot|')
for n in range(1,10,2):
    gomi,phases = order_n_trotter_approx(tau,n)
    e_trotter = phases[0]/tau ##上で得られたargsはλ×τなのでλを求めるためにτで割る。
    print( f"{n}, {e_trotter:.10f}, {abs(exact_eigenvalue - e_trotter):.3e}" )

実行結果は以下のようになります。

N, E_trot, |exact_eig - E_trot|
1, -0.8602760326, 4.842e-04
3, -0.8607068561, 5.342e-05
5, -0.8607410548, 1.922e-05
7, -0.8607504700, 9.804e-06
9, -0.8607543437, 5.931e-06

このように次数が大きくなるほど近似精度が上がっており、真のエネルギー固有値を化学的精度（ 1.6×10−3 Ha）と呼ばれる、化学計算において必要とされる精度で近似するには今回はN＝１で十分であることが分かります。

##2. 制御時間発展演算子を量子コンピュータで計算できるサイズに分解し実装する。
量子コンピュータ上で制御時間発展演算子 $\Lambda \left( \left( U_{\textrm{Trot}}^{(N)} \right)^{2^k} \right)$ を実装するためには、これを簡単な量子ゲートに分解する必要があります。
今回の例では$ U_{\text{Trot}}^{(N)} $に含まれる、

• $\Lambda(R_Z(\theta))$
• $\Lambda(R_{XX}(\theta))$
• $\Lambda(R_{YY}(\theta))$

という制御回転ゲートを実装できれば良いです。

以下のコードでは、 Qulacs で制御時間発展演算子 $\Lambda \left( \left( U_{\textrm{Trot}}^{(N)} \right)^{2^k} \right)$ の量子回路を実装し、IQPEで実行すべき回路を作っています。

def IQPE_circuit(g_list, tau, kickback_phase, k, n_trotter_step = 1):
    n_qubits = 3 ## 2 for system, and 1 for ancillary
    a_ind = 2 ##index of ancillary
    phi = -(tau/n_trotter_step) * g_list ##パウリ演算子の係数
    circuit = QuantumCircuit(n_qubits)
    ##ここから回路を構成する。まずancilaにアダマールをかける
    circuit.add_H_gate(a_ind)
    ##ancilaにkickback位相回転ゲートをかける
    circuit.add_RZ_gate(a_ind,kickback_phase)
    ##制御時間発展演算子を2^{k-1}回かける
    for _ in range(2**(k-1)):
        for _ in range(n_trotter_step):
            # CU(Z0) つまり controlled exp(i phi[0]*Z_0)
            circuit.add_RZ_gate(0,phi[0])
            circuit.add_CNOT_gate(a_ind,0)
            circuit.add_RZ_gate(0,-phi[0])
            circuit.add_CNOT_gate(a_ind,0)
            
            # CU(Y0 Y1)
            circuit.add_Sdag_gate(0)
            circuit.add_Sdag_gate(1)
            circuit.add_H_gate(0)
            circuit.add_H_gate(1)
            circuit.add_CNOT_gate(0,1)
            circuit.add_RZ_gate(1,phi[1])
            circuit.add_CNOT_gate(a_ind,1)
            circuit.add_RZ_gate(1,-phi[1])
            circuit.add_CNOT_gate(a_ind,1)
            circuit.add_CNOT_gate(0,1)
            circuit.add_H_gate(0)
            circuit.add_H_gate(1)
            circuit.add_S_gate(0)
            circuit.add_S_gate(1)
            
            # CU(Z1)
            circuit.add_RZ_gate(1, phi[2])
            circuit.add_CNOT_gate(a_ind, 1)
            circuit.add_RZ_gate(1, -phi[2])
            circuit.add_CNOT_gate(a_ind, 1)
            
           # CU(X0 X1)
            circuit.add_H_gate(0)
            circuit.add_H_gate(1)
            circuit.add_CNOT_gate(0, 1)
            circuit.add_RZ_gate(1, phi[3])
            circuit.add_CNOT_gate(a_ind, 1)
            circuit.add_RZ_gate(1, -phi[3])
            circuit.add_CNOT_gate(a_ind, 1)
            circuit.add_CNOT_gate(0, 1)
            circuit.add_H_gate(0)
            circuit.add_H_gate(1)

    ## Apply Hadamard to ancilla
    circuit.add_H_gate(a_ind)
    return circuit

##3. 基底状態と十分重なりのある初期状態を準備する。
これまでは簡単のために 𝑈 に作用する状態は固有状態と仮定してきました。しかし、入力波動関数が固有状態と十分近い場合、高い精度でその固有値を求めることができます。

今回の水素分子の基底エネルギーを求める問題の場合、Hatree-Fock状態 |𝜙⟩=|01⟩ が十分基底状態に近いので、これを用いて基底エネルギーを求めます。

##4. IQPEでエネルギー固有値を求める。

from qulacs.circuit import QuantumCircuitOptimizer

def iterative_phase_estimation(g_list, tau, n_itter, init_state, n_trotter_step=1, kickback_phase = 0.0):
    for k in reversed(range(1, n_itter+1)):
        psi = init_state.copy()
        circuit = IQPE_circuit(np.array(g_list), tau, kickback_phase, k, n_trotter_step = n_trotter_step)
        ##実行時間短縮のために回路の最適化を行う
        opt = QuantumCircuitOptimizer()
        max_block_size = 4
        opt.optimize(circuit,max_block_size)
        ##回路の実行
        circuit.update_quantum_state(psi)
        ##部分トレース
        p0 = psi.get_marginal_probability([2,2,0])
        p1 = psi.get_marginal_probability([2,2,1])
        ##位相キックバックの更新
        print(f"k={k:2d}, p0={p0:.3f}, p1={p1:.3f}")
        kth_digit = 1 if (p0 < p1) else 0
        kickback_phase = 0.5*kickback_phase + np.pi*0.5*kth_digit
    return 2 * kickback_phase

```python
n_qubits = 3 #2つは電子配置、１つは補助
g_list = [0.3593, 0.0896, -0.4826, 0.0896]
# pauli_strings = ['Z 0', 'Y 0 Y 1', 'Z 1', 'X 0 X 1']
hf_state = QuantumState(n_qubits)
hf_state.set_computational_basis(0b001) #|0>|01>
tau = 0.640

e_trotter = -0.8602760325707504 ##U_{Trot}^{(N)}の理論値
print(f"e_trotter={e_trotter:.10f}")

result_list = []
for n_itter in range(1, 13):
    iqpe_phase = iterative_phase_estimation(g_list, tau, n_itter, hf_state, n_trotter_step=5 , kickback_phase = 0.0)
    e_iqpe = -iqpe_phase/tau## U=exp(-iH*tau)なのでIQPEでは-H*tauの固有値を求めたことになる。なのでHの固有値に変換する。
    print(f"n_itter={n_itter:2d}, e_iqpe={e_iqpe:10f}, error={np.abs(e_iqpe-e_trotter):.5e}")
    result_list.append([n_itter,e_iqpe])

#結果のプロット
result_array = np.array(result_list)
plt.xlabel("# of digit", fontsize = 15)
plt.ylabel("Error", fontsize = 15)
plt.semilogy(result_array[:,0], np.abs(result_array[:,1]-e_trotter), "bo-")
plt.xlim(0,13)
plt.fill_between([0,13], 1.6e-3, color = "lightgrey") ## chemical accuracy領域を着色 

上の実行結果から12桁以上求めれば化学的精度が達成されることがわかります。

##色々試して気付いたこと

• トロッター分解の次数を増やしても精度は変わらなかった。これはQPEで固有値の位相の各桁を求める際に多数決を利用しているので、１次のままでもほぼ正確に位相の値が求められているからだと考えられる。

• n_itter(求める位相の桁数)は13くらいまでは1分程度で求められた。逆にこれ以上の桁数は求めるのにもっと時間がかかってしまう（指数的な増加）。

• なぜかn_trotter_stepを増加させると一番下の桁の精度が落ちる。

##次に試したいこと

• トロッター分解以外で時間発展演算子を実装する手法を試す。
• juliaなどほかの言語で実装して計算時間を比較する。

##参考文献
[1]P. J. J. O’Malley et al. , “Scalable Quantum Simulation of Molecular Energies“, PHYSICAL REVIEW X 6, 031007 (2016)