15. Quantum Krylov subspace method

この章では、発展的な話題として、Quantum Krylov subspace method、とくにLanczos法の量子版について説明する。

Full-CIの章でも述べたように、Krylov部分空間法は、与えられたハミルトニアン\(H\)と初期状態\(\ket{\psi}\)から生成されるKrylov部分空間

\[ \mathcal{K}_m(H, \ket{\psi}) = \text{span}\{\ket{\psi}, H\ket{\psi}, H^2\ket{\psi}, \ldots, H^{m-1}\ket{\psi}\} \]

を用いて、固有値問題を効率的に解く手法である。

量子コンピュータ上での実行に馴染む形として、ゲート\(U(t) = e^{-iHt}\)を用いて、異なる時間ステップでの時間発展演算子を適用し、 量子版のKrylov部分空間法を考えることができる。

具体的には、初期状態\(\ket{\Phi}_0\)から、Hamiltonian \(H\)による時間発展演算子\(e^{-iHt}\)を作用させて

\[ \ket{\Phi}_0, e^{-iHt_1}\ket{\Phi}_0, \ldots, e^{-iHt_{M}}\ket{\Phi}_0 \]

を順次生成し、これらの非直交基底を用いて、問題をsubspaceで解く。 つまり、

(15.1)\[\begin{split} \begin{align} N_{kl} & = \langle \Phi_k \ket{\Phi_l} = \bra{\Phi_0} e^{-i(t_l-t_k)H} \ket{\Phi_0} \\ \tilde{H}_{kl} & = \bra{\Phi_k} H \ket{\Phi_l} = \bra{\Phi_0} H e^{-i(t_l-t_k)H} \ket{\Phi_0} \end{align} \end{split}\]

のもとで、一般化固有値問題

\[ \tilde{H} \ket{\Phi} = E N \ket{\Phi} \]

を解く。QPEでは量子コンピュータ上(この資料ではもちろんsimulatorだが)でHamiltonianの時間発展に加えて、逆量子フーリエ変換が必要であるのに対し、Quantum Krylovでは前者のみを考える。

ライブラリのimportとPairingHamiltonianのクラスを定義しておく↓

Show code cell source Hide code cell source

import numpy as np
import scipy
import itertools
from itertools import combinations
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit, transpile
from qiskit.circuit.library import PauliEvolutionGate, PauliGate
from qiskit.quantum_info import SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.visualization import plot_histogram
import seaborn as sns
from tqdm import tqdm

class PairingHamiltonian:
    def __init__(self, Norb, Nocc, gval, delta_eps=1.0):
        self.Norb = Norb
        self.Nocc = Nocc
        self.delta_eps = delta_eps
        self.gval = gval
        self.basis = self.make_basis()
        self.epsilon = self.eval_epsilon()
        self.Hmat = self.eval_Hmat()

    def make_basis(self):
        self.basis = []
        for occ in combinations(range(self.Norb), self.Nocc):
            self.basis.append(occ)

        return self.basis
    
    def eval_epsilon(self):
        self.epsilon = [ 2 * i * self.delta_eps for i in range(self.Norb) ]
        return self.epsilon
    
    def eval_Hmat(self):
        dim = len(self.basis)
        self.Hmat = np.zeros((dim, dim))
        for bra_idx, bra in enumerate(self.basis):
            for ket_idx, ket in enumerate(self.basis):
                # Hamming distance
                diff = [ i for i in bra if i not in ket ]
                same = [ i for i in bra if i in ket ]
                # for SPE term
                if bra_idx == ket_idx:
                    self.Hmat[bra_idx, ket_idx] += np.sum( [self.epsilon[i] for i in same])
                    self.Hmat[bra_idx, ket_idx] += - self.gval * len(same) 
                # for pairing term
                if len(diff) == 1:
                    self.Hmat[bra_idx, ket_idx] = - self.gval

        return self.Hmat

def tuple_to_bitstring(tup, Norb, rev=True):
    bitint = 0
    for i in tup:
        bitint += 2**i
    if rev:
        bitstring = "|"+format(bitint, f'0{Norb}b')[::-1]+">"
    else:
        bitstring = "|"+format(bitint, f'0{Norb}b')+">"        
    return bitstring

def cG1(circ, c_qubit, i, j, theta):
    theta_4 = theta / 4 
    circ.cx(i,j)
    circ.ry(theta_4, i)
    circ.cx(j,i)
    circ.ry(-theta_4, i)
    circ.cx(c_qubit, i)
    circ.ry(theta_4, i)
    circ.cx(j,i)
    circ.ry(-theta_4, i)
    circ.cx(c_qubit, i)
    circ.cx(i,j)

def cG1_gate(theta):
    circ = QuantumCircuit(2)
    G(circ, 0, 1, theta)
    circ.name = "cG1"   
    circ = circ.to_gate()
    circ = circ.control(1) 
    return circ

def G(circ, i, j, theta):
    theta_2 = theta / 2 
    circ.cx(i,j)
    circ.ry(theta_2, i)
    circ.cx(j,i)
    circ.ry(-theta_2, i)
    circ.cx(j,i)
    circ.cx(i,j)  

def G_gate(theta):
    circ = QuantumCircuit(2)
    G(circ, 0, 1, theta)
    circ.name = "G"    
    return circ.to_gate()

def get_idx_ancilla_in_string(n_qubit, ancilla, Qiskit_ordering):
    idx_ancilla = None
    if ancilla != None:
        if Qiskit_ordering:
            idx_ancilla = n_qubit-1 - ancilla
        else:
            idx_ancilla = ancilla
    return idx_ancilla

def expec_Zstring(res, idx_relevant, Qiskit_ordering=True, target_qubits=[], ancilla_qubit=None):
    exp_val = exp_val_p0 = exp_val_p1 = 0.0
    n_shot = sum(res.values())
    n_qubit = len(list(res.keys())[0])
    idx_ancilla = get_idx_ancilla_in_string(n_qubit, ancilla_qubit, Qiskit_ordering)
    for bitstr, count in res.items():
        if ancilla_qubit != None and target_qubits != []:
            bitstr_target = ''.join([bitstr[k] for k in range(n_qubit) if k != idx_ancilla])
        else:
            bitstr_target = bitstr
        tmp = 1.0
        for idx in idx_relevant:
            if Qiskit_ordering:
                idx = -1 - idx
            bit = int(bitstr_target[idx])            
            tmp *= (1 - 2 * bit)
        exp_val += tmp * count
        
        if ancilla_qubit != None:
            if int(bitstr[idx_ancilla]) == 0:
                exp_val_p0 += tmp * count
            else:
                exp_val_p1 += tmp * count
    exp_val /= n_shot; exp_val_p0 /= n_shot; exp_val_p1 /= n_shot
    return exp_val, exp_val_p0, exp_val_p1


params_exact = np.array([-0.48104276, -1.03976498, -0.98963981, -1.18481738, -0.54832984])

Norb = 4
Nocc = 2
gval = 0.33  

Hamil = PairingHamiltonian(Norb, Nocc, gval)
evals, evecs = np.linalg.eigh(Hamil.Hmat)
evals = np.linalg.eigvalsh(Hamil.Hmat)
Egs_exact = evals[0]
E_HF = Hamil.Hmat[0,0]

SPEs = Hamil.epsilon
pauli_list = [ ]
obs = [ ]
coeffs = [ ]

# I term
coeff = 0.0
op = "I" * Norb
for i in range(Norb):
    coeff += 0.5 * ( SPEs[i] - Hamil.gval ) 
obs += [op]
coeffs += [coeff]
# -Zp term
for i in range(Norb):
    op = "I" * Norb
    op = op[:i] + "Z" + op[i+1:]
    coeff = -0.5 * ( SPEs[i] - Hamil.gval )

    op = op[::-1]
    obs += [op]
    coeffs += [coeff]
# XX+YY term
for i in range(Hamil.Norb):
    for j in range(i+1, Hamil.Norb):
        factor = - Hamil.gval / 2
        op = "I" * Norb
        op = op[:i] + "X" + op[i+1:j] + "X" + op[j+1:]
        op = op[::-1]
        obs += [op]
        coeffs += [ factor ]
        op = "I" * Norb
        op = op[::-1]
        op = op[:i] + "Y" + op[i+1:j] + "Y" + op[j+1:]
        obs += [op]
        coeffs += [ factor ]

hamiltonian_op = SparsePauliOp(obs, coeffs)

# print("basis:", Hamil.basis)
# print([tuple_to_bitstring(tup, Norb) for tup in Hamil.basis])
# print("eps: ", Hamil.epsilon)
# print("Hmat: ", Hamil.Hmat)
print("evals: ", evals)
# print("Egs_exact: ", Egs_exact, " E_HF", E_HF)
# print("gs evec", evecs[:,0])
# print("gs prob", evecs[:,0]**2)

evals:  [1.18985184 3.29649666 5.34       5.34       7.42853393 9.44511758]

15.1. 量子コンピュータ上での\(N_{kl}, \tilde{H}_{kl}\)の計算

式(15.1)の各要素を量子コンピュータでどう測るか、が問題となる。

このような行列要素の計算を含め、線形代数演算を量子コンピュータで行う方法についてはarXiv:1902.10394が詳しい。

最も代表的な方法は、ancilla qubitを用いて(修正0Hadamard testを用いる方法で、以下でもこの方法を説明する。 この他にも、異なるsymmetry sectorに属する状態の重ね合わせを作ることでancillaや制御ユニタリを用いずに行列要素を測定する方法も提案されている。 Cristian L. Cortes and Stephen K. Gray: Phys. Rev. A 105, 022417 (2022)

具体的な説明に入る前に、少し約束事とconventionを整理しておこう。

通常のLanczos法においても、Block Lanczosやmulti-reference Lanczosなどと呼ばれるように、 複数の初期状態からKrylov部分空間を生成する手法があるが、ここではsingle-referenceの場合を考える。

\(\ket{0}^{\otimes N}\)で初期化されたターゲット量子ビットに対して、 reference state作成を行うためのユニタリ演算子を\(U_{p}\)とし、 ハミルトニアンの時間発展演算子を\(U_H(t)\)と書くことにして、 これらを組み合わせたユニタリ演算子を\(U_i \equiv U_H(t_i) U_p\)とする。

15.1.1. overlap \(\langle \Psi_a | \Psi_b \rangle\) の計算

まず

\[\begin{split} \begin{align} \ket{\Psi} & = \frac{1}{\sqrt{2}} \left( \ket{0} \otimes \ket{x} + \ket{1} \otimes \ket{y} \right) \\ \end{align} \end{split}\]

という状態が作れたとして、\(\langle x | y \rangle\)を測るにはどうすればよいかを考える。 初めのqubitに対してHadamardゲートを作用させると、

\[\begin{split} \begin{align} \ket{\Psi} & \to \frac{1}{2} \left( \ket{0} \otimes (\ket{x} + \ket{y}) + \ket{1} \otimes (\ket{x} - \ket{y}) \right) \\ \end{align} \end{split}\]

となる。このとき、ancillaを測定し、\(\ket{0}\)が得られる確率は

\[ \begin{align} P(0) & = \frac{1}{2} \left( \bra{x} + \bra{y} \right) \left( \ket{x} + \ket{y} \right) = \frac{1}{2} \left( 1 + \mathrm{Re} \langle x | y \rangle \right) \end{align} \]

となる。また、ancillaに\(S^\dagger\)-ゲートに相当するphase-rotationを施し

\[\begin{split} \begin{align} \ket{\Psi} & \to \frac{1}{2} \left( \ket{0} \otimes \ket{x} - i \ket{1} \otimes \ket{y} \right) \\ \end{align} \end{split}\]

Hadamardゲートを作用させると、

\[\begin{split} \begin{align} \ket{\Psi} & \to \frac{1}{2} \left( \ket{0} \otimes (\ket{x} - i \ket{y}) + \ket{1} \otimes (\ket{x} + i \ket{y}) \right) \\ \end{align} \end{split}\]

となり、\(\ket{0}\)が得られる確率は

\[ \begin{align} P(0) & = \frac{1}{2} \left( \bra{x} + i \bra{y} \right) \left( \ket{x} - i \ket{y} \right) = \frac{1}{2} \left( 1 + \mathrm{Im} \langle x | y \rangle \right) \end{align} \]

となる。このように、\(\langle x | y \rangle\)の実部と虚部を測定することができる。 このことは、\(\langle \psi | X | \psi \rangle = \langle \psi | H | \psi \rangle, \langle \psi | Y |\psi \rangle = \langle \psi | S X S^\dagger | \psi \rangle = \langle \psi | S H Z H S^\dagger | \psi \rangle\)から、 \(X\)と\(Y\)成分の測定を行うことで、overlap行列要素を測定することができることを意味する。

実際の回路は以下のとおりである。２つの状態の重ね合わせを作り、ancilla qubitを用いてoverlapを測定する。

15.1.1.1. 計算例

ターゲット量子ビットが4つあるとして、適当な初期状態から、これまた適当なオペレータで時間発展させた場合のオーバーラップを計算してみよう。

検算しやすいように、

\[\begin{split} \begin{align} \ket{\Phi_0} & = \frac{\ket{1100} + \ket{0011} }{\sqrt{2}} \\ H & = 0.123 Z_0Z_1 \end{align} \end{split}\]

とでもしておこう。

\(H \ket{\Phi_0} = 0.123 \ket{\Phi_0}\)なので、\(e^{-iHt} \ket{\Phi_0}\)を計算すると、

\[\begin{split} \begin{align} \ket{\Psi_a} & = e^{-iHt_a} \ket{\Phi_0} = e^{- 0.123 i t_a} \ket{\Phi_0} \\ \ket{\Psi_b} & = e^{-iHt_b} \ket{\Phi_0} = e^{- 0.123 i t_b} \ket{\Phi_0} \\ \langle \Psi_a | \Psi_b \rangle & = e^{0.123 i (t_a - t_b)} = \cos(0.123 (t_a - t_b)) + i \sin(0.123 (t_a - t_b)) \end{align} \end{split}\]

となる。

状態作成のための回路を適当に作ってみよう。

ancilla_qubits = [0]
target_qubits = list(range(1,Norb+1))

operator = ["ZZII"]
coeff = [0.123]
hamil_op = SparsePauliOp(operator, coeff)

t_i = 0.12
t_j = 0.56

Re_Exact = np.cos( np.sum(coeff) * (t_i-t_j))
Im_Exact = np.sin( np.sum(coeff) * (t_i-t_j))
print("Exact: ", Re_Exact, Im_Exact)

Exact:  0.9985358702188959 -0.05409358451972206

qc_prep = QuantumCircuit(Norb)
qc_prep.h(0)
qc_prep.cx(0,1)
qc_prep.x(0)
qc_prep.cx(0,2)
qc_prep.cx(2,3)
qc_prep.x(0)
qc_prep.draw('mpl')

state_vector = Statevector.from_instruction(qc_prep)
state_vector.draw("latex")

\[\frac{\sqrt{2}}{2} |0011\rangle+\frac{\sqrt{2}}{2} |1100\rangle\]

これをユニタリ演算子\(U_p\)としておく。

Uprep = qc_prep.to_gate()
Uprep.name = "$U_p$"
Updag = Uprep.inverse()

qc_prep = QuantumCircuit(Norb)
qc_prep.append(Uprep, range(Norb))
state_vector = Statevector.from_instruction(qc_prep)
display(state_vector.draw("latex"))

# 逆演算U^{\dagger}_pを作用させて|0>に戻るかチェック
qc_prep.append(Updag, range(Norb))
state_vector = Statevector.from_instruction(qc_prep)
display(state_vector.draw("latex"))

\[\frac{\sqrt{2}}{2} |0011\rangle+\frac{\sqrt{2}}{2} |1100\rangle\]

\[ |0000\rangle\]

\(U_i \equiv e^{- i H t_i}\)と\(U_j \equiv e^{-i H t_j}\)、これらのcontrolled-Uを用意しておく:

trotter_steps = 1 
Ui = PauliEvolutionGate(hamil_op, t_i, synthesis=SuzukiTrotter(order=1,reps=trotter_steps))
Uj = PauliEvolutionGate(hamil_op, t_j, synthesis=SuzukiTrotter(order=1,reps=trotter_steps))

circuit_Ui = QuantumCircuit(Norb)
circuit_Ui.append(Ui, range(Norb))
circuit_Ui.name = '$U(t_i)$'
gate_Ui = circuit_Ui.to_gate()

circuit_Uj = QuantumCircuit(Norb)
circuit_Uj.append(Uj, range(Norb))
circuit_Uj.name = '$U(t_j)$'
gate_Uj = circuit_Uj.to_gate()

# c-Ui: controlled U = exp(-iHt_i)
circuit_cUi = QuantumCircuit(Norb)
circuit_cUi.append(Uprep, range(Norb))
circuit_cUi.append(Ui, range(Norb))
circuit_cUi.name = 'c-$U(t_i) U_p$'
gate_cUi = circuit_cUi.to_gate().control(1)

# C-Uj 
circuit_cUj = QuantumCircuit(Norb)
circuit_cUj.append(Uprep, range(Norb))
circuit_cUj.append(Uj, range(Norb))
circuit_cUj.name = 'c-$U(t_j) U_p$'
gate_cUj = circuit_cUj.to_gate().control(1)

まず \(\ket{\Psi_i} = U_i U_p |0\rangle^{\otimes 4}, \ket{\Psi_j} = U_j U_p |0\rangle^{\otimes 4}\)を作り、overlapをstatevectorで計算してみよう。

qc_psi_i = QuantumCircuit(Norb)
qc_psi_i.append(Uprep, range(Norb))
qc_psi_i.append(gate_Ui, range(Norb))
psi_i = Statevector.from_instruction(qc_psi_i)

qc_psi_j = QuantumCircuit(Norb)
qc_psi_j.append(Uprep, range(Norb))
qc_psi_j.append(gate_Uj, range(Norb))
psi_j = Statevector.from_instruction(qc_psi_j)

overlap = np.dot(psi_i.data.conj(), psi_j.data)
Re_overlap = np.real(overlap)
Im_overlap = np.imag(overlap)

print("psi_i"); display(psi_i.draw("latex"))
print("psi_j"); display(psi_j.draw("latex"))
print("<psi_i|psi_j>", overlap)
print("cos(t_i-t_j)", np.cos( np.sum(coeff) * (t_i-t_j) ), "sin(t_i-t_j)",  np.sin( np.sum(coeff) * (t_i-t_j) ) )

psi_i

\[(0.7070297583 - 0.0104365171 i) |0011\rangle+(0.7070297583 - 0.0104365171 i) |1100\rangle\]

psi_j

\[(0.7054300263 - 0.0486670107 i) |0011\rangle+(0.7054300263 - 0.0486670107 i) |1100\rangle\]

<psi_i|psi_j> (0.9985358702188957-0.05409358451972205j)
cos(t_i-t_j) 0.9985358702188959 sin(t_i-t_j) -0.05409358451972206

# circuit to measure Re <psi_i|psi_j>
qc_re = QuantumCircuit(1+Norb,1)
qc_re.h(0)
qc_re.append(gate_cUi, ancilla_qubits + target_qubits)
qc_re.x(0)
qc_re.append(gate_cUj, ancilla_qubits+target_qubits)
qc_re.h(0)
qc_re.measure(0,0)
display(qc_re.draw('mpl'))
qc_re = qc_re.decompose()

# circuit to measure Im <psi_i|psi_j>
qc_im = QuantumCircuit(1+Norb,1)
qc_im.h(0)
qc_im.append(gate_cUi, ancilla_qubits + target_qubits)
qc_im.x(0)
qc_im.append(gate_cUj, ancilla_qubits+target_qubits)
qc_im.sdg(0)
qc_im.h(0)
qc_im.measure(0,0)
display(qc_im.draw('mpl'))
qc_im = qc_im.decompose()

from qiskit_aer.primitives import SamplerV2
sampler = SamplerV2()
num_shot = 10**5

job = sampler.run([qc_re, qc_im], shots=num_shot)
results  = job.result()

# Real part
prob = results[0].data.c.get_counts()
p0 = prob.get('0',0) / num_shot
p1 = prob.get('1',0) / num_shot
print(f"Estimated Re<psi_i|psi_j> = {p0 - p1:9.5f}  Exact: {Re_overlap:9.5f} Error: {np.abs(p0-p1-Re_overlap):6.2e}")

# Imaginary part
prob = results[1].data.c.get_counts()
p0 = prob.get('0',0) / num_shot
p1 = prob.get('1',0) / num_shot
print(f"Estimated Im<psi_i|psi_j> = {p0 - p1:9.5f}  Exact: {Im_overlap:9.5f} Error: {np.abs(p0-p1-Im_overlap):6.2e}")

Estimated Re<psi_i|psi_j> =   0.99846  Exact:   0.99854 Error: 7.59e-05
Estimated Im<psi_i|psi_j> =  -0.05182  Exact:  -0.05409 Error: 2.27e-03

15.1.2. \(\langle \Psi_a | H | \Psi_b \rangle \)の計算

次に、\(\langle \Psi_a | H | \Psi_b \rangle\)の計算を考える。

対角成分は、状態作成の回路を用いて、\( \langle \hat{H} \rangle = \sum_n \langle \Psi_i | \hat{H}_n | \Psi_i \rangle\)を測定すれば良い。 ここで、Hamiltonianの各項は、\(\hat{H}_n = G_0 \otimes G_1 \otimes G_2 \otimes \ldots, G \equiv \{I, X, Y, Z\}\)と、各パウリ演算子の積で表現してある。

15.1.3. ハミルトニアンの各項の期待値測定のための回路

"IIXIXII"のようなpauli string(的な変数)に対して、必要なbasis rotationを足す関数を作っておこう。 必要な操作は\(X\)と\(Y\)の測定のための回路の追加と、そのZ基底での測定である。

Show code cell source Hide code cell source

def additional_qc(qc_in, pauli_str, register_target, Qiskit_order=True):
    pauli_str = str(pauli_str)
    if Qiskit_order:
        pauli_str = pauli_str[::-1]

    for i in range(len(pauli_str)):
        if pauli_str[i] == "X":
            qc_in.h(register_target[i])
        elif pauli_str[i] == "Y":
            qc_in.sdg(register_target[i]); qc_in.h(register_target[i])
        elif pauli_str[i] == "Z" or pauli_str[i] == "I":
            pass
        else:
            raise ValueError("Invalid Pauli string: ", pauli_str)

def get_idx_to_measure(pauli_str, Qiskit_order=True):
    idxs = [ idx for idx, p in enumerate(pauli_str) if p != "I"]
    if Qiskit_order:
        idxs = [ len(pauli_str) - 1 - idx for idx in idxs]
    return idxs

hamil = ["XXII", "IYIY"]
for h in hamil:
    idx_relevant = get_idx_to_measure(h)
    print(h) # idx_relevant)
    qc = QuantumCircuit(Norb)
    qc.append(Uprep, range(Norb))
    additional_qc(qc, h, range(Norb))    
    display(qc.draw('mpl', scale=0.3))

XXII

IYIY

既にQiskitの章で作った、pairing Hamiltonianの状態作成回路などを用いて期待値計算が正しく出来ることを確認しておこう。

def ansatz(params, method="FCI"):
    qc = QuantumCircuit(Norb)
    # HF
    for i in range(Nocc):
        qc.x(i)
    # Additional gates
    if method == "FCI":
        if Norb != 4 or Nocc != 2:
            print("You are using parameters to reproduce the FCI wavefunction for Norb=4 and Nocc=2")
            print("Please make sure whether this is correct for your case!")
        qc.append(G_gate(params[0]), [1, 2])
        qc.append(G_gate(params[1]), [2, 3])
        qc.append(cG1_gate(params[2]), [2, 0, 1])
        qc.append(cG1_gate(params[3]), [3, 0, 1])
        qc.append(cG1_gate(params[4]), [3, 1, 2])        
    return qc

Hamil_coeffs = hamiltonian_op.coeffs

idxs_relevant = [ ]
qcs = [ ]
verbose = False
using_statevector = True
for h_term in list(hamiltonian_op.paulis):
    h_term = h_term.to_label()
    idx_relevant = get_idx_to_measure(h_term)
    idxs_relevant.append(idx_relevant)
    qc = ansatz(params_exact)
    additional_qc(qc, h_term, range(Norb))
    if verbose:
        print(h_term, idx_relevant)
        display(qc.draw('mpl', scale=0.5))
    qc = qc.decompose().decompose()
    if not(using_statevector):
        qc.measure_all()
    qcs.append(qc)


if using_statevector:
    results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs]
else:
    job = sampler.run(qcs, shots=num_shot)
    results  = job.result()

Esum = 0.0
for idx in range(len(qcs)):
    h_term = hamiltonian_op.paulis[idx]
    idx_relevant = idxs_relevant[idx]
    if using_statevector:
        res = results[idx]
    else:
        res = results[idx].data.meas.get_counts()

    expval, dummy, dummy_ = expec_Zstring(res, idx_relevant)
    print(h_term.to_label(), Hamil_coeffs[idx], expval)
    Esum += Hamil_coeffs[idx] * expval
print("Esum", Esum, "Egs_exact", Egs_exact)

IIII (5.34+0j) 1.0
IIIZ (0.165+0j) -0.9719921725083537
IIZI (-0.835+0j) -0.9138806897044239
IZII (-1.835+0j) 0.9138801737253329
ZIII (-2.835+0j) 0.9719926884874447
IIXX (-0.165+0j) 0.048212081880989244
YYII (-0.165+0j) 0.04821177091461068
IXIX (-0.165+0j) 0.19420741548699116
YIYI (-0.165+0j) 0.19420491023445183
XIIX (-0.165+0j) 0.13006151287556628
YIIY (-0.165+0j) 0.13006151287556622
IXXI (-0.165+0j) 0.35567940424674865
IYYI (-0.165+0j) 0.35567940424674854
XIXI (-0.165+0j) 0.19420491023445188
IYIY (-0.165+0j) 0.19420741548699086
XXII (-0.165+0j) 0.048211770914610645
IIYY (-0.165+0j) 0.04821208188098913
Esum (1.1898518352304353+0j) Egs_exact 1.1898518351360725

無事、期待値を計算するための回路追加(≒basis rotation)と自作関数が正しく動いたことが確認できた。

15.1.4. 対角成分 \(\langle \Psi_a | H | \Psi_a \rangle\) の計算

まず対角成分の計算をする。

Hamiltonianの固有状態にすればエネルギーの推定値が得られるはずで、チェックになる。 加えて、適当な時間発展した状態の場合もあわせて考えよう。その際は時間発展は適当な値にしておく。

state = "exact"
#state = "Uprep"

Hamil_coeffs = hamiltonian_op.coeffs
Hamil_paulis = hamiltonian_op.paulis
qcs = []
verbose = False
using_statevector = True
for idx_H in range(len(Hamil_paulis)):
    op_string = Hamil_paulis[idx_H].to_label()
    idx_relevant = get_idx_to_measure(op_string)
    if state == "exact":
        qc_exact = ansatz(params_exact)
        qc_reH_D = qc_exact.copy()
    elif state == "Uprep":
        qc_reH_D = QuantumCircuit(Norb)
        qc_reH_D.append(Uprep, range(Norb))    
        qc_reH_D.append(gate_Ui, range(Norb))
    additional_qc(qc_reH_D, op_string, range(Norb))
    if verbose:
        print(op_string, Hamil_coeffs[idx_H])
        display(qc_reH_D.draw('mpl', scale=0.3))
    if not(using_statevector):
        qc_reH_D.measure_all()
    qc_reH_D = qc_reH_D.decompose().decompose()
    qcs.append(qc_reH_D)

if using_statevector:
    results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs]
else:
    job = sampler.run(qcs, shots=num_shot)
    results  = job.result()

Hsum = 0.0
for idx_H in range(len(Hamil_paulis)):
    op_string = Hamil_paulis[idx_H].to_label()
    idx_relevant = get_idx_to_measure(op_string)
    if using_statevector:
        res = results[idx_H]
    else:
        res = results[idx_H].data.meas.get_counts()
    expval, dummy, dummy_ = expec_Zstring(res, idx_relevant)
    Hsum += Hamil_coeffs[idx_H] * expval
    print("operator: ", op_string, "coeff: ", Hamil_coeffs[idx_H], "exp. val: ",  expval)

print("using statevector", using_statevector, "Hsum", Hsum)
if state == "exact":
    print("Hsum should be close to Egs_exact", Egs_exact, "diff.", Hsum - Egs_exact)

operator:  IIII coeff:  (5.34+0j) exp. val:  1.0
operator:  IIIZ coeff:  (0.165+0j) exp. val:  -0.9719921725083537
operator:  IIZI coeff:  (-0.835+0j) exp. val:  -0.9138806897044239
operator:  IZII coeff:  (-1.835+0j) exp. val:  0.9138801737253329
operator:  ZIII coeff:  (-2.835+0j) exp. val:  0.9719926884874447
operator:  IIXX coeff:  (-0.165+0j) exp. val:  0.048212081880989244
operator:  YYII coeff:  (-0.165+0j) exp. val:  0.04821177091461068
operator:  IXIX coeff:  (-0.165+0j) exp. val:  0.19420741548699116
operator:  YIYI coeff:  (-0.165+0j) exp. val:  0.19420491023445183
operator:  XIIX coeff:  (-0.165+0j) exp. val:  0.13006151287556628
operator:  YIIY coeff:  (-0.165+0j) exp. val:  0.13006151287556622
operator:  IXXI coeff:  (-0.165+0j) exp. val:  0.35567940424674865
operator:  IYYI coeff:  (-0.165+0j) exp. val:  0.35567940424674854
operator:  XIXI coeff:  (-0.165+0j) exp. val:  0.19420491023445188
operator:  IYIY coeff:  (-0.165+0j) exp. val:  0.19420741548699086
operator:  XXII coeff:  (-0.165+0j) exp. val:  0.048211770914610645
operator:  IIYY coeff:  (-0.165+0j) exp. val:  0.04821208188098913
using statevector True Hsum (1.1898518352304353+0j)
Hsum should be close to Egs_exact 1.1898518351360725 diff. (9.43627398442004e-11+0j)

非対角成分については、先ほどのoverlapを計算するための回路の後半の制御ユニタリ(状態作成+時間発展)の後に、Hamiltonianの各パウリ演算子をターゲット量子ビットに対しても測ってやれば良い。

時間発展させると検算が大変なので、適当な状態を二個\(\ket{A}, \ket{B}\)作って、\(\langle A | \hat{H} | B \rangle\)を計算してみよう。 \(\ket{\psi_i} = \ket{\mathrm{HF}}\), \(\ket{\psi_j} = \ket{\mathrm{FCI}}\)とすると、 \( H_{ij} = \langle \psi_i | \hat{H} | \psi_j \rangle = E \langle HF | FCI \rangle\)となるので検算ができる。

qc = ansatz(params_exact)
gs_state = Statevector.from_instruction(qc)
display(gs_state.draw("latex"))
print("Egs_exact", Egs_exact)

\[0.9712139094 |0011\rangle+0.1819391345 |0101\rangle+0.0981735082 |0110\rangle+0.0981721942 |1001\rangle+0.0636069907 |1010\rangle+0.0178892924 |1100\rangle\]

Egs_exact 1.1898518351360725

HFとのoverlapは0.9712...., E=1.1898...なので、\(H_{ij} = 1.1556...\)となってほしい。 また、"IIII"の項は単純にoverlapを計算すればよいので、期待値が0.9712...になることもチェックに使える。

# circuit to measure Re <psi_i|H|psi_j> (non-diagonal part)
using_statevector =True
qcs = []
for idx_H in range(len(Hamil_paulis)):
    op_string = Hamil_paulis[idx_H].to_label()
    qc_reH_nonD = QuantumCircuit(1+Norb)
    qc_reH_nonD.h(0)
    # Controlled-Ui
    gate_HF = ansatz(params_exact, method="HF").control(1)
    gate_HF.name = "HF"
    qc_reH_nonD.append(gate_HF, ancilla_qubits + target_qubits)
    qc_reH_nonD.x(0)
    # Controlled-Uj
    gate_FCI = ansatz(params_exact, method="FCI").control(1)
    gate_FCI.name = "FCI"
    qc_reH_nonD.append(gate_FCI, ancilla_qubits+target_qubits)

    # To measure Real part
    qc_reH_nonD.h(0)
    additional_qc(qc_reH_nonD, op_string, target_qubits)
    if not(using_statevector):
        qc_reH_nonD.measure_all()
    #display(qc_reH_nonD.draw('mpl', scale=0.6))
    qc_reH = qc_reH_nonD.decompose().decompose()
    qcs.append(qc_reH)

if using_statevector:
    results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs]
else:
    num_shot = 10**6
    job = sampler.run(qcs, shots=num_shot)
    results  = job.result()

Re_H_ij = 0.0
for idx_H in range(len(Hamil_paulis)):
    op_string = Hamil_paulis[idx_H].to_label()
    idx_relevant = get_idx_to_measure(op_string)
    if using_statevector:
        res = results[idx_H]
    else:
        res = results[idx_H].data.meas.get_counts()   
    dummy_e, p0, p1 = expec_Zstring(res, idx_relevant, target_qubits=range(len(target_qubits)), ancilla_qubit=0)
    expval = p0 - p1
    Re_H_ij += Hamil_coeffs[idx_H] * expval
    print("operator: ", op_string, "coeff: ", Hamil_coeffs[idx_H], "exp. val: ",  expval)
print("Re H_{ij}", Re_H_ij)

operator:  IIII coeff:  (5.34+0j) exp. val:  0.9712139094394295
operator:  IIIZ coeff:  (0.165+0j) exp. val:  -0.9712139094394293
operator:  IIZI coeff:  (-0.835+0j) exp. val:  -0.9712139094394294
operator:  IZII coeff:  (-1.835+0j) exp. val:  0.9712139094394293
operator:  ZIII coeff:  (-2.835+0j) exp. val:  0.9712139094394294
operator:  IIXX coeff:  (-0.165+0j) exp. val:  -6.938893903907228e-18
operator:  YYII coeff:  (-0.165+0j) exp. val:  1.0408340855860843e-17
operator:  IXIX coeff:  (-0.165+0j) exp. val:  0.09817350815322777
operator:  YIYI coeff:  (-0.165+0j) exp. val:  0.09817219419757485
operator:  XIIX coeff:  (-0.165+0j) exp. val:  0.06360699065731032
operator:  YIIY coeff:  (-0.165+0j) exp. val:  0.06360699065731042
operator:  IXXI coeff:  (-0.165+0j) exp. val:  0.18193913447081592
operator:  IYYI coeff:  (-0.165+0j) exp. val:  0.18193913447081603
operator:  XIXI coeff:  (-0.165+0j) exp. val:  0.09817219419757496
operator:  IYIY coeff:  (-0.165+0j) exp. val:  0.09817350815322778
operator:  XXII coeff:  (-0.165+0j) exp. val:  -1.214306433183765e-17
operator:  IIYY coeff:  (-0.165+0j) exp. val:  -1.734723475976807e-18
Re H_{ij} (1.1556023355807892+0j)

時間発展を考えてない部分的チェックだが、大丈夫そうだ。

15.2. Quantum Krylov法の実装

iterationごとに、\(n\)回目の状態\(\ket{\Psi_n} = \exp(-iHn \delta t ) \ket{\Phi_0}\)を作成し、\(n-1\)回目overlap行列の要素を計算する。

検算のために、固有状態\(\ket{\Phi}\)の時間発展を考えよう。

\(H \ket{\Phi} = E \ket{\Phi}\)なので、\(e^{-iHt} \ket{\Phi} = e^{-iEt} \ket{\Phi}\)から、 \(\ket{\psi_0} = \ket{\Phi}, \ket{\psi_1} = e^{-iH\delta t}\ket{\Phi}\)とすると、

\[\begin{split} \begin{align} N_{01} & \equiv \langle \psi_0 | \psi_1 \rangle = \bra{\Phi} e^{-i\hat{H}\delta t} \ket{\Phi} = e^{-iE\delta t} \braket{\Phi | \Phi} = e^{-iE\delta t} \nonumber \\ H_{00} & \equiv \langle \psi_0 | \hat{H} | \psi_0 \rangle = E \nonumber \\ H_{01} & \equiv \langle \psi_0 | \hat{H} | \psi_1 \rangle = \bra{\Phi} \hat{H} e^{-i\hat{H}\delta t} \ket{\Phi} = E e^{-iE\delta t} \nonumber \end{align} \end{split}\]

状態作成を厳密解に取ることでこれらの値は、shot errorの範囲内で再現していることがわかる。

以下では、ランダムな試行関数から始めて固有値の推定ができているかをチェックしよう。

Show code cell source Hide code cell source

# def make_U_and_cU(i, delta_t, trotter_steps, hamiltonian_op, ancilla_qubits, target_qubits, Uprep):
#     Ntar = len(target_qubits)
#     time = i * delta_t
#     circuit_U = QuantumCircuit(Ntar)
#     expiHt = PauliEvolutionGate(hamiltonian_op, time, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
#     circuit_U.append(expiHt, range(Ntar))
#     qc_U = circuit_U.decompose().decompose()
    
#     qc_cU = QuantumCircuit(Ntar)
#     qc_cU.append(Uprep, range(Ntar))
#     qc_cU.append(expiHt, range(Ntar))
#     qc_cU = qc_cU.decompose().decompose()
#     qc_cU = qc_cU.to_gate().control(1)

#     return qc_U, qc_cU

# def make_overlap_qc(Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, using_statevector):
#     qc_re = QuantumCircuit(1+Ntar, 1)
#     qc_re.h(0)
#     qc_re.append(gate_cUi, ancilla_qubits + target_qubits)
#     qc_re.x(0)
#     qc_re.append(gate_cUj, ancilla_qubits+target_qubits)
#     qc_re.h(0)
#     if not using_statevector:
#         qc_re.measure(0,0)
#     qc_re = qc_re.decompose()

#     # Overlap, Im part 
#     qc_im = QuantumCircuit(1+Ntar,1)
#     qc_im.h(0)
#     qc_im.append(gate_cUi, ancilla_qubits + target_qubits)
#     qc_im.x(0)
#     qc_im.append(gate_cUj, ancilla_qubits + target_qubits)
#     qc_im.sdg(0)
#     qc_im.h(0)
#     if not using_statevector:
#         qc_im.measure(0,0)
#     qc_im = qc_im.decompose()

#     return qc_re, qc_im

# def measure_overlap(num_shot, Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, sampler, using_statevector,
#                     do_simulation=True):
#     qc_re, qc_im = make_overlap_qc(Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, using_statevector)
#     if do_simulation:
#         if using_statevector:
#             results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in [qc_re, qc_im]]
#             prob_Re = results[0]
#             prob_Im = results[1] 
#         else:
#             job = sampler.run([qc_re, qc_im], shots=num_shot)
#             results  = job.result()
#             prob_Re = results[0].data.c.get_counts()
#             prob_Im = results[1].data.c.get_counts()

#         p0 = np.sum( [ count for bitstr, count in prob_Re.items() if bitstr[-1] == '0' ] ) / np.sum(list(prob_Re.values()))
#         p1 = np.sum( [ count for bitstr, count in prob_Re.items() if bitstr[-1] == '1' ] ) / np.sum(list(prob_Re.values()))
#         ReN = p0 - p1

#         p0 = np.sum( [ count for bitstr, count in prob_Im.items() if bitstr[-1] == '0' ] ) / np.sum(list(prob_Im.values()))
#         p1 = np.sum( [ count for bitstr, count in prob_Im.items() if bitstr[-1] == '1' ] ) / np.sum(list(prob_Im.values()))
#         ImN = p0 - p1

#         U_ij = ReN + 1j * ImN
#         return U_ij
#     else: #only resource estimation
#         print("qc_re:", dict(qc_re.count_ops()))
#         print("qc_im:", dict(qc_im.count_ops()))
#         return None

# def make_cU(Uprep, Ui, Ntar):
#     circuit_cUi = QuantumCircuit(Ntar)
#     circuit_cUi.append(Uprep, range(Ntar))
#     circuit_cUi.append(Ui, range(Ntar))
#     circuit_cUi = circuit_cUi.decompose().decompose()
#     return circuit_cUi.to_gate().control(1)

# def make_nonD_H_qc(op_string, Ntar, ancilla_qubits, target_qubits, gate_cUi, gate_cUj, qcs_re, qcs_im, using_statevector):
#     # real part
#     qc_reH_nonD = QuantumCircuit(1+Ntar)
#     qc_reH_nonD.h(0)
#     qc_reH_nonD.append(gate_cUi, ancilla_qubits + target_qubits)
#     qc_reH_nonD.x(0)
#     qc_reH_nonD.append(gate_cUj, ancilla_qubits + target_qubits)
#     qc_reH_nonD.h(0)
#     additional_qc(qc_reH_nonD, op_string, target_qubits)
#     if not using_statevector:
#         qc_reH_nonD.measure_all()
#     qc_reH = qc_reH_nonD.decompose().decompose()
#     qcs_re.append(qc_reH)

#     # imaginary part
#     qc_imH_nonD = QuantumCircuit(1+Ntar)
#     qc_imH_nonD.h(0)
#     qc_imH_nonD.append(gate_cUi, ancilla_qubits + target_qubits)
#     qc_imH_nonD.x(0)
#     qc_imH_nonD.append(gate_cUj, ancilla_qubits + target_qubits)
#     qc_imH_nonD.sdg(0)
#     qc_imH_nonD.h(0)
#     additional_qc(qc_imH_nonD, op_string, target_qubits)
#     if not using_statevector:
#         qc_imH_nonD.measure_all()
#     qc_imH = qc_imH_nonD.decompose().decompose()
#     qcs_im.append(qc_imH)

#     return None

# def plot_convergence(ws, evals):
#     cols = sns.color_palette("deep", n_colors=len(evals))
#     fig = plt.figure(figsize=(8, 4))
#     ax = fig.add_subplot(1, 1, 1)
#     for it in range(len(ws)):
#         x = it + 1
#         for n in range(len(ws[it])):
#             y = ws[it][n]
#             ax.scatter(x, y, marker="D", color=cols[n])
#     for n in range(len(evals)):
#         ax.axhline(y=evals[n], color='gray', linestyle='dotted')
#     ax.set_xlabel("iteration")
#     ax.set_ylabel("overlap")
#     ax.set_ylim(0, np.max(evals)*1.05)
#     #plt.show()
#     if using_statevector:
#         fig.savefig("QKrylov_Norb"+str(Norb)+"Nocc"+str(Nocc)+"_statevector.pdf")
#     else:
#         fig.savefig("QKrylov_Norb"+str(Norb)+"Nocc"+str(Nocc)+"_shot.pdf")
#     plt.close()

# def QuantumKrylov(Uprep: QuantumCircuit, # circuit to prepare a reference state
#                  hamiltonian_op: SparsePauliOp, # Hamiltonian operator
#                  sampler,
#                  delta_t=0.01,
#                  max_iterations=10,
#                  trotter_steps=1, 
#                  Trotter_order=1,
#                  num_shot=10**4, 
#                  ancilla_qubits=[], target_qubits=[],
#                  using_statevector=False,
#                  do_simulation=True,
#                  ):
#     Ntar = len(target_qubits)
#     N = np.zeros((max_iterations, max_iterations), dtype=np.complex128)
#     H = np.zeros((max_iterations, max_iterations), dtype=np.complex128)
#     ws = []
#     Unitaries = [ ]
#     if True:
#         Ui = PauliEvolutionGate(hamiltonian_op, delta_t, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
#         gate_cUi = make_cU(Uprep, Ui, Ntar)
#         Uj = PauliEvolutionGate(hamiltonian_op, 2*delta_t, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
#         gate_cUj = make_cU(Uprep, Uj, Ntar)

#         qc_Ui = QuantumCircuit(1+Ntar)
#         qc_Ui.append(gate_cUi, range(Ntar+1))
#         qc_Ui = qc_Ui.decompose(reps=5)
#         print("Circuit for c-UiUp:", dict(qc_Ui.count_ops()))
#         U_ij = measure_overlap(num_shot, Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, sampler, using_statevector, do_simulation)
#         num_of_Hamil_term = len(hamiltonian_op.paulis)
#         print("num of Hamil term: ", num_of_Hamil_term)
#     if not(do_simulation):
#         return None
    
#     for it in tqdm(range(max_iterations)):
#         print("iteration: ", it)
#         ## make controlled U = exp(-iHδt * it)
#         Ui = PauliEvolutionGate(hamiltonian_op, it*delta_t, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
#         gate_cUi = make_cU(Uprep, Ui, Ntar)
#         Unitaries.append(gate_cUi)
#         N[it, it] = 1.0
#         ## evaluate overlap to previous states
#         for j in range(it-1, -1, -1):
#             gate_cUj = Unitaries[j]
#             U_ij = measure_overlap(num_shot, Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, sampler, using_statevector, do_simulation)
#             N[it, j] = U_ij
#             N[j, it] = np.conj(U_ij)
#         ### evaluate H_ii no need ancilla qubit
#         qcs = []  
#         for idx_H in range(len(Hamil_paulis)):
#             op_string = Hamil_paulis[idx_H].to_label()
#             idx_relevant = get_idx_to_measure(op_string)
#             qc_reH_D = QuantumCircuit(Ntar)
#             qc_reH_D.append(Uprep, range(Ntar))
#             qc_reH_D.append(Ui, range(Ntar))
#             additional_qc(qc_reH_D, op_string, range(Ntar))
#             if not(using_statevector):
#                 qc_reH_D.measure_all()
#             qc_reH_D = qc_reH_D.decompose().decompose()
#             qcs.append(qc_reH_D)
#         if using_statevector:
#             results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs]
#         else:
#             job = sampler.run(qcs, shots=num_shot)
#             results  = job.result()

#         Hsum = 0.0
#         for idx_H in range(len(Hamil_paulis)):
#             op_string = Hamil_paulis[idx_H].to_label()
#             idx_relevant = get_idx_to_measure(op_string)
#             if using_statevector:
#                 res = results[idx_H]
#             else:
#                 res = results[idx_H].data.meas.get_counts()
#             expval, dummy, dummy_ = expec_Zstring(res, idx_relevant)
#             Hsum += Hamil_coeffs[idx_H] * expval
#             #print("operator: ", op_string, "coeff: ", Hamil_coeffs[idx_H], "exp. val: ",  expval)
#         #print(f"H[{it},{it}] = {Hsum}")
#         H[it, it] = Hsum

#         ### evaluate H_ij ancilla qubit is needed
#         for j in range(it-1, -1, -1):
#             gate_cUj = Unitaries[j]
#             qcs_re = []
#             qcs_im = []
#             for idx_H in range(len(Hamil_paulis)):
#                 op_string = Hamil_paulis[idx_H].to_label()
#                 make_nonD_H_qc(op_string, Ntar, ancilla_qubits, target_qubits, gate_cUi, gate_cUj, qcs_re, qcs_im, using_statevector)
#             # Re part
#             if using_statevector:
#                 results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs_re]
#             else:
#                 job = sampler.run(qcs_re, shots=num_shot)
#                 results  = job.result()
#             Re_H_ij = 0.0
#             for idx_H in range(len(Hamil_paulis)):
#                 op_string = Hamil_paulis[idx_H].to_label()
#                 idx_relevant = get_idx_to_measure(op_string)
#                 if using_statevector:
#                     res = results[idx_H]
#                 else:
#                     res = results[idx_H].data.meas.get_counts()
#                 dummy_e, p0, p1 = expec_Zstring(res, idx_relevant, target_qubits=range(len(target_qubits)), ancilla_qubit=0)
#                 expval = p0 - p1
#                 Re_H_ij += Hamil_coeffs[idx_H] * expval
#                 #print("operator: ", op_string, "coeff: ", Hamil_coeffs[idx_H], "exp. val: ",  expval)
#             print(f"Re H[{it}, {j}]", Re_H_ij)

#             # Im part
#             if using_statevector:
#                 results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs_im]
#             else:
#                 job = sampler.run(qcs_im, shots=num_shot)
#                 results  = job.result() 
#             Im_H_ij = 0.0
#             for idx_H in range(len(Hamil_paulis)):
#                 op_string = Hamil_paulis[idx_H].to_label()
#                 idx_relevant = get_idx_to_measure(op_string)
#                 if using_statevector:
#                     res = results[idx_H]
#                 else:
#                     res = results[idx_H].data.meas.get_counts()   
#                 dummy_e, p0, p1 = expec_Zstring(res, idx_relevant, target_qubits=range(len(target_qubits)), ancilla_qubit=0)
#                 expval = p0 - p1
#                 Im_H_ij += Hamil_coeffs[idx_H] * expval
#             print("Im H_{ij}", Im_H_ij)

#             H[it, j] = Re_H_ij + 1j * Im_H_ij
#             H[j, it] = Re_H_ij - 1j * Im_H_ij

#         # solve generalized eigenvalue problem
#         Nsub = N[:it+1, :it+1]
#         Hsub = H[:it+1, :it+1]        
#         lam, v = scipy.linalg.eigh(Nsub) 
#         # truncate orthogonal basis with small eigenvalues
#         cols = [ i for i in range(it+1) if lam[i] >= 1.e-9]
#         r = len(cols)
#         Ur = v[:,cols]
#         sq_Sigma_inv = np.diag(lam[cols]**(-0.5))

#         X = Ur @ sq_Sigma_inv @ Ur.conj().T
#         Xdag = X.conj().T
#         tildeH = X @ Hsub @ Xdag
#         w, v = scipy.linalg.eigh(tildeH)
        
#         w_r = w[-r:]
#         ws += [w_r]

#         print("eigs of N: ", lam, "cond", np.linalg.cond(Nsub), "r:", r)
#         print(f"w: {w_r}")
#         print("")
#     return H[:it+1,:it+1], N[:it+1,:it+1], ws

# initial = "exact"
# initial = ""

# if initial == "exact":
#     Uprep = ansatz(params_exact) # exact ground state
# else:
#     np.random.seed(0)
#     params_random = np.random.rand(5)
#     Uprep = ansatz(params_random) # random state

# Hamil_coeffs = hamiltonian_op.coeffs
# Hamil_paulis = hamiltonian_op.paulis
# dt = 1.23

# using_statevector = True
# #using_statevector = False
# Hmat, Nmat, Ens = QuantumKrylov(Uprep, hamiltonian_op, sampler, delta_t=dt, max_iterations=min([10,len(evecs)]), trotter_steps=40, 
#                            num_shot=10**4, 
#                            ancilla_qubits=[0], target_qubits=list(range(1,Norb+1)),
#                            using_statevector=using_statevector)
# print("Eexact", evals)
# plot_convergence(Ens, evals)

# if initial == "exact":
#     print("<Ψ0|Ψ1>", np.exp(- 1j * Egs_exact * dt) , "N[0,1]", Nmat[0,1])
#     print("<Ψ0|H|Ψ0>", Egs_exact, "H[0,0]", Hmat[0,0])
#     print("<Ψ0|H|Ψ1>", np.exp(- 1j * Egs_exact * dt) * Egs_exact, "H[0,1]", Hmat[0,1])
#     print("<Ψ0|Ψ2>", np.exp(- 1j * Egs_exact * dt * 2), "N[0,2]", Nmat[0,2])
#     print("<Ψ0|H|Ψ2>", np.exp(- 1j * Egs_exact * dt * 2) * Egs_exact, "H[0,2]", Hmat[0,2])

# print("Nmat:")
# for i in range(Nmat.shape[0]):
#     print(" ".join([f"{Nmat[i,j]:9.5f}" for j in range(Nmat.shape[1])]))
# print("Hmat:")
# for i in range(Hmat.shape[0]):
#     print(" ".join([f"{Hmat[i,j]:9.5f}" for j in range(Hmat.shape[1])]))

# print("Nmat:")
# for i in range(Nmat.shape[0]):
#     print(" ".join([f"{Nmat[i,j]:9.5f}" for j in range(Nmat.shape[1])]))

# Nhand = Nmat.copy()
# Hhand = Hmat.copy()
# for i in range(1, Nhand.shape[0]-1):
#     for j in range(i+1, Nhand.shape[1]):
#         # Assuming Toeplitz structure
#         Hhand[j,i] = Hhand[i, j] = Hhand[i-1, j-1] 
#         Nhand[j,i] = Nhand[i, j] = Nhand[i-1, j-1]

# print("Nmat:")
# for i in range(Nmat.shape[0]):
#     print(" ".join([f"{Nhand[i,j]:9.5f}" for j in range(Nmat.shape[1])]))

# lam, v = scipy.linalg.eigh(Nhand)
# # truncate orthogonal basis with small eigenvalues
# cols = [ i for i in range(len(lam)) if lam[i] >= 1.e-9]
# r = len(cols)
# Ur = v[:,cols]
# sq_Sigma_inv = np.diag(lam[cols]**(-0.5))

# X = Ur @ sq_Sigma_inv @ Ur.conj().T
# Xdag = X.conj().T
# tildeH = X @ Hhand @ Xdag
# w, v = scipy.linalg.eigh(tildeH)

# print("v", w)
# print("Eexact", evals)

\(N_\mathrm{orb}=4, N_\mathrm{occ}=2\)の場合は次元が6なので、iterationが6になった時点で理想的には厳密解と一致するはずである。 ちなみに、3番目の状態は縮退しているので、緑と赤のシンボルが近づいていくのは誤りではない。

もう少しサイズの大きな\(N_\mathrm{orb}=6, N_\mathrm{occ}=2\)の場合の結果(次元=15)も示しておく。

statevectorを用いた場合は上のように、厳密解へと収束する様子が見て取れるが、 shotで測る場合は、ノイズなしのsimulator(shot errorのみ)でも、\(H,N\)の精度により疑似固有値が得られてしまう。 \(N_\mathrm{orb}=4, N_\mathrm{occ}=2\)で、\(10^4\) shotで計算した結果を示しておく。

今の実装では、overlapを対角化し小さな固有値を捨てることで、一般化固有値問題をより安定的に解いているものの、それでもshot測定の場合は反復に対して計算が不安定になってしまう。実機のノイズがあると、overlapやHが悪条件になり計算がさらに不安定になることが予想されるため、応用の上ではさらなる検討が必要であろう。

15.3. 関連する話題

multi-referenceによるQuantum Krylov法は、J. Chem. Theory Comput. 2020, 16, 2326が参考になる。 overlap行列の計算については、Hadamard test以外にも、上で述べた方法など制御ユニタリを必要としない方法も提案されていて、様々な進展があるようだ。

また、回路等を工夫してIBM実機を用いてKrylov diagonalizationを達成した近年の発展なども非常に興味深い。N. Yoshioka et al., Nature Comm. 16, 5014 (2025)

Hの測定回路を減らしたバージョン

def make_U_and_cU(i, delta_t, trotter_steps, hamiltonian_op, ancilla_qubits, target_qubits, Uprep):
    Ntar = len(target_qubits)
    time = i * delta_t
    circuit_U = QuantumCircuit(Ntar)
    expiHt = PauliEvolutionGate(hamiltonian_op, time, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
    circuit_U.append(expiHt, range(Ntar))
    qc_U = circuit_U.decompose().decompose()
    
    qc_cU = QuantumCircuit(Ntar)
    qc_cU.append(Uprep, range(Ntar))
    qc_cU.append(expiHt, range(Ntar))
    qc_cU = qc_cU.decompose().decompose()
    qc_cU = qc_cU.to_gate().control(1)

    return qc_U, qc_cU

def make_overlap_qc(Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, using_statevector):
    qc_re = QuantumCircuit(1+Ntar, 1)
    qc_re.h(0)
    qc_re.append(gate_cUi, ancilla_qubits + target_qubits)
    qc_re.x(0)
    qc_re.append(gate_cUj, ancilla_qubits+target_qubits)
    qc_re.h(0)
    if not using_statevector:
        qc_re.measure(0,0)
    qc_re = qc_re.decompose()

    # Overlap, Im part 
    qc_im = QuantumCircuit(1+Ntar,1)
    qc_im.h(0)
    qc_im.append(gate_cUi, ancilla_qubits + target_qubits)
    qc_im.x(0)
    qc_im.append(gate_cUj, ancilla_qubits + target_qubits)
    qc_im.sdg(0)
    qc_im.h(0)
    if not using_statevector:
        qc_im.measure(0,0)
    qc_im = qc_im.decompose()

    return qc_re, qc_im

def measure_overlap(num_shot, Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, sampler, using_statevector,
                    do_simulation=True):
    qc_re, qc_im = make_overlap_qc(Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, using_statevector)
    if do_simulation:
        if using_statevector:
            results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in [qc_re, qc_im]]
            prob_Re = results[0]
            prob_Im = results[1] 
        else:
            job = sampler.run([qc_re, qc_im], shots=num_shot)
            results  = job.result()
            prob_Re = results[0].data.c.get_counts()
            prob_Im = results[1].data.c.get_counts()

        p0 = np.sum( [ count for bitstr, count in prob_Re.items() if bitstr[-1] == '0' ] ) / np.sum(list(prob_Re.values()))
        p1 = np.sum( [ count for bitstr, count in prob_Re.items() if bitstr[-1] == '1' ] ) / np.sum(list(prob_Re.values()))
        ReN = p0 - p1

        p0 = np.sum( [ count for bitstr, count in prob_Im.items() if bitstr[-1] == '0' ] ) / np.sum(list(prob_Im.values()))
        p1 = np.sum( [ count for bitstr, count in prob_Im.items() if bitstr[-1] == '1' ] ) / np.sum(list(prob_Im.values()))
        ImN = p0 - p1

        U_ij = ReN + 1j * ImN
        return U_ij
    else: #only resource estimation
        print("qc_re:", dict(qc_re.count_ops()))
        print("qc_im:", dict(qc_im.count_ops()))
        return None

def make_cU(Uprep, Ui, Ntar):
    circuit_cUi = QuantumCircuit(Ntar)
    circuit_cUi.append(Uprep, range(Ntar))
    circuit_cUi.append(Ui, range(Ntar))
    circuit_cUi = circuit_cUi.decompose().decompose()
    return circuit_cUi.to_gate().control(1)

def make_Circ_forNondiagH(term_types, Ntar, ancilla_qubits, target_qubits, 
                   gate_cUi, gate_cUj, qcs_re, qcs_im, using_statevector):
    
    for idx_term in range(len(term_types)):
        term = term_types[idx_term]

        qc = QuantumCircuit(1+Ntar)
        qc.h(0)
        qc.append(gate_cUi, ancilla_qubits + target_qubits)
        qc.x(0)
        qc.append(gate_cUj, ancilla_qubits + target_qubits)
        qc_im = qc.copy() # copy here
        qc.h(0)

        qc_im.sdg(0)
        qc_im.h(0)

        if term == "IZ":
            pass
        elif term == "XX":
            qc.h(target_qubits)
            qc_im.h(target_qubits)
        elif term == "YY":
            qc.sdg(target_qubits)
            qc.h(target_qubits)
            qc_im.sdg(target_qubits)
            qc_im.h(target_qubits)
        else:
            raise ValueError(f"Unsupported term type: {term}")
        if not(using_statevector):
            qc.measure_all()
            qc_im.measure_all()

        qc = qc.decompose(reps=5)
        qcs_re.append(qc)

        qc_im = qc_im.decompose(reps=5)
        qcs_im.append(qc_im)

    return None

def plot_convergence(ws, evals):
    cols = sns.color_palette("deep", n_colors=len(evals))
    fig = plt.figure(figsize=(8, 4))
    ax = fig.add_subplot(1, 1, 1)
    for it in range(len(ws)):
        x = it + 1
        for n in range(len(ws[it])):
            y = ws[it][n]
            ax.scatter(x, y, marker="D", color=cols[n])
    for n in range(len(evals)):
        ax.axhline(y=evals[n], color='gray', linestyle='dotted')
    ax.set_xlabel("iteration")
    ax.set_ylabel("overlap")
    ax.set_ylim(0, np.max(evals)*1.05)
    #plt.show()
    if using_statevector:
        fig.savefig("QKrylov_Norb"+str(Norb)+"Nocc"+str(Nocc)+"_statevector.pdf")
    else:
        fig.savefig("QKrylov_Norb"+str(Norb)+"Nocc"+str(Nocc)+"_shot.pdf")
    plt.close()

def get_idx_circuit(op_string, term_types):
    idx_circuit = None
    if set(op_string) == {'I', 'Z'} or set(op_string) == {'I'}:                            
        for i in range(len(term_types)):
            if set(term_types[i]) == {'I', 'Z'} or set(term_types[i]) == {'I'}:
                idx_circuit = i
                break
        if idx_circuit is None:
            raise ValueError(f"Corresponding circuit for {op_string} not found in term_types: {term_types}")
    elif set(op_string) == {'X', 'I'} or set(op_string) == {'X'}:
        for i in range(len(term_types)):
            if set(term_types[i]) == {'X', 'I'} or set(term_types[i]) == {'X'}:
                idx_circuit = i
                break
        if idx_circuit is None:
            raise ValueError(f"Corresponding circuit for {op_string} not found in term_types: {term_types}")
    elif set(op_string) == {'Y', 'I'} or set(op_string) == {'Y'}:
        for i in range(len(term_types)):
            if set(term_types[i]) == {'Y', 'I'} or set(term_types[i]) == {'Y'}:
                idx_circuit = i
                break
        if idx_circuit is None:
            raise ValueError(f"Corresponding circuit for {op_string} not found in term_types: {term_types}")
    else:
        raise ValueError(f"Unexpected operator string: {op_string}. Supported types are 'IZ', 'XX', 'YY' for now.")
    return idx_circuit            

def QuantumKrylov(Uprep: QuantumCircuit, # circuit to prepare a reference state
                 hamiltonian_op: SparsePauliOp, # Hamiltonian operator
                 sampler,
                 delta_t=0.01,
                 max_iterations=10,
                 trotter_steps=1, 
                 num_shot=10**4, 
                 ancilla_qubits=[], target_qubits=[],
                 using_statevector=False,
                 do_simulation=True,
                 ):
    Ntar = len(target_qubits)
    N = np.zeros((max_iterations, max_iterations), dtype=np.complex128)
    H = np.zeros((max_iterations, max_iterations), dtype=np.complex128)
    ws = []
    Unitaries = [ ]
    if True:
        Ui = PauliEvolutionGate(hamiltonian_op, delta_t, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
        gate_cUi = make_cU(Uprep, Ui, Ntar)
        Uj = PauliEvolutionGate(hamiltonian_op, 2*delta_t, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
        gate_cUj = make_cU(Uprep, Uj, Ntar)

        qc_Ui = QuantumCircuit(1+Ntar)
        qc_Ui.append(gate_cUi, range(Ntar+1))
        qc_Ui = qc_Ui.decompose(reps=5)
        print("Circuit for c-UiUp:", dict(qc_Ui.count_ops()))
        U_ij = measure_overlap(num_shot, Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, sampler, using_statevector, do_simulation)
        num_of_Hamil_term = len(hamiltonian_op.paulis)
        print("num of Hamil term: ", num_of_Hamil_term)
    if not(do_simulation):
        return None
    
    for it in tqdm(range(max_iterations)):
        print("iteration: ", it)
        ## make controlled U = exp(-iHδt * it)
        Ui = PauliEvolutionGate(hamiltonian_op, it*delta_t, synthesis=SuzukiTrotter(order=Trotter_order,reps=trotter_steps))
        gate_cUi = make_cU(Uprep, Ui, Ntar)
        Unitaries.append(gate_cUi)
        N[it, it] = 1.0
        ## evaluate overlap to previous states
        for j in range(it-1, -1, -1):
            gate_cUj = Unitaries[j]
            U_ij = measure_overlap(num_shot, Ntar, gate_cUi, gate_cUj, ancilla_qubits, target_qubits, sampler, using_statevector, do_simulation)
            N[it, j] = U_ij
            N[j, it] = np.conj(U_ij)
            print(f"N[{it}, {j}]", N[it,j])
        ### evaluate H_ii no need ancilla qubit
        # ansatz, ansatz+All-H, ansatz+All-SdagH
        qcs = [ ]
        term_types = { 0:"IZ", 1:"XX", 2:"YY"}  # general case => "XXYZ..."
        for idx_term in range(len(term_types)):
            term = term_types[idx_term]
            qc = QuantumCircuit(Ntar)
            qc.append(Uprep, range(Ntar))
            qc.append(Ui, range(Ntar))
            if term == "IZ":
                pass
            elif term == "XX":
                qc.h(range(Ntar))
            elif term == "YY":
                qc.sdg(range(Ntar))
                qc.h(range(Ntar))
            else:
                raise ValueError(f"Unsupported term type: {term}")
            if not(using_statevector):
                qc.measure_all()
            qcs.append(qc.decompose(reps=3))
                      
        if using_statevector:
            results = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs]
        else:
            job = sampler.run(qcs, shots=num_shot)
            results  = job.result()

        Hsum = 0.0
        for idx_H in range(len(Hamil_paulis)):
            op_string = Hamil_paulis[idx_H].to_label()
            idx_relevant = get_idx_to_measure(op_string)

            idx_circuit = get_idx_circuit(op_string, term_types)
            if using_statevector:
                res = results[idx_circuit]
            else:
                res = results[idx_circuit].data.meas.get_counts()
            expval, dummy, dummy_ = expec_Zstring(res, idx_relevant)
            Hsum += Hamil_coeffs[idx_H] * expval
            #print("operator: ", op_string, "coeff: ", Hamil_coeffs[idx_H], "exp. val: ",  expval)
        print(f"H[{it}, {it}]", Hsum)
        H[it, it] = Hsum

        ### evaluate H_ij ancilla qubit is needed
        for j in range(it-1, -1, -1):
            gate_cUj = Unitaries[j]
            qcs_re = []
            qcs_im = []
            make_Circ_forNondiagH(term_types, Ntar, ancilla_qubits, target_qubits,
                                  gate_cUi, gate_cUj, qcs_re, qcs_im, using_statevector)

            # Re part
            if using_statevector:
                results_Re = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs_re]
                results_Im = [ Statevector.from_instruction(qc).probabilities_dict() for qc in qcs_im]
            else:
                job = sampler.run(qcs_re, shots=num_shot)
                results_Re = job.result()
                job = sampler.run(qcs_im, shots=num_shot)
                results_Im = job.result()

            Re_H_ij = Im_H_ij = 0.0
            for idx_H in range(len(Hamil_paulis)):
                op_string = Hamil_paulis[idx_H].to_label()
                idx_circuit = get_idx_circuit(op_string, term_types)

                idx_relevant = get_idx_to_measure(op_string)
                if using_statevector:
                    res_Re = results_Re[idx_circuit]
                    res_Im = results_Im[idx_circuit]
                else:
                    res_Re = results_Re[idx_circuit].data.meas.get_counts()
                    res_Im = results_Im[idx_circuit].data.meas.get_counts()
                dummy_e, p0, p1 = expec_Zstring(res_Re, idx_relevant, target_qubits=range(len(target_qubits)), ancilla_qubit=0)
                expval = p0 - p1
                Re_H_ij += Hamil_coeffs[idx_H] * expval

                dummy_e, p0, p1 = expec_Zstring(res_Im, idx_relevant, target_qubits=range(len(target_qubits)), ancilla_qubit=0)
                expval = p0 - p1
                Im_H_ij += Hamil_coeffs[idx_H] * expval
                #print("operator: ", op_string, "coeff: ", Hamil_coeffs[idx_H], "exp. val: ",  expval)
            print(f"H[{it},{j}]", Re_H_ij + 1j * Im_H_ij)
            H[it, j] = Re_H_ij + 1j * Im_H_ij
            H[j, it] = Re_H_ij - 1j * Im_H_ij

        # solve generalized eigenvalue problem
        Nsub = N[:it+1, :it+1]
        Hsub = H[:it+1, :it+1]        
        lam, v = scipy.linalg.eigh(Nsub) 
        # truncate orthogonal basis with small eigenvalues
        cols = [ i for i in range(it+1) if lam[i] >= 1.e-9]
        r = len(cols)
        Ur = v[:,cols]
        sq_Sigma_inv = np.diag(lam[cols]**(-0.5))

        X = Ur @ sq_Sigma_inv @ Ur.conj().T
        Xdag = X.conj().T
        tildeH = X @ Hsub @ Xdag
        w, v = scipy.linalg.eigh(tildeH)
        
        w_r = w[-r:]
        ws += [w_r]

        print("eigs of N: ", lam, "cond", np.linalg.cond(Nsub), "r:", r)
        print(f"w: {w_r}")
        print("")
    return H[:it+1,:it+1], N[:it+1,:it+1], ws

using_statevector = True # False
Trotter_order = 4; trortter_steps = 10

Trotter_order = 2; trortter_steps = 40
sampler = None
initial = "exact"
initial = ""

if initial == "exact":
    Uprep = ansatz(params_exact) # exact ground state
else:
    np.random.seed(0)
    params_random = np.random.rand(5)
    Uprep = ansatz(params_random) # random state

Hamil_coeffs = hamiltonian_op.coeffs
Hamil_paulis = hamiltonian_op.paulis
dt = 0.345
dt = 1.234


Hmat, Nmat, Ens = QuantumKrylov(Uprep, hamiltonian_op, sampler, delta_t=dt, 
                                max_iterations=6, 
                                trotter_steps=trortter_steps,
                                ancilla_qubits=[0], target_qubits=list(range(1,Norb+1)),
                                using_statevector=using_statevector)
print("Eexact", evals)

print("Nmat:")
print(Nmat)
print("Hmat:")
print(Hmat)

Circuit for c-UiUp: {'u': 39105, 'cx': 21720}
num of Hamil term:  17

  0%|          | 0/6 [00:00<?, ?it/s]

iteration:  0

 17%|█▋        | 1/6 [00:01<00:05,  1.17s/it]

H[0, 0] (1.611274845675447+0j)
eigs of N:  [1.] cond 1.0 r: 1
w: [1.61127485]

iteration:  1
N[1, 0] (-0.02205787422925065+0.684044848669389j)
H[1, 1] (1.6113755425167746+0j)

 33%|███▎      | 2/6 [01:03<02:27, 36.97s/it]

H[1,0] (-0.27429739568458983+0.5443589163485355j)
eigs of N:  [0.3155996 1.6844004] cond 5.337143606763782 r: 2
w: [1.22663682 3.41188057]

iteration:  2
N[2, 1] (-0.023003707506494442+0.6840866153366303j)
N[2, 0] (-0.840645843683734+0.3245200318034092j)
H[2, 2] (1.611712704125866+0j)
H[2,1] (-0.2756979095112377+0.5443975449090435j)

 50%|█████     | 3/6 [03:07<03:50, 76.82s/it]

H[2,0] (-1.0453061009598457+0.7045661378227202j)
eigs of N:  [0.01404922 0.48310804 2.50284273] cond 178.14813133123457 r: 3
w: [1.21112679 3.3478967  6.7463286 ]

iteration:  3
N[3, 2] (-0.02484021806364589+0.6840705215050461j)
N[3, 1] (-0.8415089439262474+0.32156774843884756j)
N[3, 0] (-0.0912994834005122-0.8343352082690061j)
H[3, 3] (1.6120094699157808+0j)
H[3,2] (-0.27801311488626534+0.5441737547266671j)
H[3,1] (-1.047147748265603+0.7003074541127555j)

 67%|██████▋   | 4/6 [06:17<04:03, 121.71s/it]

H[3,0] (0.1902825478425333-1.073441171094444j)
eigs of N:  [0.00787958 0.03117733 0.63485453 3.32608855] cond 422.1148260406033 r: 4
w: [1.19286222 3.3066403  5.90774424 7.53783249]

iteration:  4
N[4, 3] (-0.027578792561831378+0.6840557525047684j)
N[4, 2] (-0.8425808791189628+0.31703613496159855j)
N[4, 1] (-0.08571409149638098-0.8356925151576265j)
N[4, 0] (0.6188257511702022-0.41245926149483925j)
H[4, 4] (1.6117622126714968+0j)
H[4,3] (-0.2815335748620287+0.5435340914884537j)
H[4,2] (-1.0490655023351396+0.6946118726838197j)
H[4,1] (0.19807674156556204-1.0746523765286686j)

 83%|████████▎ | 5/6 [10:31<02:49, 169.22s/it]

H[4,0] (0.46415548786291233-0.7122906789336276j)
eigs of N:  [2.31162508e-03 1.44724402e-02 4.04391377e-02 8.39601586e-01
 4.10317521e+00] cond 1775.0176051239803 r: 5
w: [1.18985194 3.29649668 5.34       7.42853285 9.4450987 ]

iteration:  5
N[5, 4] (-0.031361364250796775+0.6840856167569126j)
N[5, 3] (-0.8444005067938032+0.3100522030934049j)
N[5, 2] (-0.07711874953831793-0.8374510373659039j)
N[5, 1] (0.6233431403322208-0.4040209834358724j)
N[5, 0] (0.3945672814475031+0.9033814324633849j)
H[5, 5] (1.6117325948397814+0j)
H[5,4] (-0.2847194686864231+0.542466670601386j)
H[5,3] (-1.0516106370384422+0.6866162143131196j)
H[5,2] (0.2078989634019971-1.0751399414715586j)
H[5,1] (0.47140502480682817-0.703395206974318j)

100%|██████████| 6/6 [15:46<00:00, 157.75s/it]

H[5,0] (0.4806128383055382+1.4935442823877352j)
eigs of N:  [1.42705988e-08 2.59549683e-03 1.64094847e-02 4.70779651e-02
 9.72900191e-01 4.96101685e+00] cond 347639003.8197853 r: 6
w: [1.18985184 3.29649666 5.34       5.3400002  7.42853393 9.44511758]

Eexact [1.18985184 3.29649666 5.34       5.34       7.42853393 9.44511758]
Nmat:
[[ 1.        +0.j         -0.02205787-0.68404485j -0.84064584-0.32452003j
  -0.09129948+0.83433521j  0.61882575+0.41245926j  0.39456728-0.90338143j]
 [-0.02205787+0.68404485j  1.        +0.j         -0.02300371-0.68408662j
  -0.84150894-0.32156775j -0.08571409+0.83569252j  0.62334314+0.40402098j]
 [-0.84064584+0.32452003j -0.02300371+0.68408662j  1.        +0.j
  -0.02484022-0.68407052j -0.84258088-0.31703613j -0.07711875+0.83745104j]
 [-0.09129948-0.83433521j -0.84150894+0.32156775j -0.02484022+0.68407052j
   1.        +0.j         -0.02757879-0.68405575j -0.84440051-0.3100522j ]
 [ 0.61882575-0.41245926j -0.08571409-0.83569252j -0.84258088+0.31703613j
  -0.02757879+0.68405575j  1.        +0.j         -0.03136136-0.68408562j]
 [ 0.39456728+0.90338143j  0.62334314-0.40402098j -0.07711875-0.83745104j
  -0.84440051+0.3100522j  -0.03136136+0.68408562j  1.        +0.j        ]]
Hmat:
[[ 1.61127485+0.j         -0.2742974 -0.54435892j -1.0453061 -0.70456614j
   0.19028255+1.07344117j  0.46415549+0.71229068j  0.48061284-1.49354428j]
 [-0.2742974 +0.54435892j  1.61137554+0.j         -0.27569791-0.54439754j
  -1.04714775-0.70030745j  0.19807674+1.07465238j  0.47140502+0.70339521j]
 [-1.0453061 +0.70456614j -0.27569791+0.54439754j  1.6117127 +0.j
  -0.27801311-0.54417375j -1.0490655 -0.69461187j  0.20789896+1.07513994j]
 [ 0.19028255-1.07344117j -1.04714775+0.70030745j -0.27801311+0.54417375j
   1.61200947+0.j         -0.28153357-0.54353409j -1.05161064-0.68661621j]
 [ 0.46415549-0.71229068j  0.19807674-1.07465238j -1.0490655 +0.69461187j
  -0.28153357+0.54353409j  1.61176221+0.j         -0.28471947-0.54246667j]
 [ 0.48061284+1.49354428j  0.47140502-0.70339521j  0.20789896-1.07513994j
  -1.05161064+0.68661621j -0.28471947+0.54246667j  1.61173259+0.j        ]]