격자 해밀토니안의 Krylov 양자 대각화

사용 시간 예상: Heron 또는 Nighthawk 프로세서에서 70분 (참고: 이는 예상치이며, 실제 실행 시간은 다를 수 있습니다.)

학습 목표

• 스펙트럼 필터 역할을 하는 유한 해밀토니안 함수를 학습하는 것으로 Krylov 양자 대각화(KQD)를 해석하는 방법

• 확장 스왑 테스트 측정을 통해 투영된 해밀토니안 및 겹침 행렬을 구축하는 방법

• 결과로 나온 일반화 고유값 문제(GEVP)를 풀고 격자 해밀토니안의 바닥 상태 에너지 추정치를 구하는 방법

사전 준비 사항

• Qiskit primitives

• VQE를 사용한 Heisenberg 체인의 바닥 상태 에너지 추정

• 양자 대각화 알고리즘 과정

배경

이 튜토리얼은 Qiskit 패턴의 맥락에서 Krylov 양자 대각화(KQD) 알고리즘을 구현하는 방법을 설명합니다. 먼저 알고리즘의 이론적 배경을 학습한 후, QPU에서의 실행 데모를 확인할 수 있습니다.

다체 해밀토니안의 저에너지 특성을 추정하는 것은 양자 시뮬레이션에서 핵심적인 작업입니다. 예를 들어, 바닥 상태 에너지와 저준위 여기는 화학적 안정성, 자기 정렬, 양자 상전이, 물질 반응과 직접적으로 관련이 있습니다. 고전 컴퓨터에서는 힐베르트 공간의 차원이 궤도나 스핀의 수에 따라 지수적으로 증가하므로, 직접 대각화는 곧 실용적이지 않게 됩니다.

이 문제에 대한 여러 양자 컴퓨팅 접근 방식이 있습니다. 변분 양자 고유값 해법(VQE)과 같은 근시일 변분 방법은 비교적 얕은 매개변수화된 Circuit을 사용하지만, 많은 양자 Circuit 평가가 필요한 비선형 고전 최적화 루프가 필요합니다. 반면 양자 위상 추정(QPE)은 엄격한 보장과 함께 고유값 추정에 대한 더 직접적인 경로를 제공하지만, 표준 QPE는 긴 결맞음 Circuit이 필요하며 주로 오류 내성 양자 컴퓨터에 적합합니다. KQD는 이 두 접근 방식 사이에 있습니다. 위상 추정 기반 알고리즘처럼 실시간 해밀토니안 진화를 사용하지만, 완전한 위상 추정 대신 고전적으로 풀 수 있는 소형 투영 고유값 문제로 대체합니다.

$n$-qubit 해밀토니안 $H$와 기준 상태 $\lvert \psi_{0}\rangle$를 생각해봅시다. KQD 방법은 실시간 시간 진화된 상태로부터 Krylov 부분공간을 구성합니다.

$$\begin{equation*} \lvert \psi_\ell\rangle = e^{-i\ell\Delta t H}\lvert \psi_{0}\rangle, \qquad \ell = 0,1,\ldots,r-1, \end{equation*}$$

여기서 $r$은 Krylov 차원이고 $\Delta t$는 시간 스텝입니다. Krylov 부분공간의 모든 상태는 이러한 기저 상태의 선형 결합으로 표현됩니다.

$$\begin{equation*} \lvert \psi(\mathbf{c})\rangle = \frac{\sum_{\ell=0}^{r-1} c_\ell \lvert \psi_\ell\rangle} {\left\|\sum_{\ell=0}^{r-1} c_\ell \lvert \psi_\ell\rangle\right\|}, \end{equation*}$$

where the denominator normalizes the state.

간단한 대수를 통해 해당하는 에너지가 레일리 몫으로 표현됨을 알 수 있습니다.

$$\begin{equation*} E(\mathbf{c}) =\langle \psi(\mathbf{c})|H|\psi(\mathbf{c}) \rangle= \frac{ \sum_{k,\ell} c_k^* c_\ell \langle \psi_k\vert H\vert\psi_\ell\rangle }{ \sum_{k,\ell} c_k^* c_\ell \langle \psi_k\vert\psi_\ell\rangle } = \frac{ \mathbf{c}^{\dagger}\mathcal{H}\mathbf{c} }{ \mathbf{c}^{\dagger}\mathcal{S}\mathbf{c} }. \end{equation*}$$

여기서 행렬 $\mathcal{S}$와 $\mathcal{H}$는

$$\begin{equation*} \mathcal{S}_{k\ell}=\langle \psi_k\vert\psi_\ell\rangle, \qquad \mathcal{H}_{k\ell}=\langle \psi_k\vert H\vert\psi_\ell\rangle \end{equation*}$$

투영된 겹침 행렬과 해밀토니안 행렬을 정의합니다. 그 항목들은 양자 Circuit 측정을 사용하여 추정됩니다.

최소 $E(\mathbf{c})$를 주는 계수 $\mathbf{c}$를 구하고자 합니다.

$$\begin{equation*} \min_{\mathbf{c}\neq \mathbf{0}} E(\mathbf{c}). \end{equation*}$$

레일리-리츠 정리에 따르면, 이 최소화는 일반화 고유값 문제(GEVP)를 푸는 것과 동일합니다.

$$\begin{equation*} \mathcal{H}\mathbf{c}=E \mathcal{S}\mathbf{c}. \end{equation*}$$

차원 $r$은 고전 컴퓨터가 GEVP를 풀 수 있을 만큼 충분히 작을 수 있다는 점에 유의하세요.

이는 고전적 부분공간 대각화에서 사용되는 것과 동일한 변분 원리이지만, 여기서는 기저 상태가 양자 시간 진화에 의해 생성됩니다. VQE와 비교했을 때, KQD는 실시간 진화에 의존하기 때문에 일반적으로 더 깊은 Circuit이 필요합니다. 그 대신 KQD는 비선형 매개변수 최적화와 반복적인 양자 하드웨어 실행을 피하며, 투영된 부분공간이 확장됨에 따라 체계적으로 개선됩니다. 이 알고리즘은 기존 양자 하드웨어에서 대규모로 시연되었으며[2], 그 성능은 입증 가능한 보장과 함께 분석될 수 있습니다[1].

요구 사항

이 튜토리얼을 시작하기 전에 다음이 설치되어 있는지 확인하세요:

• Qiskit SDK v2.3 이상(시각화 지원 포함)

• Qiskit Runtime v0.22 이상 ( pip install qiskit-ibm-runtime )

• SciPy (pip install scipy)

• Matplotlib (pip install matplotlib)

• Pandas (pip install pandas)

하드웨어 실행에는 qiskit-ibm-runtime과 IBM Quantum® 계정에 대한 액세스가 필요합니다.

설정

설정 셀은 워크플로우에 필요한 모듈을 가져오고 헬퍼 함수를 정의합니다.

1. 하이젠베르크 해밀토니안 구축하기;

2. 임계값이 적용된 GEVP 풀기;

3. 학습된 Krylov 필터 평가하기;

4. 필터 값을 스펙트럼 가중치로 변환하기;

5. 기준 및 필터링된 에너지 분포 플로팅하기.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy pandas qiskit qiskit-ibm-runtime scipy

from __future__ import annotations

import warnings

import numpy as np
import pandas as pd
import scipy.linalg as la
import matplotlib.pyplot as plt

from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.primitives import StatevectorEstimator
from qiskit.quantum_info import Operator, SparsePauliOp
from qiskit.synthesis import LieTrotter, SuzukiTrotter
from qiskit.transpiler import PassManager, Layout
from qiskit.transpiler.passes import CommutativeOptimization
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2, Batch
from qiskit_ibm_runtime.fake_provider import FakeMarrakesh

warnings.filterwarnings("ignore")

def make_heisenberg_hamiltonian(
    num_qubits: int,
    coupling: float = 1.0,
) -> SparsePauliOp:
    """Make a Heisenberg Hamiltonian for a 1D chain of qubits with nearest-neighbor interactions."""
    terms: list[tuple[str, complex]] = []

    def append_term(q0: int, q1: int, pauli: str):
        label = ["I"] * num_qubits
        label[num_qubits - 1 - q0] = pauli[0]
        label[num_qubits - 1 - q1] = pauli[1]
        terms.append(("".join(label), coupling))

    for pauli in ("XX", "YY", "ZZ"):
        for q in range(num_qubits - 1):
            append_term(q, q + 1, pauli)
    return SparsePauliOp.from_list(terms).simplify()

def _basis_state_transition_amplitude_sparse(
    hamiltonian: SparsePauliOp,
    bra_state: int,
    ket_state: int,
) -> complex:
    """Evaluate <bra_state|H|ket_state> for computational-basis states."""
    num_qubits = hamiltonian.num_qubits
    amplitude = 0.0 + 0.0j
    for pauli, coeff in zip(hamiltonian.paulis, hamiltonian.coeffs):
        new_state = ket_state
        phase = 1.0 + 0.0j
        for q in range(num_qubits):
            x = bool(pauli.x[q])
            z = bool(pauli.z[q])
            if not x and not z:
                continue
            bit = (new_state >> q) & 1
            if x and z:
                # Y|0> = i|1>, Y|1> = -i|0>
                phase *= 1j if bit == 0 else -1j
                new_state ^= 1 << q
            elif x:
                new_state ^= 1 << q
            else:
                # Z|0> = |0>, Z|1> = -|1>
                if bit:
                    phase *= -1
        if new_state == bra_state:
            amplitude += coeff * phase
    return amplitude

def basis_state_expectation_sparse(
    hamiltonian: SparsePauliOp,
    bitstring: str,
) -> complex:
    """Evaluate <bitstring|H|bitstring>."""
    state = int(bitstring, 2)
    return _basis_state_transition_amplitude_sparse(
        hamiltonian,
        bra_state=state,
        ket_state=state,
    )

def diagonalize_single_1_subspace(
    hamiltonian: SparsePauliOp,
) -> np.ndarray:
    """Diagonalize the Hamiltonian projected onto the single-excitation subspace."""
    num_qubits = hamiltonian.num_qubits

    # Integer basis states |...010...>, with the excitation at qubit k.
    basis = [1 << k for k in range(num_qubits)]

    h_single = np.empty((num_qubits, num_qubits), dtype=complex)

    for row, bra_state in enumerate(basis):
        for col, ket_state in enumerate(basis):
            h_single[row, col] = _basis_state_transition_amplitude_sparse(
                hamiltonian,
                bra_state=bra_state,
                ket_state=ket_state,
            )
    # Remove floating-point-level asymmetry.
    h_single = 0.5 * (h_single + h_single.conj().T)
    evals, _ = np.linalg.eigh(h_single)
    return np.real(evals)

def simple_transpilation(circuit: QuantumCircuit) -> QuantumCircuit:
    """Transpilation to simplify the circuit"""
    pm = PassManager(
        [
            CommutativeOptimization(),
        ]
    )
    circuit = transpile(circuit, optimization_level=3)
    circuit = pm.run(circuit)
    return circuit

def summarize_circuit(circuit: QuantumCircuit) -> dict[str, int | str]:
    """Summarize the circuit with depth, size, and 2-qubit gate information."""
    two_qubit_total = sum(
        inst.operation.num_qubits == 2 for inst in circuit.data
    )
    two_qubit_depth = circuit.depth(lambda x: x[0].num_qubits == 2)
    return {
        "depth": circuit.depth(),
        "size": circuit.size(),
        "2q gates": two_qubit_total,
        "2q depth": two_qubit_depth,
    }

def solve_thresholded_gevp(
    h_matrix: np.ndarray,
    s_matrix: np.ndarray,
    threshold: float = 1e-10,
) -> tuple[float, np.ndarray, int]:
    """Solve H c = E S c using canonical orthogonalization of S."""
    s_vals, s_vecs = la.eigh(s_matrix)

    valid = s_vals > threshold
    if not np.any(valid):
        raise ValueError(
            "All overlap eigenvalues were removed by thresholding."
        )

    keep = valid

    orthogonalizer = s_vecs[:, keep] @ np.diag(1.0 / np.sqrt(s_vals[keep]))
    h_orth = orthogonalizer.conj().T @ h_matrix @ orthogonalizer
    h_orth = 0.5 * (h_orth + h_orth.conj().T)

    eigvals, eigvecs = la.eigh(h_orth)

    coeffs = orthogonalizer @ eigvecs[:, 0]
    normalization = np.sqrt(np.real(coeffs.conj().T @ s_matrix @ coeffs))
    coeffs /= normalization

    return float(np.real(eigvals[0])), coeffs, int(np.sum(keep))

이 튜토리얼의 첫 번째 부분에서는 로컬 상태 벡터 시뮬레이터를 사용하여 KQD 방법을 시연합니다. 이후에는 실제 양자 Backend를 사용하여 유틸리티 규모 문제를 다룹니다.

또한 Backend별 Transpile을 시연하고 결과 Circuit을 검사하기 위해 가짜 Backend를 정의합니다.

try:
    service = QiskitRuntimeService()
except Exception:
    QiskitRuntimeService.save_account(
        token="<api_token>", instance="<instance>", overwrite=True
    )
    service = QiskitRuntimeService()

backend = FakeMarrakesh()

소규모 시뮬레이터 예제

1단계: 고전적 입력을 양자 문제로 매핑하기

해밀토니안과 기준 상태

이 예제에서는 12-qubit의 열린 경계 하이젠베르크 사슬($n=12$)을 사용하며,

$$\begin{equation*} H=\sum_{i=0}^{n-2}\left(X_iX_{i+1}+Y_iY_{i+1}+Z_iZ_{i+1}\right), \end{equation*}$$

단일 여기 곱 상태를

$$\begin{equation*} |\psi_{0}\rangle=|000001000000\rangle \end{equation*}$$

기준 상태로 사용합니다. 위에서 정의한 하이젠베르크 해밀토니안은 총 여기 수를 보존하므로, 기준 상태는 qubit 수에 따라 선형적으로만 차원이 증가하는 단일 여기 부분공간에 머무릅니다. 따라서 해당 부분공간으로 제한된 해밀토니안을 대각화하여 정확한 바닥 상태 에너지를 효율적으로 계산할 수 있으며, 이를 순전히 KQD 추정치에 대한 진단 벤치마크로 사용합니다. KQD 워크플로우 자체는 Qiskit primitives를 사용하여 투영된 행렬 요소를 추정하고 결과로 나온 투영된 문제를 고전적으로 풉니다.

# Problem definition for the simulator example
num_qubits = 12
hamiltonian = make_heisenberg_hamiltonian(num_qubits=num_qubits, coupling=1.0)
ref_bitstring = "000001000000"
ref_energy = basis_state_expectation_sparse(hamiltonian, ref_bitstring)

print("Hamiltonian:")
print(hamiltonian)
print(f"Reference state: |{ref_bitstring}>")
print("Reference energy: ", ref_energy)

Hamiltonian:
SparsePauliOp(['IIIIIIIIIIXX', 'IIIIIIIIIXXI', 'IIIIIIIIXXII', 'IIIIIIIXXIII', 'IIIIIIXXIIII', 'IIIIIXXIIIII', 'IIIIXXIIIIII', 'IIIXXIIIIIII', 'IIXXIIIIIIII', 'IXXIIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIIYYI', 'IIIIIIIIYYII', 'IIIIIIIYYIII', 'IIIIIIYYIIII', 'IIIIIYYIIIII', 'IIIIYYIIIIII', 'IIIYYIIIIIII', 'IIYYIIIIIIII', 'IYYIIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIIZZI', 'IIIIIIIIZZII', 'IIIIIIIZZIII', 'IIIIIIZZIIII', 'IIIIIZZIIIII', 'IIIIZZIIIIII', 'IIIZZIIIIIII', 'IIZZIIIIIIII', 'IZZIIIIIIIII', 'ZZIIIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])
Reference state: |000001000000>
Reference energy:  (7+0j)

알고리즘의 매개변수 설정

해밀토니안 노름의 상한을 기반으로, 참고 문헌 [1]에서는 시간 스텝 $\Delta t$를 $\pi/\|H\|$로 어림잡아 제안합니다. 스펙트럼 노름 $\|H\|$는 계산하기 어려우므로, 대신 그 상한을 사용합니다.

$$\begin{equation*} \|H\| \le \sum_{i=0}^{n-2}\underbrace{\|X_i X_{i+1}+Y_{i} Y_{i+1}+Z_{i} Z_{i+1}\|}_{4 \times 4 \text{ matrix, easy to calculate the norm}} = 3(n-1). \end{equation*}$$

Krylov 차원을 $r=10$으로, 시간 스텝당 Trotter 스텝 수를 $5$로 설정합니다. 이는 가장 깊은 Circuit($t_{\max}=(r-1)\Delta t$)을 감당할 수 있는 수준으로 유지하면서 저준위 스펙트럼을 해석하기에 충분히 큰 Krylov 공간이며, 그 가장 깊은 Circuit에서 이산화 오류를 작게 유지하기에 충분한 Trotter 스텝 수입니다.

dt = np.pi / (3 * (num_qubits - 1))
print("dt in Krylov basis:  ", dt)

krylov_dim = 10
num_trotter_steps = 5

dt in Krylov basis:   0.09519977738150888

Circuit 구축하기

여기서는 행렬 요소 $\mathcal{H}_{k\ell}$와 $\mathcal{S}_{k\ell}$를 추정하기 위한 Circuit을 구축합니다. $H$의 모든 거듭제곱은 서로 교환되므로, 다음과 같습니다.

$$\begin{equation*} \mathcal{H}_{k\ell} = \langle\psi_{0}|H e^{-i(\ell-k)\Delta tH}|\psi_{0}\rangle=\mathcal{H}_{0,\ell-k}, \qquad \mathcal{S}_{k\ell} = \langle\psi_{0}|e^{-i(\ell-k)\Delta tH}|\psi_{0}\rangle=\mathcal{S}_{0,\ell-k}. \end{equation*}$$

이처럼 원소가 인덱스 차이에만 의존하는 행렬을 퇴플리츠(Toeplitz) 행렬이라고 하며, $d=\ell-k$로 인덱싱되는 첫 번째 행의 원소들로부터 전체를 재구성할 수 있습니다.

여기서는 확장 스왑 테스트라고 하는 Circuit을 제시하는데, 이는 다음을 준비합니다.

$$\begin{equation*} |\Phi_{0d}\rangle= \frac{|0\rangle|\psi_0\rangle+|1\rangle|\psi_d\rangle}{\sqrt{2}}, \end{equation*}$$

여기서 $|\psi_d\rangle=e^{-id\Delta tH}|\psi_{0}\rangle$입니다.

참조 상태

참조 상태 $|\psi_0\rangle$를 준비합니다.

qc_ref = QuantumCircuit(num_qubits)
for i, b in enumerate(reversed(ref_bitstring)):
    if b == "1":
        qc_ref.x(i)
display(qc_ref.draw("mpl", scale=0.5))

시간 진화

해밀토니안에 의해 생성되는 시간 진화 연산자를 간단한 Lie-Trotterization으로 근사하여 구현합니다.

t = Parameter("t")

evol_gate = PauliEvolutionGate(
    hamiltonian,
    time=t,
    synthesis=LieTrotter(reps=num_trotter_steps),
    label="U(t)",
)

# Synthesize U(t) first and then control the synthesized circuit.
# This makes the controlled structure visible in the circuit drawer.
evolution_circuit = QuantumCircuit(num_qubits, name="U(t)")
evolution_circuit.append(evol_gate, range(num_qubits))
evolution_circuit = simple_transpilation(evolution_circuit)

# Make a controlled version of the evolution circuit.
controlled_evolution_gate = evolution_circuit.to_gate(label="U(t)").control(
    1, label="C-U(t)"
)
display(
    evolution_circuit.assign_parameters({t: 1.5}).draw(
        "mpl", scale=0.5, fold=-1
    )
)

확장된 스왑 테스트 회로 [3]

회로는 먼저 시스템 레지스터에 참조 상태를 준비하는 동안 앤실라는 $|0\rangle$에 남아 있습니다:

$$\begin{equation*} |0\rangle |0\rangle^{\otimes n} \longrightarrow |0\rangle |\psi_{0}\rangle . \end{equation*}$$

그런 다음, 앤실라에 아다마르 게이트를 적용하면 두 가지 branch의 코히런트 중첩이 생성됩니다:

$$\begin{equation*} |0\rangle |\psi_0\rangle \longrightarrow \frac{|0\rangle + |1\rangle}{\sqrt{2}} |\psi_0\rangle = \frac{|0\rangle|\psi_0\rangle + |1\rangle|\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

마지막으로, 제어된 시간 진화 게이트가 적용됩니다:

$$\begin{equation*} U(t=d\Delta t) = e^{-id\Delta t H} \end{equation*}$$

앤실라가 $|1\rangle$ branch에 있을 때만 적용됩니다. 따라서,

$$\begin{equation*} \frac{|0\rangle|\psi_0\rangle + |1\rangle|\psi_0\rangle}{\sqrt{2}} \longrightarrow \frac{|0\rangle|\psi_0\rangle + |1\rangle U(d\Delta t)|\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

다음 코드 블록에서 우리는 다음을 구현합니다:

$$\begin{equation*} |\Psi(t)\rangle = \frac{|0\rangle|\psi_0\rangle + |1\rangle U(t)|\psi_0\rangle}{\sqrt{2}}, \end{equation*}$$

이는 실행 단계에서 $d=1,\cdots r-1$에 대해 $t=d\Delta t$로 할당됩니다.

ancilla = 0
system_qubits = list(range(1, num_qubits + 1))

extended_swap_test = QuantumCircuit(num_qubits + 1)

# Append state preparation part
extended_swap_test = extended_swap_test.compose(qc_ref, system_qubits)

# Prepare the coherent branch label, (|0> + |1>) / sqrt(2).
extended_swap_test.h(ancilla)

# Apply U(t) only to the |1> branch of the ancilla.
extended_swap_test.append(
    controlled_evolution_gate, [ancilla] + system_qubits
)

# Decompose once more for visualization so that control bullets are visible.
display(extended_swap_test.draw("mpl", fold=-1))

관측 가능량

임의의 에르미트 시스템 관측 가능량 $O$에 대해, 여기서는 계산할 관측 가능량을 다음과 같이 설정합니다:

$$\begin{equation*} z=\langle \psi_0|O|\psi_d\rangle . \end{equation*}$$

이는 $O=I$가 겹침 요소 $\mathcal{S}_{0d}=\langle\psi_0|\psi_d\rangle$를 제공하는 반면, $O=H$는 해밀토니안 요소 $\mathcal{H}_{0d}=\langle\psi_0|H|\psi_d\rangle$를 제공하기 때문입니다.

$X=|0\rangle\langle 1|+|1\rangle\langle 0|$를 사용하면, 다음을 얻습니다

$$\begin{align*} \langle \Psi_{0d}|X\otimes O| \Psi_{0d} \rangle &= \left(\frac{\langle0|\langle\psi_0| + \langle1| \langle\psi_d|}{\sqrt{2}}\right)(X\otimes O)\left(\frac{|0\rangle|\psi_0\rangle + |1\rangle |\psi_d\rangle}{\sqrt{2}}\right) \\ &= \frac{1}{2} \left( \langle \psi_0|O|\psi_d\rangle + \langle \psi_d|O|\psi_0\rangle \right) \\ &= \frac{z+z^*}{2} = \operatorname{Re} z. \end{align*}$$

마찬가지로, $Y=-i|0\rangle\langle 1|+i|1\rangle\langle 0|$를 사용하면,

$$\begin{align*} \langle \Psi_{0d}|Y\otimes O| \Psi_{0d} \rangle &= \left(\frac{\langle0|\langle\psi_0| + \langle1| \langle\psi_d|}{\sqrt{2}}\right)(Y\otimes O)\left(\frac{|0\rangle|\psi_0\rangle + |1\rangle |\psi_d\rangle}{\sqrt{2}}\right) \\ &= \frac{1}{2} \left( -i\langle \psi_0|O|\psi_d\rangle + i\langle \psi_d|O|\psi_0\rangle \right) \\ &= \frac{-iz+iz^*}{2} = \operatorname{Im} z. \end{align*}$$

따라서, 다음을 얻습니다

$$\begin{equation*} \langle \Phi_{0d}|X\otimes O|\Phi_{0d}\rangle = \operatorname{Re}\langle\psi_0|O|\psi_d\rangle , \quad \langle \Phi_{0d}|Y\otimes O|\Phi_{0d}\rangle = \operatorname{Im}\langle\psi_0|O|\psi_d\rangle . \end{equation*}$$

마지막으로, 각 상태 $|\Psi_{0d}\rangle$에 대해, 다음을 측정해야 합니다:

$$\begin{align*} X \otimes I \text{ for } \operatorname{Re}\mathcal{S}_{0d},\\ Y \otimes I \text{ for } \operatorname{Im}\mathcal{S}_{0d},\\ X \otimes H \text{ for } \operatorname{Re}\mathcal{H}_{0d},\\ Y \otimes H \text{ for } \operatorname{Im}\mathcal{H}_{0d}.\\ \end{align*}$$

n_qubits = hamiltonian.num_qubits

observable_labels = [
    "Re S_0d",
    "Im S_0d",
    "Re H_0d",
    "Im H_0d",
]

# X ⊗ I and Y ⊗ I.
# Qiskit's Pauli-label convention places qubit 0 on the rightmost character,
# so the ancilla Pauli is appended to the right.
obs_x_identity = SparsePauliOp("I" * n_qubits + "X")
obs_y_identity = SparsePauliOp("I" * n_qubits + "Y")

# X ⊗ H and Y ⊗ H.
obs_x_hamiltonian = SparsePauliOp.from_list(
    [
        (label + "X", coeff)
        for label, coeff in zip(
            hamiltonian.paulis.to_labels(),
            hamiltonian.coeffs,
        )
    ]
)

obs_y_hamiltonian = SparsePauliOp.from_list(
    [
        (label + "Y", coeff)
        for label, coeff in zip(
            hamiltonian.paulis.to_labels(),
            hamiltonian.coeffs,
        )
    ]
)

observables = [
    obs_x_identity,
    obs_y_identity,
    obs_x_hamiltonian,
    obs_y_hamiltonian,
]

for obs, label in zip(observables, observable_labels):
    print(f"Observable: {label}")
    print(obs)
    print()

Observable: Re S_0d
SparsePauliOp(['IIIIIIIIIIIIX'],
              coeffs=[1.+0.j])

Observable: Im S_0d
SparsePauliOp(['IIIIIIIIIIIIY'],
              coeffs=[1.+0.j])

Observable: Re H_0d
SparsePauliOp(['IIIIIIIIIIXXX', 'IIIIIIIIIXXIX', 'IIIIIIIIXXIIX', 'IIIIIIIXXIIIX', 'IIIIIIXXIIIIX', 'IIIIIXXIIIIIX', 'IIIIXXIIIIIIX', 'IIIXXIIIIIIIX', 'IIXXIIIIIIIIX', 'IXXIIIIIIIIIX', 'XXIIIIIIIIIIX', 'IIIIIIIIIIYYX', 'IIIIIIIIIYYIX', 'IIIIIIIIYYIIX', 'IIIIIIIYYIIIX', 'IIIIIIYYIIIIX', 'IIIIIYYIIIIIX', 'IIIIYYIIIIIIX', 'IIIYYIIIIIIIX', 'IIYYIIIIIIIIX', 'IYYIIIIIIIIIX', 'YYIIIIIIIIIIX', 'IIIIIIIIIIZZX', 'IIIIIIIIIZZIX', 'IIIIIIIIZZIIX', 'IIIIIIIZZIIIX', 'IIIIIIZZIIIIX', 'IIIIIZZIIIIIX', 'IIIIZZIIIIIIX', 'IIIZZIIIIIIIX', 'IIZZIIIIIIIIX', 'IZZIIIIIIIIIX', 'ZZIIIIIIIIIIX'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

Observable: Im H_0d
SparsePauliOp(['IIIIIIIIIIXXY', 'IIIIIIIIIXXIY', 'IIIIIIIIXXIIY', 'IIIIIIIXXIIIY', 'IIIIIIXXIIIIY', 'IIIIIXXIIIIIY', 'IIIIXXIIIIIIY', 'IIIXXIIIIIIIY', 'IIXXIIIIIIIIY', 'IXXIIIIIIIIIY', 'XXIIIIIIIIIIY', 'IIIIIIIIIIYYY', 'IIIIIIIIIYYIY', 'IIIIIIIIYYIIY', 'IIIIIIIYYIIIY', 'IIIIIIYYIIIIY', 'IIIIIYYIIIIIY', 'IIIIYYIIIIIIY', 'IIIYYIIIIIIIY', 'IIYYIIIIIIIIY', 'IYYIIIIIIIIIY', 'YYIIIIIIIIIIY', 'IIIIIIIIIIZZY', 'IIIIIIIIIZZIY', 'IIIIIIIIZZIIY', 'IIIIIIIZZIIIY', 'IIIIIIZZIIIIY', 'IIIIIZZIIIIIY', 'IIIIZZIIIIIIY', 'IIIZZIIIIIIIY', 'IIZZIIIIIIIIY', 'IZZIIIIIIIIIY', 'ZZIIIIIIIIIIY'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

$\mathcal{H}$ 행렬 요소의 추정을 위해서는, Pauli 항의 개수가 $\mathcal{S}$ 행렬 요소보다 훨씬 많습니다.

이제 시프팅 기법 [4]을 사용하여 측정해야 하는 해밀토니안 항의 개수를 줄일 수 있습니다. 해밀토니안을 다음과 같이 분리합니다:

$$\begin{equation*} H = (H-T) + T, \end{equation*}$$

여기서 $T$는 참조 상태가 그 고유 상태가 되도록 선택됩니다,

$$\begin{equation*} T|\psi_0\rangle = \tau |\psi_0\rangle . \end{equation*}$$

Then,

$$\begin{align*} \mathcal{H}_{0d} &= \langle\psi_0|H|\psi_d\rangle \\ &= \langle\psi_0|(H-T)|\psi_d\rangle + \langle\psi_0|T|\psi_d\rangle \\ &= \langle\psi_0|(H-T)|\psi_d\rangle + \tau \langle\psi_0|\psi_d\rangle \\ &= \widetilde{\mathcal{H}}_{0d} + \tau \mathcal{S}_{0d}. \end{align*}$$

Here,

$$\begin{equation*} \widetilde{\mathcal{H}}_{0d} = \langle\psi_0|(H-T)|\psi_d\rangle \end{equation*}$$

는 시프트된 해밀토니안 행렬 요소입니다. 따라서, $X\otimes (H-T)$와 $Y\otimes (H-T)$만 측정하면 됩니다. $T$로부터의 기여는 이미 측정된 겹침 행렬 요소 $\mathcal{S}_{0d}$를 사용하여 고전적으로 재구성됩니다.

이 예제에서, 자연스러운 선택은 하이젠베르크 해밀토니안의 대각 부분입니다,

$$\begin{equation*} T = \sum_{i=0}^{n-2} Z_iZ_{i+1}. \end{equation*}$$

참조 상태는 계산 기저 상태이므로, 모든 $Z_iZ_{i+1}$ 항의 고유 상태입니다.

그러나, 더 유리한 선택은 대각 $ZZ$ 항뿐만 아니라 참조 상태를 소멸시키는 $XX+YY$ 항도 포함하는 것입니다.

각 인접 쌍에 대해, 연산자 $XX+YY$는 다음을 만족합니다

$$\begin{equation*} (XX+YY)|00\rangle = 0, \qquad (XX+YY)|11\rangle = 0, \end{equation*}$$

and

$$\begin{equation*} (XX+YY)|01\rangle = 2|10\rangle, \qquad (XX+YY)|10\rangle = 2|01\rangle . \end{equation*}$$

따라서, $XX+YY$ 항은 두 인접 큐비트가 참조 비트열에서 서로 다른 점유 상태를 가질 때만 기여합니다. 두 큐비트가 모두 $0$이거나 모두 $1$이면, 이 항은 참조 상태를 소멸시키므로 마찬가지로 시프트되어 제거될 수 있습니다.

$|\psi_{\rm ref}\rangle = |z_0\cdots z_{n-1}\rangle$라고 하고, 여기서 $z_i \in \{0,1\}$입니다. 따라서, 다음을 선택할 수 있습니다

$$\begin{equation*} T=\sum_{i=0}^{n-2} Z_iZ_{i+1} + \sum_{z_i={z_{i+1}}} \left( X_iX_{i+1} + Y_iY_{i+1} \right). \end{equation*}$$

이 연산자는 여전히 다음을 만족합니다

$$\begin{equation*} T|\psi_0\rangle = \tau |\psi_0\rangle , \end{equation*}$$

이는 $ZZ$ 항이 $|\psi_0\rangle$에 대각적으로 작용하는 반면, 시프트된 $XX+YY$ 항은 0을 주기 때문입니다. 따라서 해당 고유값은 $ZZ$ 항에 의해서만 결정됩니다,

$$\begin{equation*} \tau = \sum_{i=0}^{n-2} (-1)^{z_i} (-1)^{z_{i+1}} . \end{equation*}$$

이 선택으로, 시프트된 해밀토니안은 다음과 같이 됩니다:

$$\begin{equation*} H-T = \sum_{z_i\ne z_{i+1}} \left( X_iX_{i+1} + Y_iY_{i+1} \right). \end{equation*}$$

결과적으로, 참조 상태에서 서로 다른 점유 상태를 갖는 엣지만 측정하면 됩니다. 모든 $ZZ$ 항과 모든 비활성 $XX+YY$ 항은 겹침 기여 $\tau \mathcal{S}_{0d}$를 통해 재구성되거나, 구성상 0의 기여를 줍니다.

이는 대각 부분만 시프트하는 것보다 더 작은 관측 가능량을 제공합니다. 특히, 국소화된 여기를 갖는 계산 기저 참조 상태의 경우, 여기에 인접한 엣지만 $H-T$에 남습니다. 따라서, $X\otimes(H-T)$와 $Y\otimes(H-T)$의 Pauli 항 개수는 상당히 줄어들 수 있는 반면, 재구성된 행렬 요소

$$\begin{equation*} \widetilde{\mathcal{H}}_{0d} + \tau \mathcal{S}_{0d} \end{equation*}$$

는 정확히 동일하게 유지됩니다.

def make_reduced_heisenberg_observables(
    ref_bitstring: str,
    coupling: float = 1.0,
) -> tuple[SparsePauliOp, SparsePauliOp, float]:
    """Build X⊗(H-T), Y⊗(H-T), and tau."""
    n_qubits = len(ref_bitstring)

    shifted_terms: list[tuple[str, complex]] = []
    tau = 0.0

    def bit(q: int) -> str:
        return ref_bitstring[n_qubits - 1 - q]

    def append_term(q0: int, q1: int, pauli: str):
        label = ["I"] * n_qubits
        label[n_qubits - 1 - q0] = pauli[0]
        label[n_qubits - 1 - q1] = pauli[1]
        shifted_terms.append(("".join(label), coupling))

    for q in range(n_qubits - 1):
        same_occupation = bit(q) == bit(q + 1)

        # ZZ contribution to tau
        tau += coupling * (1.0 if same_occupation else -1.0)

        # XX + YY survives only for opposite occupations.
        if not same_occupation:
            append_term(q, q + 1, "XX")
            append_term(q, q + 1, "YY")

    if shifted_terms:
        obs_x_shifted_hamiltonian = SparsePauliOp.from_list(
            [(label + "X", coeff) for label, coeff in shifted_terms]
        )
        obs_y_shifted_hamiltonian = SparsePauliOp.from_list(
            [(label + "Y", coeff) for label, coeff in shifted_terms]
        )
    else:
        obs_x_shifted_hamiltonian = SparsePauliOp(
            "I" * n_qubits + "X", coeffs=[0.0]
        )
        obs_y_shifted_hamiltonian = SparsePauliOp(
            "I" * n_qubits + "Y", coeffs=[0.0]
        )

    return obs_x_shifted_hamiltonian, obs_y_shifted_hamiltonian, tau

obs_x_shifted_hamiltonian, obs_y_shifted_hamiltonian, shift_tau = (
    make_reduced_heisenberg_observables(ref_bitstring)
)

print("Observable: Re shifted H_0d")
print(obs_x_shifted_hamiltonian)
print()

print("Observable: Im shifted H_0d")
print(obs_y_shifted_hamiltonian)
print()

print("tau =", shift_tau)

observables = [
    obs_x_identity,
    obs_y_identity,
    obs_x_shifted_hamiltonian,
    obs_y_shifted_hamiltonian,
]
observable_labels = [
    "Re S_0d",
    "Im S_0d",
    "Re shifted H_0d",
    "Im shifted H_0d",
]

Observable: Re shifted H_0d
SparsePauliOp(['IIIIIXXIIIIIX', 'IIIIIYYIIIIIX', 'IIIIXXIIIIIIX', 'IIIIYYIIIIIIX'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

Observable: Im shifted H_0d
SparsePauliOp(['IIIIIXXIIIIIY', 'IIIIIYYIIIIIY', 'IIIIXXIIIIIIY', 'IIIIYYIIIIIIY'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

tau = 7.0

Step 2: 양자 하드웨어 실행을 위한 문제 최적화

이제 추상적인 확장 스왑 테스트 회로를 하드웨어 지향 템플릿으로 전환합니다. 그 전에, 먼저 추상적인 수준에서 회로를 더 최적화합니다.

해밀토니안 항 순서 비교

먼저, 해밀토니안 시뮬레이션을 위해 하이젠베르크 해밀토니안에서 Pauli 항의 여러 순서를 비교합니다. 해밀토니안 자체는 변경되지 않지만, 순서는 곱 공식 회로가 생성되는 방식과 회로에서 구조가 병렬화될 수 있는 정도에 영향을 미칩니다. 예를 들어, 단순한 순서는 모든 최근접 이웃 $XX$ 항을 나열한 다음, 모든 $YY$ 항, 그다음 모든 $ZZ$ 항을 나열합니다. 이는 $(0,1)$과 $(1,2)$와 같은 인접한 엣지를 서로 옆에 배치하므로, 병렬로 실행할 수 없습니다. 짝수-후-홀수 순서는 먼저 분리된 짝수 엣지를 방문한 다음 홀수 엣지를 방문하여, 병렬 2-큐비트 레이어를 노출합니다. 짝수-홀수 엣지 그룹화 순서는 한 단계 더 나아갑니다: 각 엣지에 대해, 로컬 $XX$, $YY$, $ZZ$ 항을 함께 유지하면서도 짝수 엣지를 홀수 엣지보다 먼저 방문합니다. 짝수-후-홀수 및 짝수-홀수 엣지 그룹화 순서는 병렬 2-큐비트 레이어를 노출시켜 회로 깊이를 줄일 것으로 예상되며, 엣지 그룹화 순서는 동일한 엣지에서의 로컬 2-큐비트 상호작용이 컴팩트한 블록으로 처리되므로 Trotter 오차를 추가로 줄일 것으로 예상됩니다.

해밀토니안 순서는 Trotter 오차에도 영향을 미칩니다. 교환하지 않는 항을 서로 옆에 배치하면, 기저 전환이 더 자주 발생하여 더 많은 Trotter 오차를 유발합니다. 동일한 Pauli 기저 변환이 필요한 항을 그룹화함으로써, 불필요한 기저 변경을 피할 수 있습니다.

여기서, 비교는 동일한 트랜스파일 조건을 사용하여 첫 번째 행 Krylov 추정값에 나타나는 가장 큰 진화 시간, $t_{\max} = (r-1)\Delta t$를 사용합니다.

Trotter 오차를 측정하기 위해, Trotter화된 회로와 정확한 해밀토니안 진화 사이의 **프로세스 불충실도(process infidelity)**를 사용합니다,

$$\begin{equation*} \text{Infidelity}(U,V) = 1-F(U,V) = 1 - \frac{|\operatorname{Tr}(U^\dagger V)|^2}{d^2}, \end{equation*}$$

where $d=2^n$ is the dimension of the Hilbert space.

이 진단은 밀집 행렬을 사용하므로, 이 작은 12-큐비트 예제에는 적합하지만 확장 가능한 서브루틴으로 의도된 것은 아닙니다.

# The comparison uses the largest time that appears in the first-row Krylov estimates.
# Circuit depth does not depend on this numeric value, but the Trotter error does.
comparison_time = (krylov_dim - 1) * dt

def make_heisenberg_hamiltonian_ordered(
    num_qubits: int,
    ordering: str,
    coupling: float = 1.0,
) -> SparsePauliOp:
    """Return the same Heisenberg Hamiltonian with a specified term ordering."""
    terms: list[tuple[str, complex]] = []
    even_edges = [(q, q + 1) for q in range(0, num_qubits - 1, 2)]
    odd_edges = [(q, q + 1) for q in range(1, num_qubits - 1, 2)]

    def append_term(q0: int, q1: int, pauli: str):
        label = ["I"] * num_qubits
        label[num_qubits - 1 - q0] = pauli[0]
        label[num_qubits - 1 - q1] = pauli[1]
        terms.append(("".join(label), coupling))

    if ordering == "naive":
        for pauli in ("XX", "YY", "ZZ"):
            for q in range(num_qubits - 1):
                append_term(q, q + 1, pauli)
    elif ordering == "even-then-odd":
        for pauli in ("XX", "YY", "ZZ"):
            for q0, q1 in even_edges + odd_edges:
                append_term(q0, q1, pauli)
    elif ordering == "even-odd edge-grouped":
        for q0, q1 in even_edges + odd_edges:
            for pauli in ("XX", "YY", "ZZ"):
                append_term(q0, q1, pauli)
    else:
        raise ValueError(f"Unknown ordering: {ordering}")

    return SparsePauliOp.from_list(terms).simplify()

def build_numeric_evolution_circuit(
    hamiltonian: SparsePauliOp,
    synthesis,
    time_value: float,
    **synthesis_kwargs,
) -> QuantumCircuit:
    """Build a numeric circuit for exp(-i H t) with a chosen synthesis rule."""
    evolution_gate = PauliEvolutionGate(
        hamiltonian,
        time=time_value,
        synthesis=synthesis(**synthesis_kwargs),
    )
    circuit = QuantumCircuit(hamiltonian.num_qubits)
    circuit.append(evolution_gate, range(hamiltonian.num_qubits))
    return circuit

def process_infidelity(
    circuit: QuantumCircuit,
    exact_matrix: np.ndarray,
) -> float:
    """Return 1 - |Tr(U_circuit† U_exact) / d|²."""
    circuit_matrix = np.asarray(Operator(circuit).data)
    dim = circuit_matrix.shape[0]
    normalized_trace = np.vdot(circuit_matrix, exact_matrix) / dim
    fidelity = np.abs(normalized_trace) ** 2
    return float(np.clip(1.0 - fidelity, 0.0, 1.0))

hamiltonians_by_ordering = {
    ordering: make_heisenberg_hamiltonian_ordered(
        num_qubits, ordering, coupling=1.0
    )
    for ordering in ["naive", "even-then-odd", "even-odd edge-grouped"]
}

print(f"Comparison time: {comparison_time}\n")

for order_name, ham_ordered in hamiltonians_by_ordering.items():
    print(f"{order_name}:")
    print([op for op, _ in ham_ordered.to_list()])
    print()

print("Precomputing the exact evolution operator... ", end="")
exact_matrix = la.expm(-1j * comparison_time * hamiltonian.to_matrix())
print("Done")

Comparison time: 0.8567979964335799

naive:
['IIIIIIIIIIXX', 'IIIIIIIIIXXI', 'IIIIIIIIXXII', 'IIIIIIIXXIII', 'IIIIIIXXIIII', 'IIIIIXXIIIII', 'IIIIXXIIIIII', 'IIIXXIIIIIII', 'IIXXIIIIIIII', 'IXXIIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIIYYI', 'IIIIIIIIYYII', 'IIIIIIIYYIII', 'IIIIIIYYIIII', 'IIIIIYYIIIII', 'IIIIYYIIIIII', 'IIIYYIIIIIII', 'IIYYIIIIIIII', 'IYYIIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIIZZI', 'IIIIIIIIZZII', 'IIIIIIIZZIII', 'IIIIIIZZIIII', 'IIIIIZZIIIII', 'IIIIZZIIIIII', 'IIIZZIIIIIII', 'IIZZIIIIIIII', 'IZZIIIIIIIII', 'ZZIIIIIIIIII']

even-then-odd:
['IIIIIIIIIIXX', 'IIIIIIIIXXII', 'IIIIIIXXIIII', 'IIIIXXIIIIII', 'IIXXIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIXXIII', 'IIIIIXXIIIII', 'IIIXXIIIIIII', 'IXXIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIYYII', 'IIIIIIYYIIII', 'IIIIYYIIIIII', 'IIYYIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIYYI', 'IIIIIIIYYIII', 'IIIIIYYIIIII', 'IIIYYIIIIIII', 'IYYIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIZZII', 'IIIIIIZZIIII', 'IIIIZZIIIIII', 'IIZZIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIZZI', 'IIIIIIIZZIII', 'IIIIIZZIIIII', 'IIIZZIIIIIII', 'IZZIIIIIIIII']

even-odd edge-grouped:
['IIIIIIIIIIXX', 'IIIIIIIIIIYY', 'IIIIIIIIIIZZ', 'IIIIIIIIXXII', 'IIIIIIIIYYII', 'IIIIIIIIZZII', 'IIIIIIXXIIII', 'IIIIIIYYIIII', 'IIIIIIZZIIII', 'IIIIXXIIIIII', 'IIIIYYIIIIII', 'IIIIZZIIIIII', 'IIXXIIIIIIII', 'IIYYIIIIIIII', 'IIZZIIIIIIII', 'XXIIIIIIIIII', 'YYIIIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIIIYYI', 'IIIIIIIIIZZI', 'IIIIIIIXXIII', 'IIIIIIIYYIII', 'IIIIIIIZZIII', 'IIIIIXXIIIII', 'IIIIIYYIIIII', 'IIIIIZZIIIII', 'IIIXXIIIIIII', 'IIIYYIIIIIII', 'IIIZZIIIIIII', 'IXXIIIIIIIII', 'IYYIIIIIIIII', 'IZZIIIIIIIII']

Precomputing the exact evolution operator... Done

먼저 합성 규칙을 하나의 1차 Trotter 스텝으로 고정하고 Pauli 항 순서만 변경합니다. 이 비교의 목표는 주로 분리된 최근접 이웃 엣지를 트랜스파일러에 노출함으로써 회로 깊이와 2-큐비트 비용을 얼마나 줄일 수 있는지 확인하는 것입니다.

ordering_comparison_rows = []
infidelity_reps = [1, 2, 4, 8]

for ordering, ham_ordered in hamiltonians_by_ordering.items():
    for reps in infidelity_reps:
        circuit = build_numeric_evolution_circuit(
            ham_ordered,
            LieTrotter,
            comparison_time,
            reps=reps,
        )
        decomposed_circuit = simple_transpilation(circuit)
        infidelity = process_infidelity(decomposed_circuit, exact_matrix)
        ordering_comparison_rows.append(
            {
                "ordering": ordering,
                "synthesis": f"LieTrotter(reps={reps})",
                "infidelity": infidelity,
                **summarize_circuit(decomposed_circuit),
            }
        )
        if reps == 1:
            print(f"Circuit for {ordering} ordering:")
            print([op for op, _ in ham_ordered.to_list()])
            display(decomposed_circuit.draw("mpl", fold=-1, scale=0.6))

ordering_comparison_df = pd.DataFrame(ordering_comparison_rows)
display(ordering_comparison_df)

hamiltonian_for_synthesis = hamiltonians_by_ordering["even-odd edge-grouped"]

Circuit for naive ordering:
['IIIIIIIIIIXX', 'IIIIIIIIIXXI', 'IIIIIIIIXXII', 'IIIIIIIXXIII', 'IIIIIIXXIIII', 'IIIIIXXIIIII', 'IIIIXXIIIIII', 'IIIXXIIIIIII', 'IIXXIIIIIIII', 'IXXIIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIIYYI', 'IIIIIIIIYYII', 'IIIIIIIYYIII', 'IIIIIIYYIIII', 'IIIIIYYIIIII', 'IIIIYYIIIIII', 'IIIYYIIIIIII', 'IIYYIIIIIIII', 'IYYIIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIIZZI', 'IIIIIIIIZZII', 'IIIIIIIZZIII', 'IIIIIIZZIIII', 'IIIIIZZIIIII', 'IIIIZZIIIIII', 'IIIZZIIIIIII', 'IIZZIIIIIIII', 'IZZIIIIIIIII', 'ZZIIIIIIIIII']

Circuit for even-then-odd ordering:
['IIIIIIIIIIXX', 'IIIIIIIIXXII', 'IIIIIIXXIIII', 'IIIIXXIIIIII', 'IIXXIIIIIIII', 'XXIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIXXIII', 'IIIIIXXIIIII', 'IIIXXIIIIIII', 'IXXIIIIIIIII', 'IIIIIIIIIIYY', 'IIIIIIIIYYII', 'IIIIIIYYIIII', 'IIIIYYIIIIII', 'IIYYIIIIIIII', 'YYIIIIIIIIII', 'IIIIIIIIIYYI', 'IIIIIIIYYIII', 'IIIIIYYIIIII', 'IIIYYIIIIIII', 'IYYIIIIIIIII', 'IIIIIIIIIIZZ', 'IIIIIIIIZZII', 'IIIIIIZZIIII', 'IIIIZZIIIIII', 'IIZZIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIZZI', 'IIIIIIIZZIII', 'IIIIIZZIIIII', 'IIIZZIIIIIII', 'IZZIIIIIIIII']

Circuit for even-odd edge-grouped ordering:
['IIIIIIIIIIXX', 'IIIIIIIIIIYY', 'IIIIIIIIIIZZ', 'IIIIIIIIXXII', 'IIIIIIIIYYII', 'IIIIIIIIZZII', 'IIIIIIXXIIII', 'IIIIIIYYIIII', 'IIIIIIZZIIII', 'IIIIXXIIIIII', 'IIIIYYIIIIII', 'IIIIZZIIIIII', 'IIXXIIIIIIII', 'IIYYIIIIIIII', 'IIZZIIIIIIII', 'XXIIIIIIIIII', 'YYIIIIIIIIII', 'ZZIIIIIIIIII', 'IIIIIIIIIXXI', 'IIIIIIIIIYYI', 'IIIIIIIIIZZI', 'IIIIIIIXXIII', 'IIIIIIIYYIII', 'IIIIIIIZZIII', 'IIIIIXXIIIII', 'IIIIIYYIIIII', 'IIIIIZZIIIII', 'IIIXXIIIIIII', 'IIIYYIIIIIII', 'IIIZZIIIIIII', 'IXXIIIIIIIII', 'IYYIIIIIIIII', 'IZZIIIIIIIII']

ordering           synthesis  infidelity  depth  size  \
0                   naive  LieTrotter(reps=1)    0.999917     15    33   
1                   naive  LieTrotter(reps=2)    0.805998     21    66   
2                   naive  LieTrotter(reps=4)    0.271389     33   132   
3                   naive  LieTrotter(reps=8)    0.074590     57   264   
4           even-then-odd  LieTrotter(reps=1)    0.999917      6    33   
5           even-then-odd  LieTrotter(reps=2)    0.805998     12    66   
6           even-then-odd  LieTrotter(reps=4)    0.271389     24   132   
7           even-then-odd  LieTrotter(reps=8)    0.074590     48   264   
8   even-odd edge-grouped  LieTrotter(reps=1)    0.998432      6    33   
9   even-odd edge-grouped  LieTrotter(reps=2)    0.653843     12    66   
10  even-odd edge-grouped  LieTrotter(reps=4)    0.181529     24   132   
11  even-odd edge-grouped  LieTrotter(reps=8)    0.045244     48   264   

    2q gates  2q depth  
0         33        15  
1         66        21  
2        132        33  
3        264        57  
4         33         6  
5         66        12  
6        132        24  
7        264        48  
8         33         6  
9         66        12  
10       132        24  
11       264        48

reps=1에서, 짝수-후-홀수 및 짝수-홀수 엣지 그룹화 순서 모두 병렬 2-큐비트 레이어를 노출시켜 깊이를 15에서 6으로 축소하는 반면, 2-큐비트 게이트 개수는 세 가지 순서 모두에서 동일하게 유지됨을 관찰합니다. 그러나, 세 가지 순서 모두 불충실도가 1에 가까우므로, 순서를 더 명확하게 구분하기 위해 Trotter 반복 횟수를 스윕합니다. 반복 횟수가 증가함에 따라, even-odd edge-grouped 순서의 불충실도가 다른 두 순서보다 더 빠르게 감소하여, reps=8에서 0.045에 도달하는 반면 단순 순서와 짝수-후-홀수 순서는 0.075입니다.

곱 공식 합성 비교

다음으로, 해밀토니안 순서를 짝수-홀수 엣지 그룹화 순서로 고정한 상태에서 Trotter화의 다양한 고급 설정을 살펴봅니다. 1차 Lie-Trotter, 2차 Suzuki-Trotter, 4차 Suzuki-Trotter를 고려합니다.

synthesis_comparison_rows = []

for num_trotter_steps in [1, 2, 3, 4, 5]:
    synthesis_cases = [
        ("LieTrotter", LieTrotter, {"reps": num_trotter_steps}),
        (
            "SuzukiTrotter(order=2)",
            SuzukiTrotter,
            {"order": 2, "reps": num_trotter_steps},
        ),
        (
            "SuzukiTrotter(order=4)",
            SuzukiTrotter,
            {"order": 4, "reps": num_trotter_steps},
        ),
    ]
    for label, synthesis, kwargs in synthesis_cases:
        circuit = build_numeric_evolution_circuit(
            hamiltonian_for_synthesis,
            synthesis,
            comparison_time,
            **kwargs,
        )
        decomposed_circuit = simple_transpilation(circuit)
        synthesis_comparison_rows.append(
            {
                "synthesis": label,
                "reps": kwargs["reps"],
                "infidelity": process_infidelity(
                    decomposed_circuit, exact_matrix
                ),
                **summarize_circuit(decomposed_circuit),
            }
        )

synthesis_comparison_df = pd.DataFrame(synthesis_comparison_rows)
display(synthesis_comparison_df.sort_values(["2q gates", "2q depth"]))

# For memory free
exact_matrix = None

synthesis  reps    infidelity  depth  size  2q gates  \
0               LieTrotter     1  9.984324e-01      6    33        33   
1   SuzukiTrotter(order=2)     1  9.733399e-01      9    51        51   
3               LieTrotter     2  6.538427e-01     12    66        66   
4   SuzukiTrotter(order=2)     2  2.522533e-01     15    84        84   
6               LieTrotter     3  3.197242e-01     18    99        99   
7   SuzukiTrotter(order=2)     3  4.804050e-02     21   117       117   
9               LieTrotter     4  1.815291e-01     24   132       132   
10  SuzukiTrotter(order=2)     4  1.453103e-02     27   150       150   
12              LieTrotter     5  1.161770e-01     30   165       165   
2   SuzukiTrotter(order=4)     1  2.884402e-01     33   183       183   
13  SuzukiTrotter(order=2)     5  5.803455e-03     33   183       183   
5   SuzukiTrotter(order=4)     2  1.641162e-03     63   348       348   
8   SuzukiTrotter(order=4)     3  2.907076e-05     93   513       513   
11  SuzukiTrotter(order=4)     4  2.791061e-06    123   678       678   
14  SuzukiTrotter(order=4)     5  4.736685e-07    153   843       843   

    2q depth  
0          6  
1          9  
3         12  
4         15  
6         18  
7         21  
9         24  
10        27  
12        30  
2         33  
13        33  
5         63  
8         93  
11       123  
14       153

Lie-Trotter는 가장 얕은 회로를 제공하지만 가장 큰 오차를 가지는 반면, 4차 Suzuki-Trotter는 더 정확하지만 회로 깊이를 증가시킵니다. 이 튜토리얼의 나머지 부분에서는, 1차 공식에 비해 Trotter 오차를 상당히 줄이면서도 작은 깊이의 회로를 제공하기 때문에 2차 Suzuki-Trotter를 선택합니다.

제어된 시간 진화 게이트 제거

확장된 스왑 테스트에서, 단일 앤실라 큐비트로 시간 진화를 제어하려면 앤실라가 시스템 전체에 걸쳐 많은 게이트를 제어해야 합니다. 이는 상당한 라우팅 오버헤드를 초래할 수 있으며, 최악의 경우 사실상 all-to-one 연결성이 필요합니다. 이를 피하기 위해, 해밀토니안의 대칭성을 활용하여 제어된 시간 진화 게이트를 제어가 없는 버전으로 대체함으로써 추가적인 최적화가 가능합니다. 다음 회로를 살펴봅시다.

여기서, $B_{\rm ref}$는 참조 상태를 준비하며, $B_{\rm ref}|0^n\rangle = |\psi_{0}\rangle$입니다.

먼저 $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)|\psi_0\rangle$을 준비한 다음 $|1\rangle$ branch에서만 $U(t)$를 적용하는 대신, 회로는 두 branch를 다음과 같이 직접 준비합니다

$$\begin{equation*} |0\rangle|\psi_0\rangle \quad \text{and} \quad |1\rangle|\psi_d\rangle , \end{equation*}$$

where

$$\begin{equation*} |\psi_d\rangle = U(d\Delta t)|\psi_0\rangle . \end{equation*}$$

회로는 먼저 앤실라에 아다마르 게이트를 적용하고 $|1\rangle$ branch에서만 참조 상태를 준비합니다:

$$\begin{equation*} |0\rangle |0^n\rangle \longrightarrow \frac{|0\rangle |0^n\rangle + |1\rangle B_{\rm ref}|0^n\rangle}{\sqrt{2}} = \frac{|0\rangle |0^n\rangle + |1\rangle |\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

그런 다음 제어되지 않은 시간 진화 연산자가 두 branch 모두에 적용됩니다:

$$\begin{equation*} \frac{|0\rangle |0^n\rangle + |1\rangle |\psi_0\rangle}{\sqrt{2}} \longrightarrow \frac{|0\rangle U(t)|0^n\rangle + |1\rangle U(t)|\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

해밀토니안이 여기 개수를 보존하므로, $|0^n\rangle$가 그 고유 상태임을 알 수 있으며, 따라서 진화 연산자는 해밀토니안 하에서 위상만 축적합니다:

$$\begin{equation*} U(t)|0^n\rangle = e^{-iE_{\rm vac}t}|0^n\rangle . \end{equation*}$$

Therefore,

$$\begin{equation*} \frac{|0\rangle U(t)|0^n\rangle + |1\rangle U(t)|\psi_0\rangle}{\sqrt{2}} = \frac{e^{-iE_{\rm vac}t}|0\rangle |0^n\rangle + |1\rangle U(t)|\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

다음으로, $B_{\rm ref}$는 $|0\rangle$ branch에서만 적용됩니다.

$$\begin{equation*} \frac{e^{-iE_{\rm vac}t}|0\rangle |\psi_0\rangle + |1\rangle U(t)|\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

이 시점에서, 두 branch는 추가적인 상대 위상을 가집니다. 이를 제거하기 위해, 앤실라 위상 게이트를 적용합니다

$$\begin{equation*} P(-E_{\rm vac}t) = \begin{pmatrix} 1 & 0 \\ 0 & e^{-iE_{\rm vac}t} \end{pmatrix}. \end{equation*}$$

이는 상태를 다음과 같이 변환합니다

$$\begin{equation*} \frac{e^{-iE_{\rm vac}t}|0\rangle |\psi_0\rangle + |1\rangle U(t)|\psi_0\rangle}{\sqrt{2}} \longrightarrow \frac{e^{-iE_{\rm vac}t}|0\rangle |\psi_0\rangle + e^{-iE_{\rm vac}t}|1\rangle U(t)|\psi_0\rangle}{\sqrt{2}} . \end{equation*}$$

따라서, 관련 없는 전역 위상까지 감안하면, 마침내 다음을 준비했습니다

$$\begin{equation*} |\Phi_{0d}\rangle = \frac{ |0\rangle|\psi_0\rangle + |1\rangle U(t)|\psi_0\rangle }{\sqrt{2}} . \end{equation*}$$

다음 코드 블록에서, 이 제어가 없는 회로를 구현합니다.

controlled_extended_swap_test = extended_swap_test

# Reuse the ordering and synthesis rule selected by the comparison above.
evol_gate_optimized = PauliEvolutionGate(
    hamiltonian_for_synthesis,
    time=t,
    synthesis=SuzukiTrotter(order=2, reps=num_trotter_steps),
)

uncontrolled_evolution = QuantumCircuit(num_qubits, name="U_ST2(t)")
uncontrolled_evolution.append(evol_gate_optimized, range(num_qubits))

controlled_state_prep = QuantumCircuit(num_qubits + 1, name="C-Prep")
for q, bit in enumerate(reversed(ref_bitstring)):
    if bit == "1":
        controlled_state_prep.cx(ancilla, system_qubits[q])

vacuum_bitstring = "0" * num_qubits
vacuum_energy = basis_state_expectation_sparse(hamiltonian, "0" * num_qubits)

optimized_extended_swap_test = QuantumCircuit(num_qubits + 1)
optimized_extended_swap_test.h(ancilla)

# Prepare |psi_ref> only on the |1> branch.
optimized_extended_swap_test.compose(controlled_state_prep, inplace=True)
optimized_extended_swap_test.barrier()

# Apply the Trotterized time evolution without control.
optimized_extended_swap_test.compose(
    uncontrolled_evolution,
    qubits=system_qubits,
    inplace=True,
)
optimized_extended_swap_test.barrier()

# Map |0>|0...0> to |0>|psi_ref>, leaving the |1> branch unchanged.
optimized_extended_swap_test.x(ancilla)
optimized_extended_swap_test.compose(controlled_state_prep, inplace=True)
optimized_extended_swap_test.x(ancilla)

# Cancel the known vacuum phase so that the same X/Y observables can be used.
optimized_extended_swap_test.p(-vacuum_energy * t, ancilla)
optimized_extended_swap_test = simple_transpilation(
    optimized_extended_swap_test
)

print(f"Vacuum energy E_vac = {vacuum_energy:.1f}")
display(
    optimized_extended_swap_test.assign_parameters({t: 1.0}).draw(
        "mpl", scale=0.5, fold=26
    )
)

Vacuum energy E_vac = 11.0+0.0j

트랜스파일

이제, 하드웨어에서 실행 가능하도록 제어된 회로와 제어가 없는 회로를 트랜스파일합니다. 트랜스파일된 회로의 결과를 비교해 봅시다.

pass_manager = generate_preset_pass_manager(
    backend=backend,
    optimization_level=3,
)

isa_controlled_extended_swap_test = pass_manager.run(
    controlled_extended_swap_test
)

pass_manager = generate_preset_pass_manager(
    backend=backend, optimization_level=3, routing_method="none"
)
isa_optimized_extended_swap_test = pass_manager.run(
    optimized_extended_swap_test
)

transpilation_result = [
    {
        "label": "abstract controlled U(t)",
        **summarize_circuit(isa_controlled_extended_swap_test),
    },
    {
        "label": "optimized non-controlled U(t)",
        **summarize_circuit(isa_optimized_extended_swap_test),
    },
]

display(pd.DataFrame(transpilation_result))

def filter_qubits_from_layout(layout):
    q_layout = Layout(
        {
            physical: virtual
            for physical, virtual in layout.get_physical_bits().items()
            if virtual._register.name == "q"
        }
    )
    return q_layout

print(
    filter_qubits_from_layout(
        isa_optimized_extended_swap_test.layout.initial_layout
    )
)

isa_observables = [
    op.apply_layout(isa_optimized_extended_swap_test.layout)
    for op in observables
]

label  depth   size  2q gates  2q depth
0       abstract controlled U(t)  15457  23786      4686      4580
1  optimized non-controlled U(t)    261   1716       307        57

Layout({
18: <Qubit register=(13, "q"), index=0>,
5: <Qubit register=(13, "q"), index=1>,
6: <Qubit register=(13, "q"), index=2>,
7: <Qubit register=(13, "q"), index=3>,
8: <Qubit register=(13, "q"), index=4>,
9: <Qubit register=(13, "q"), index=5>,
10: <Qubit register=(13, "q"), index=6>,
11: <Qubit register=(13, "q"), index=7>,
12: <Qubit register=(13, "q"), index=8>,
13: <Qubit register=(13, "q"), index=9>,
14: <Qubit register=(13, "q"), index=10>,
15: <Qubit register=(13, "q"), index=11>,
19: <Qubit register=(13, "q"), index=12>
})

Step 3: Qiskit 프리미티브를 사용한 실행

다음 단계는 여러 $d$ 값에 대해 동일한 매개변수화된 회로를 제출하는 것입니다. 각 $d$에 대해, 네 개의 기댓값을 추정합니다: $X\otimes I$, $Y\otimes I$, $X\otimes (H-T)$, $Y\otimes (H-T)$. 이 네 개의 숫자는 그런 다음 복소수 첫 번째 행 요소 $\mathcal{S}_{0d}$와 $\mathcal{H}_{0d}$로 결합됩니다.

여기서, $\mathcal{S}_{00}=1$이고 $\mathcal{H}_{00}=\langle\psi_{0}|(H-T)|\psi_{0}\rangle=\langle\psi_{0}|H|\psi_{0}\rangle-\tau$는 $\lvert\psi_0\rangle$가 희소하므로 고전적으로 계산할 수 있으므로, $d=0$ 경우는 건너뜁니다.

pub_list = []
d_values = list(range(1, krylov_dim))

# Exact local statevector estimator.
estimator = StatevectorEstimator()

# We use the circuit before the transpilation for the local simulator,
# but we will use the transpiled circuit for the real backend.
for d in d_values:
    parameter_values = [d * dt]

    for ob in observables:
        pub_list.append(
            (
                optimized_extended_swap_test,
                ob,
                parameter_values,
            )
        )

job = estimator.run(pub_list)

# Local PrimitiveJob does not provide Runtime-style job inputs,
# so preserve the inputs directly.
inputs = pub_list
result = job.result()

print(f"Number of Krylov basis states: r = {len(d_values)}")
print(f"Number of PUBs: {len(pub_list)}")
print(f"Each d uses observables: {observable_labels}")

Number of Krylov basis states: r = 9
Number of PUBs: 36
Each d uses observables: ['Re S_0d', 'Im S_0d', 'Re shifted H_0d', 'Im shifted H_0d']

4단계: 후처리 및 원하는 고전적 형식으로 결과 반환

투영된 행렬을 추정한 후, GEVP를 정규화하고 풉니다

$$\begin{equation*} \mathcal{H}\mathbf{c} = E\mathcal{S}\mathbf{c}. \end{equation*}$$

가장 작은 일반화된 고유값이 바닥 상태 에너지의 KQD 추정치를 제공합니다.

h_shifted_row_est = np.zeros(krylov_dim, dtype=complex)
s_row_est = np.zeros(krylov_dim, dtype=complex)

h_shifted_row_est[0] = ref_energy - shift_tau
s_row_est[0] = 1.0

for idx, (pub_input, pub_result) in enumerate(zip(inputs, result)):
    d_index, obs_index = divmod(idx, len(observables))
    ev = np.asarray(pub_result.data.evs).reshape(-1)[0]
    std = np.asarray(pub_result.data.stds).reshape(-1)[0]

    if obs_index == 0:
        s_row_est[d_index + 1] = ev
    elif obs_index == 1:
        s_row_est[d_index + 1] += 1j * ev
    elif obs_index == 2:
        h_shifted_row_est[d_index + 1] = ev
    elif obs_index == 3:
        h_shifted_row_est[d_index + 1] += 1j * ev

# H_0d = shifted_H_0d + tau * S_0d.
h_row_est = h_shifted_row_est + shift_tau * s_row_est

h_matrix_est = la.toeplitz(h_row_est.conj(), h_row_est)
s_matrix_est = la.toeplitz(s_row_est.conj(), s_row_est)

s_eigvals = la.eigvalsh(0.5 * (s_matrix_est + s_matrix_est.conj().T))
positive_s_eigvals = s_eigvals[s_eigvals > 1e-12]
s_condition_number = (
    positive_s_eigvals[-1] / positive_s_eigvals[0]
    if len(positive_s_eigvals) > 0
    else np.inf
)

with np.printoptions(precision=3, suppress=True):
    print("Estimated first row of S:")
    print(s_row_est)
    print()
    print("Estimated first row of H:")
    print(h_row_est)
    print()
    print("Eigenvalues of the estimated overlap matrix S:")
    print(s_eigvals)
    print(f"Condition number above 1e-12: {s_condition_number:.3e}")
    print()

Estimated first row of S:
[ 1.   +0.j     0.758-0.596j  0.203-0.836j -0.291-0.636j -0.444-0.229j
 -0.276+0.053j -0.044+0.051j  0.006-0.121j -0.155-0.218j -0.345-0.101j]

Estimated first row of H:
[ 7.   +0.j     4.842-4.76j   0.044-6.185j -3.791-3.653j -4.137+0.39j
 -1.495+2.654j  1.331+1.777j  1.85 -0.763j -0.032-2.278j -2.222-1.372j]

Eigenvalues of the estimated overlap matrix S:
[-0.     0.     0.     0.     0.     0.     0.01   0.355  3.526  6.109]
Condition number above 1e-12: 1.904e+12

이제 회로 추정값에서 재구성된 행렬을 사용하여 일반화된 고유값 문제를 풉니다. 정확한 실시간 진화를 사용하는 이상적인 상태벡터 계산에서는, 이는 정확한 투영된 결과를 재현해야 합니다. 실제로는, 편차가 Trotter화, 샘플링 오차, 그리고 중첩 행렬의 수치적 불안정성으로부터 발생할 수 있습니다.

Krylov 부분 공간의 차원을 증가시킴에 따라 에너지가 어떻게 수렴하는지 관찰합니다.

exact_evals = diagonalize_single_1_subspace(hamiltonian)
exact_ground = min(exact_evals)
print("exact ground state energy:  ", exact_ground)

threshold = 1e-12
energy_convergence = []

for r in range(1, krylov_dim + 1):
    energy_est_kqd, coeffs_est, retained_est = solve_thresholded_gevp(
        h_matrix_est[:r, :r],
        s_matrix_est[:r, :r],
        threshold=threshold,
    )
    energy_convergence.append(energy_est_kqd)
    print(
        f"Krylov ground state energy (dim={r}, retained={retained_est}): ",
        energy_est_kqd,
    )

exact ground state energy:   3.136296694843727
Krylov ground state energy (dim=1, retained=1):  7.0
Krylov ground state energy (dim=2, retained=2):  4.184510657551266
Krylov ground state energy (dim=3, retained=3):  3.5539074630394136
Krylov ground state energy (dim=4, retained=4):  3.3366270761341044
Krylov ground state energy (dim=5, retained=5):  3.252017453225087
Krylov ground state energy (dim=6, retained=6):  3.2300275138879186
Krylov ground state energy (dim=7, retained=7):  3.2299154099085685
Krylov ground state energy (dim=8, retained=7):  3.2298063744216776
Krylov ground state energy (dim=9, retained=7):  3.2296778282872456
Krylov ground state energy (dim=10, retained=8):  3.223647515867734

def plot_energy_convergence(energy_convergence, exact_ground, krylov_dim):
    fig, ax = plt.subplots(figsize=(7, 4.5))

    ax.plot(
        range(1, krylov_dim + 1),
        energy_convergence,
        marker="o",
        label="KQD estimate",
    )

    ax.axhline(
        exact_ground,
        linestyle="--",
        label=f"Exact ground energy = {exact_ground:.6f}",
    )

    ax.set_xlabel("Krylov dimension")
    ax.set_ylabel("Ground-state energy")
    ax.grid(True, alpha=0.3)
    ax.legend()
    plt.tight_layout()
    plt.show()

plot_energy_convergence(energy_convergence, exact_ground, krylov_dim)

대규모 하드웨어 예제

이전 섹션에서는 상태벡터 시뮬레이션을 진단으로 사용할 수 있도록 12-큐비트 모델을 사용했습니다. 이제 동일한 KQD 워크플로를 30-큐비트 하이젠베르크 체인으로 확장하고 IBM Quantum 하드웨어에서 실행하기 위한 워크로드를 준비합니다.

단계 1-4를 하나의 코드 블록으로 압축

이제 이러한 모든 세부 사항을 더 큰 규모의 단일 워크플로로 통합하여 실제 양자 하드웨어에서 실행합니다. 이 섹션에서는, 결과의 신뢰성을 향상시키기 위해 현실적인 오류 완화 설정을 적용합니다. 서로 다른 $d$ 값에 해당하는 행렬 요소는 병렬로 평가될 수 있으므로, 이를 효율적으로 실행하기 위해 Batch 모드를 사용합니다.

# -------------------------Step 1-------------------------

# Map the classical problem to quantum circuits and observables.
# Problem and KQD parameters.
large_num_qubits = 30
large_krylov_dim = 7
large_num_trotter_steps = 3
large_dt = np.pi / (3 * (large_num_qubits - 1))
large_t = Parameter("t_large")

# Use a single excitation near the center of the chain.
large_excitation_qubit = large_num_qubits // 2
large_ref_label = ["0"] * large_num_qubits
large_ref_label[large_num_qubits - 1 - large_excitation_qubit] = "1"
large_ref_bitstring = "".join(large_ref_label)

# Use the ordering and product formula selected in the preceding section.
large_hamiltonian = make_heisenberg_hamiltonian_ordered(
    large_num_qubits,
    ordering="even-odd edge-grouped",
    coupling=1.0,
)
large_ref_energy = float(
    np.real(
        basis_state_expectation_sparse(large_hamiltonian, large_ref_bitstring)
    )
)
large_vacuum_energy = float(
    np.real(
        basis_state_expectation_sparse(
            large_hamiltonian, "0" * large_num_qubits
        )
    )
)

large_evolution_gate = PauliEvolutionGate(
    large_hamiltonian,
    time=large_t,
    synthesis=SuzukiTrotter(order=2, reps=large_num_trotter_steps),
)
large_uncontrolled_evolution = QuantumCircuit(
    large_num_qubits,
    name="U_ST2_large(t)",
)
large_uncontrolled_evolution.append(
    large_evolution_gate,
    range(large_num_qubits),
)

# Build the control-free extended-swap-test circuit.
large_ancilla = 0
large_system_qubits = list(range(1, large_num_qubits + 1))

large_controlled_state_prep = QuantumCircuit(
    large_num_qubits + 1,
    name="C-Prep-large",
)
large_controlled_state_prep.cx(
    large_ancilla,
    large_system_qubits[large_excitation_qubit],
)

large_extended_swap_test = QuantumCircuit(large_num_qubits + 1)
large_extended_swap_test.h(large_ancilla)
large_extended_swap_test.compose(large_controlled_state_prep, inplace=True)
large_extended_swap_test.compose(
    large_uncontrolled_evolution,
    qubits=large_system_qubits,
    inplace=True,
)
large_extended_swap_test.x(large_ancilla)
large_extended_swap_test.compose(
    large_controlled_state_prep.inverse(),
    inplace=True,
)
large_extended_swap_test.x(large_ancilla)
large_extended_swap_test.p(
    -large_vacuum_energy * large_t,
    large_ancilla,
)

# Reuse the Hamiltonian-shifting construction from the preceding section.
(
    large_obs_x_shifted_hamiltonian,
    large_obs_y_shifted_hamiltonian,
    large_shift_tau,
) = make_reduced_heisenberg_observables(large_ref_bitstring)

large_observables = [
    SparsePauliOp("I" * large_num_qubits + "X"),
    SparsePauliOp("I" * large_num_qubits + "Y"),
    large_obs_x_shifted_hamiltonian,
    large_obs_y_shifted_hamiltonian,
]
large_observable_labels = [
    "Re S_0d",
    "Im S_0d",
    "Re shifted H_0d",
    "Im shifted H_0d",
]

# -------------------------Step 2-------------------------

# Optimize the problem for quantum execution.
# Select a real backend and transpile the parameterized circuit to ISA form.
large_backend = service.backend("ibm_boston")
large_pass_manager = generate_preset_pass_manager(
    backend=large_backend, optimization_level=3, routing_method="none"
)
large_isa_circuit = large_pass_manager.run(large_extended_swap_test)
large_isa_observables = [
    observable.apply_layout(large_isa_circuit.layout)
    for observable in large_observables
]

large_two_qubit_gate_count = sum(
    instruction.operation.num_qubits == 2
    for instruction in large_isa_circuit.data
)

print(f"Backend: {large_backend.name}")
print(f"System qubits: {large_num_qubits}")
print(f"Total circuit qubits: {large_isa_circuit.num_qubits}")
print(f"Krylov dimension: {large_krylov_dim}")
print(f"Time step: {large_dt:.6f}")
print(f"Shift tau: {large_shift_tau}")
print(f"ISA circuit depth: {large_isa_circuit.depth()}")
print(f"ISA two-qubit gates: {large_two_qubit_gate_count}")
print(
    f"ISA two-qubit depth: {large_isa_circuit.depth(lambda x: x[0].num_qubits == 2)}"
)

# -------------------------Step 3-------------------------

# Execute on quantum hardware with Qiskit Runtime primitives.
# Submit one job per d, with all four observables in that job.
large_d_values = list(range(1, large_krylov_dim))
retrieve_batch_id = None
large_jobs = []

if retrieve_batch_id is None:
    large_estimator_options = {
        "default_shots": 8192,
        "dynamical_decoupling": {
            "enable": True,
            "sequence_type": "XpXm",
        },
        "resilience": {
            "measure_mitigation": True,
            "measure_noise_learning": {
                "num_randomizations": 32,
                "shots_per_randomization": 256,
            },
            "layer_noise_learning": {
                "max_layers_to_learn": 4,
                "layer_pair_depths": [0, 1, 2, 4, 16, 32],
                "num_randomizations": 32,
                "shots_per_randomization": 128,
            },
            "zne_mitigation": True,
            "zne": {
                "amplifier": "pea",
                "noise_factors": [1.0, 1.5, 2.0],
                "extrapolator": ("exponential", "linear"),
            },
        },
        "twirling": {
            "enable_gates": True,
            "enable_measure": True,
            "num_randomizations": 32,
            "shots_per_randomization": 256,
            "strategy": "active-accum",
        },
    }

    with Batch(backend=large_backend) as large_batch:
        large_batch_id = large_batch.session_id
        large_estimator = EstimatorV2(
            mode=large_batch,
            options=large_estimator_options,
        )
        # Krylov quantum diagonalization of lattice Hamiltonians -> TUT_KQDOLH.
        large_estimator.options.environment.job_tags = ["TUT_KQDOLH"]

        for d in large_d_values:
            parameter_values = [d * large_dt]
            pubs_for_d = [
                (large_isa_circuit, observable, parameter_values)
                for observable in large_isa_observables
            ]
            job = large_estimator.run(pubs_for_d)
            large_jobs.append(job)
            print(
                f"Submitted d={d}: job_id={job.job_id()}, "
                f"PUBs={len(pubs_for_d)}"
            )

    print(f"Batch ID: {large_batch_id}")
else:
    large_batch_id = retrieve_batch_id
    large_jobs = service.jobs(
        session_id=large_batch_id,
        limit=None,
        descending=False,
    )

large_job_ids_by_d = {
    d: job.job_id() for d, job in zip(large_d_values, large_jobs)
}
print(f"Job IDs by d: {large_job_ids_by_d}")

# -------------------------Step 4-------------------------

# Post-process the quantum results and solve the classical GEVP.
# Reconstruct the first rows of S and the shifted Hamiltonian matrix.
large_s_row_est = np.zeros(large_krylov_dim, dtype=complex)
large_h_shifted_row_est = np.zeros(large_krylov_dim, dtype=complex)
large_s_row_est[0] = 1.0
large_h_shifted_row_est[0] = large_ref_energy - large_shift_tau

for d, job in zip(large_d_values, large_jobs):
    job_result = job.result()

    if len(job_result) != len(large_observables):
        raise RuntimeError(
            f"Expected {len(large_observables)} PUB results for d={d}, "
            f"but received {len(job_result)}."
        )

    expectation_values = [
        np.asarray(pub_result.data.evs).reshape(-1)[0]
        for pub_result in job_result
    ]

    large_s_row_est[d] = expectation_values[0] + 1j * expectation_values[1]
    large_h_shifted_row_est[d] = (
        expectation_values[2] + 1j * expectation_values[3]
    )

# H_0d = shifted_H_0d + tau * S_0d.
large_h_row_est = large_h_shifted_row_est + large_shift_tau * large_s_row_est

large_s_matrix_est = la.toeplitz(
    large_s_row_est.conj(),
    large_s_row_est,
)
large_h_matrix_est = la.toeplitz(
    large_h_row_est.conj(),
    large_h_row_est,
)

large_exact_gnd = min(diagonalize_single_1_subspace(large_hamiltonian))

large_energy_convergence = []

with np.printoptions(precision=5, suppress=True):
    print("Estimated first row of S:")
    print(large_s_row_est)
    print("Estimated first row of H:")
    print(large_h_row_est)
    print("exact ground state energy:  ", large_exact_gnd)

    for r in range(1, large_krylov_dim + 1):
        energy_est_kqd, coeffs_est, retained_est = solve_thresholded_gevp(
            large_h_matrix_est[:r, :r],
            large_s_matrix_est[:r, :r],
            threshold=5e-2,
        )
        large_energy_convergence.append(energy_est_kqd)
        print(
            f"Krylov ground state energy (dim={r}, retained={retained_est}): ",
            energy_est_kqd,
        )

plot_energy_convergence(
    large_energy_convergence, large_exact_gnd, large_krylov_dim
)

Backend: ibm_boston
System qubits: 30
Total circuit qubits: 156
Krylov dimension: 7
Time step: 0.036110
Shift tau: 25.0
ISA circuit depth: 219
ISA two-qubit gates: 728
ISA two-qubit depth: 52
Job IDs by d: {1: 'd9fl4p4jeosc73fk4dg0', 2: 'd9fl4pineu4c739poecg', 3: 'd9fl4q2neu4c739poedg', 4: 'd9fl4qhhtsac739fjgn0', 5: 'd9fl4r4jeosc73fk4dhg', 6: 'd9fl4rkjeosc73fk4dj0'}
Estimated first row of S:
[ 1.     +0.j       0.42055-0.65466j -0.17883-0.73125j -0.64523-0.31601j
 -0.70781+0.55574j -0.02467+0.70525j  0.59301+0.55602j]
Estimated first row of H:
[ 25.      +0.j       10.3419 -16.5586j   -5.017  -18.03933j
 -16.44797 -6.98173j -17.17081+14.84502j   0.62722+17.81196j
  15.68886+12.74875j]
exact ground state energy:   21.021912418526902
Krylov ground state energy (dim=1, retained=1):  25.0
Krylov ground state energy (dim=2, retained=2):  24.432409110686205
Krylov ground state energy (dim=3, retained=3):  24.256856893278556
Krylov ground state energy (dim=4, retained=4):  23.727409826799715
Krylov ground state energy (dim=5, retained=4):  23.324720470780864
Krylov ground state energy (dim=6, retained=5):  21.91957579005085
Krylov ground state energy (dim=7, retained=5):  21.548331214122644

부록: 해밀토니안 함수(스펙트럼 필터) 관점

주요 워크플로는 KQD를 조작적으로 제시했습니다: 실시간 진화된 상태로부터 Krylov 기저를 구축하고, 투영된 행렬 $\mathcal{H}$와 $\mathcal{S}$를 추정하며, GEVP를 풉니다. 이 부록은 KQD가 왜 작동하는지 설명하는 보완적인 관점에서 동일한 계산을 다시 살펴봅니다: 해밀토니안 함수, 또는 스펙트럼 필터 관점 [3], [5]입니다. 위에서 이미 얻은 12-큐비트 모델, 시간 스텝 $\Delta t$, 그리고 Krylov 해를 재사용하며, 새로운 회로 실행이 필요하지 않습니다.

에너지 분포로서의 참조 상태

해밀토니안이 다음과 같은 고유분해를 갖는다고 합시다

$$\begin{equation*} H=\sum_{m} E_m\,|E_m\rangle\!\langle E_m|, \qquad E_0\le E_1\le\cdots, \end{equation*}$$

여기서 $|E_m\rangle$은 에너지 고유 상태입니다. 임의의 참조 상태는 이 고유 기저에서 전개될 수 있습니다,

$$\begin{equation*} |\psi_{0}\rangle=\sum_m a_m\,|E_m\rangle, \qquad a_m=\langle E_m|\psi_{0}\rangle, \end{equation*}$$

따라서 각 에너지 $E_m$에서 스펙트럼 가중치 $p_m=|a_m|^2$를 가집니다. 참조 에너지는 이 분포의 평균, $\langle H\rangle_{0}=\sum_m p_m E_m$입니다.

일반적인 $n$-큐비트 해밀토니안의 고유분해는 지수적으로 비용이 많이 들므로, 이 그림은 진단일 뿐이며, 알고리즘의 일부는 아닙니다. 하지만 여기서는, 동일한 $n=12$ 문제에 대해 저렴하게 계산할 수 있습니다: 하이젠베르크 해밀토니안은 총 여기 개수를 보존하며, 참조 $|000001000000\rangle$는 단일 여기를 가지므로, 그 전체 스펙트럼 내용은 차원이 $n$에 선형적으로만 증가하는 단일 여기 부분 공간에 존재합니다. 따라서 우리는 정확한 단일 여기 블록(위에서 이미 벤치마크로 사용됨)을 재사용하고 그 부분 공간 내의 참조 분포 $\{(E_m, p_m)\}$를 읽어냅니다.

# Exact single-excitation-subspace decomposition of the reference state.
# This reuses make_heisenberg_hamiltonian / _basis_state_transition_amplitude_sparse
# and the n=12 `hamiltonian` and `ref_bitstring` defined in the small-scale example.
single_excitation_states = [1 << k for k in range(num_qubits)]

h_single = np.array(
    [
        [
            _basis_state_transition_amplitude_sparse(hamiltonian, bra, ket)
            for ket in single_excitation_states
        ]
        for bra in single_excitation_states
    ]
)
h_single = 0.5 * (h_single + h_single.conj().T)

subspace_evals, subspace_evecs = la.eigh(h_single)
subspace_evals = np.real(subspace_evals)

# Reference-state coordinates inside the single-excitation subspace.
ref_position = single_excitation_states.index(int(ref_bitstring, 2))
ref_in_subspace = np.zeros(num_qubits, dtype=complex)
ref_in_subspace[ref_position] = 1.0

# Amplitudes and spectral weights of the reference in the energy eigenbasis.
ref_eigen_amplitudes = subspace_evecs.conj().T @ ref_in_subspace
ref_spectral_weights = np.abs(ref_eigen_amplitudes) ** 2

print(f"Single-excitation subspace dimension: {num_qubits}")
print(f"Subspace ground-state energy:         {subspace_evals[0]:.6f}")
print(
    f"Reference energy (sum p_m E_m):       {np.sum(ref_spectral_weights * subspace_evals):.6f}"
)
print(f"Reference weight on subspace ground:  {ref_spectral_weights[0]:.6f}")

Single-excitation subspace dimension: 12
Subspace ground-state energy:         3.136297
Reference energy (sum p_m E_m):       7.000000
Reference weight on subspace ground:  0.163827

KQD는 이 분포를 재구성하는 필터를 학습합니다

해밀토니안 함수 $f(H)$는 스펙트럼 미적분을 통해 정의됩니다,

$$\begin{equation*} f(H)=\sum_m f(E_m)\,|E_m\rangle\!\langle E_m|, \end{equation*}$$

다시 말해, 고유 투영자들의 가중합입니다. 이를 참조에 적용하면 각 스펙트럼 진폭이 재구성됩니다, $a_m \to a_m f(E_m)$:

$$\begin{equation*} f(H)\,|\psi_{0}\rangle=\sum_m a_m\,f(E_m)\,|E_m\rangle . \end{equation*}$$

만약 $f$가 가장 낮은 에너지에서 날카롭게 피크를 이루면($f(E_0)=1$이고 그 외에는 $f(E_m)\approx 0$), $f(H)$는 바닥 상태 투영자로 작용하며, 정규화된 출력은 바닥 상태에 (거의) 근접합니다. 따라서, 에너지에서 좋은 저역 통과 스펙트럼 필터가 정확히 우리가 원하는 것입니다.

KQD는 $f$를 미리 규정하지 않습니다. 대신 실시간 진화 기저에서 필터를 전개합니다,

$$\begin{equation*} f_{\rm KQD}(E)=\sum_{\ell=0}^{r-1} c_\ell\, e^{-i\ell\Delta t E}, \end{equation*}$$

에너지의 삼각함수이며, 그 계수 $\{c_\ell\}$는 정확히 위에서 풀린 GEVP 고유 벡터입니다. 따라서 레일리 몫(Rayleigh quotient) $\mathbf{c}^\dagger\mathcal{H}\mathbf{c}/\mathbf{c}^\dagger\mathcal{S}\mathbf{c}$를 최소화하는 것은 참조의 여기 상태 가중치를 가장 잘 억제하는 필터를 학습하는 것과 동일합니다. 더 큰 Krylov 차원 $r$은 필터에 더 많은 자유도와 바닥 상태 에너지에서 더 날카로운 피크를 제공합니다.

아래의 헬퍼는 이 학습된 필터를 에너지 축에서 평가합니다; 그런 다음 위에서 얻은 참조 분포에 적용합니다.

def trigonometric_krylov_filter(
    coeffs: np.ndarray,
    energies: np.ndarray,
    time_step: float,
) -> np.ndarray:
    """Evaluate the learned Krylov filter f(E) = sum_l c_l exp(-i l dt E)."""
    values = np.zeros_like(energies, dtype=complex)
    for ell, coeff in enumerate(coeffs):
        values += coeff * np.exp(-1j * ell * time_step * energies)
    return values

def filtered_spectral_weights(
    weights: np.ndarray,
    filter_values: np.ndarray,
) -> np.ndarray:
    """Reshape spectral weights by |f(E)|^2 and renormalize."""
    reshaped = weights * np.abs(filter_values) ** 2
    return reshaped / np.sum(reshaped)

# Recover the KQD coefficients from the already-estimated projected matrices.
# The shift only moves H by tau * S, so it does not change the GEVP eigenvector;
# we solve at the full Krylov dimension used in the small-scale example.
_, kqd_coeffs, _ = solve_thresholded_gevp(
    h_matrix_est,
    s_matrix_est,
    threshold=1e-12,
)

filter_on_spectrum = trigonometric_krylov_filter(
    kqd_coeffs, subspace_evals, dt
)
filtered_weights = filtered_spectral_weights(
    ref_spectral_weights, filter_on_spectrum
)

print(f"Ground-state overlap  (reference): {ref_spectral_weights[0]:.4f}")
print(f"Ground-state overlap  (filtered):  {filtered_weights[0]:.4f}")
print(
    f"Mean energy           (reference): {np.sum(ref_spectral_weights * subspace_evals):.6f}"
)
print(
    f"Mean energy           (filtered):  {np.sum(filtered_weights * subspace_evals):.6f}"
)

Ground-state overlap  (reference): 0.1638
Ground-state overlap  (filtered):  0.9638
Mean energy           (reference): 7.000000
Mean energy           (filtered):  3.164503

필터와 그 유연성 시각화

먼저 위에서 사용된 전체 Krylov 차원에서 학습된 필터를 보여준 다음, 차원 $r$이 증가함에 따라 어떻게 날카로워지는지 추적합니다.

막대는 참조 스펙트럼 가중치 $p_m$(이전)과 필터링된 가중치 $p_m|f_{\rm KQD}(E_m)|^2$(이후)를 보여주며, 연속적인 에너지 축에서 학습된 필터 강도 $|f_{\rm KQD}(E)|^2$와 함께 표시됩니다. 필터는 단일 여기 부분 공간의 가장 낮은 에너지에 가중치를 집중시킵니다 — 이는 소규모 예제에서 KQD 추정치가 수렴한 것과 동일한 에너지입니다. 이는 단일 여기 섹터 내에서의 바닥 상태이며, 이 여기 보존 참조 상태에 대한 관련 목표이지, 전역 바닥 상태가 아님에 유의하십시오.

fig, ax = plt.subplots(figsize=(8, 4))

visible = (ref_spectral_weights > 1e-4) | (filtered_weights > 1e-4)
ax.bar(
    subspace_evals[visible],
    ref_spectral_weights[visible],
    width=0.18,
    alpha=0.45,
    label="reference $p_m$",
)
ax.bar(
    subspace_evals[visible],
    filtered_weights[visible],
    width=0.14,
    alpha=0.9,
    label=r"filtered $p_m\,|f_{\rm KQD}(E_m)|^2$",
)

energy_grid = np.linspace(
    subspace_evals.min() - 0.5, subspace_evals.max() + 0.5, 800
)
filter_intensity = (
    np.abs(trigonometric_krylov_filter(kqd_coeffs, energy_grid, dt)) ** 2
)
filter_intensity /= filter_intensity.max()
ax.plot(
    energy_grid,
    filter_intensity,
    color="k",
    linewidth=2,
    label=r"$|f_{\rm KQD}(E)|^2$ (normalized)",
)

ax.axvline(
    subspace_evals[0],
    color="C3",
    linestyle="--",
    linewidth=1,
    label="subspace ground energy",
)
ax.set_xlabel("Energy eigenvalue $E_m$")
ax.set_ylabel("Spectral weight")
ax.set_ylim(0, 1)
ax.legend(loc="upper right")
plt.tight_layout()
plt.show()

Krylov 차원 증가: 학습된 함수의 유연성

학습된 필터는 $r$개의 계수를 가진 에너지에 대한 삼각 다항식임을 상기하십시오,

$$\begin{equation*} f_{\rm KQD}(E)=\sum_{\ell=0}^{r-1} c_\ell\, e^{-i\ell\Delta t E}. \end{equation*}$$

Krylov 차원 $r$은 정확히 자유 계수의 개수이므로, 함수의 유연성을 제어합니다. 작은 $r$은 낮은 여기 상태로 가중치가 새어 나가는 넓고 완만하게 변하는 필터만 생성할 수 있습니다; $r$이 증가하면, 필터는 목표 에너지에서 더 좁은 피크를 형성하고 나머지 여기 상태 가중치를 더 적극적으로 억제할 수 있습니다. 이는 소규모 예제에서 관찰된 에너지 수렴에 대응하는 스펙트럼 필터의 대응물입니다: $r$이 증가함에 따라, 필터링된 분포는 부분 공간의 바닥 상태로 붕괴되고 추정된 에너지는 그것으로 감소합니다.

위에서 이미 추정한 투영된 행렬을 재사용하고, 각 선행 $r\times r$ 블록에서 GEVP를 단순히 풀어, 해당 필터를 평가하고 플롯합니다.

# Sweep the Krylov dimension using the leading r x r blocks of the estimated matrices.
sweep_dims = [r for r in (2, 4, 6, 8, krylov_dim) if r <= krylov_dim]
sweep_dims = sorted(set(sweep_dims))

energy_grid = np.linspace(
    subspace_evals.min() - 0.5, subspace_evals.max() + 0.5, 800
)

sweep_cases = []
print(" r    retained    ground overlap    filtered energy")
print("--    --------    --------------    ---------------")
for r in sweep_dims:
    _, coeffs_r, retained_r = solve_thresholded_gevp(
        h_matrix_est[:r, :r],
        s_matrix_est[:r, :r],
        threshold=1e-12,
    )
    filter_on_spectrum_r = trigonometric_krylov_filter(
        coeffs_r, subspace_evals, dt
    )
    filtered_weights_r = filtered_spectral_weights(
        ref_spectral_weights, filter_on_spectrum_r
    )
    filtered_energy_r = float(np.sum(filtered_weights_r * subspace_evals))

    sweep_cases.append((r, coeffs_r, filtered_weights_r))
    print(
        f"{r:2d}    {retained_r:8d}    {filtered_weights_r[0]:14.4f}    {filtered_energy_r:15.6f}"
    )

print(f"\nSubspace ground-state energy (target): {subspace_evals[0]:.6f}")

# One panel per Krylov dimension: filtered spectrum (bars) + filter intensity (curve).
fig, axes = plt.subplots(
    len(sweep_cases),
    1,
    figsize=(8, 2.1 * len(sweep_cases)),
    sharex=True,
)
axes = np.atleast_1d(axes)

for idx, (ax, (r, coeffs_r, filtered_weights_r)) in enumerate(
    zip(axes, sweep_cases)
):
    visible = (ref_spectral_weights > 1e-4) | (filtered_weights_r > 1e-4)
    ax.bar(
        subspace_evals[visible],
        ref_spectral_weights[visible],
        width=0.18,
        alpha=0.35,
        color="C0",
        label="reference $p_m$" if idx == 0 else None,
    )
    ax.bar(
        subspace_evals[visible],
        filtered_weights_r[visible],
        width=0.14,
        alpha=0.9,
        color="C1",
        label="filtered weights" if idx == 0 else None,
    )

    filter_intensity_r = (
        np.abs(trigonometric_krylov_filter(coeffs_r, energy_grid, dt)) ** 2
    )
    filter_intensity_r /= filter_intensity_r.max()
    ax.plot(energy_grid, filter_intensity_r, color="k", linewidth=2)

    ax.axvline(subspace_evals[0], color="C3", linestyle="--", linewidth=1)
    ax.set_ylim(0, 1)
    ax.set_ylabel("weight")
    ax.legend(loc="upper right", title=f"$r={r}$")

axes[-1].set_xlabel("Energy eigenvalue $E_m$")
fig.suptitle(
    r"KQD-learned filter $|f_{\rm KQD}(E)|^2$ sharpening with Krylov dimension $r$",
    y=1.0,
)
fig.tight_layout()
plt.show()

r    retained    ground overlap    filtered energy
--    --------    --------------    ---------------
 2           2            0.4549           4.184510
 4           4            0.8373           3.310168
 6           6            0.9706           3.158100
 8           7            0.9701           3.158725
10           8            0.9638           3.164503

Subspace ground-state energy (target): 3.136297

다음 단계

이 작업이 흥미로웠다면, 다음 자료에 관심이 있을 수 있습니다:

추천

• 페르미온 격자 모델의 샘플 기반 Krylov 양자 대각화
• VQE를 이용한 하이젠베르크 체인의 바닥 상태 에너지 추정
• 양자 대각화 알고리즘 코스

참고 문헌

[1] E. N. Epperly, L. Lin, and Y. Nakatsukasa, A theory of quantum subspace diagonalization, SIAM Journal on Matrix Analysis and Applications 43, 1263-1290 (2022).

[2] N. Yoshioka, M. Amico, W. Kirby, et al., Diagonalization of large many-body Hamiltonians on a quantum processor, arXiv:2407.14431 (2024).

[3] R. M. Parrish and P. L. McMahon, Quantum filter diagonalization: quantum eigendecomposition without full quantum phase estimation, Physical Review Letters 122, 230401 (2019).

[4] G. Lee, S. Choi, J. Huh, and A. F. Izmaylov, Efficient strategies for reducing sampling error in quantum Krylov subspace diagonalization, Digital Discovery 4, 954-969 (2025).

[5] G. Lee, M. Kang, J. Hong, S. Fomichev and J. Huh, Filtered Quantum Phase Estimation, arXiv:2510.04294 (2025).