Krylov quantum diagonalization ng lattice Hamiltonians

Tantiya sa paggamit: 70 minuto sa Heron o Nighthawk processor (PAALALA: Ito ay tantiya lamang. Maaaring mag-iba ang iyong runtime.)

Mga resulta ng pagkatuto

• Paano bigyang-kahulugan ang Krylov quantum diagonalization (KQD) bilang pag-aaral ng isang finite Hamiltonian function na kumikilos bilang isang spectral filter.

• Paano bumuo ng projected Hamiltonian at overlap matrices gamit ang extended swap-test measurements.

• Paano lutasin ang resultang generalized eigenvalue problem (GEVP) at kunin ang isang ground-state energy estimate para sa isang lattice Hamiltonian.

Prerequisites

• Qiskit primitives

• Pagtatantya ng ground-state energy ng Heisenberg chain gamit ang VQE

• Kurso sa quantum diagonalization algorithms

Background

Ang tutorial na ito ay nagpapakita kung paano ipatupad ang Krylov quantum diagonalization (KQD) algorithm sa konteksto ng Qiskit patterns. Matututunan mo muna ang teorya sa likod ng algorithm at pagkatapos ay makikita mo ang demonstrasyon ng pagpapatupad nito sa QPU.

Ang pagtantiya sa mga low-energy property ng many-body Hamiltonians ay isang pangunahing gawain sa quantum simulation. Halimbawa, ang mga ground-state energy at low-lying excitation ay direktang kaugnay sa chemical stability, magnetic ordering, quantum phase transitions, at material response. Sa isang classical computer, lumalaki nang exponential ang Hilbert-space dimension kasabay ng bilang ng orbitals o spins, kaya't mabilis na nagiging hindi praktikal ang direktang diagonalization.

May ilang quantum-computing na paraan para sa problemang ito. Ginagamit ng mga near-term variational method, tulad ng variational quantum eigensolver (VQE), ang medyo shallow na parameterized circuits, ngunit kailangan nila ng isang nonlinear classical optimization loop na may maraming quantum-circuit evaluation. Sa kabilang banda, nagbibigay ang quantum phase estimation (QPE) ng mas direktang daan patungo sa eigenvalue estimation na may mahigpit na garantiya, ngunit nangangailangan ang standard QPE ng mahahabang coherent circuits at pangunahing angkop para sa fault-tolerant na quantum computers. Nasa gitna ng dalawang paraang ito ang KQD: ginagamit nito ang real-time Hamiltonian evolution, tulad ng sa mga phase-estimation-based na algorithm, ngunit pinapalitan ang full phase estimation ng isang compact projected eigenvalue problem na malulutas nang classical.

Isaalang-alang ang isang $n$-qubit Hamiltonian $H$ at isang reference state $\lvert \psi_{0}\rangle$. Bumubuo ang KQD method ng isang Krylov subspace mula sa real-time evolved states,

$$\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*}$$

kung saan ang $r$ ay ang Krylov dimension at ang $\Delta t$ ay ang time step. Anumang state sa Krylov subspace ay ipinapakita bilang isang linear combination ng mga basis states na ito,

$$\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.

Sa simpleng algebra, makikita natin na ang katumbas na energy ay isinusulat bilang Rayleigh quotient,

$$\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*}$$

Dito, ang mga matrices $\mathcal{S}$ at $\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*}$$

ay tumutukoy sa projected overlap at Hamiltonian matrices. Ang kanilang mga entry ay tinatantiya gamit ang quantum-circuit measurements.

Layunin nating hanapin ang coefficient $\mathbf{c}$ na nagbibigay ng minimum $E(\mathbf{c})$:

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

Ayon sa Rayleigh-Ritz theorem, katumbas ng minimization na ito ang paglutas ng generalized eigenvalue problem (GEVP),

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

Tandaan na ang dimension $r$ ay maaaring maliit nang sapat para lutasin ng isang classical computer ang GEVP.

Ito ang parehong variational principle na ginagamit sa classical subspace diagonalization, ngunit dito, nabubuo ang mga basis state sa pamamagitan ng quantum time evolution. Kung ikukumpara sa VQE, karaniwang nangangailangan ang KQD ng mas malalim na circuits dahil umaasa ito sa real-time evolution. Bilang kapalit, iniiwasan ng KQD ang nonlinear parameter optimization at iterative quantum-hardware execution, at systematic itong bumubuti habang pinapalawak ang projected subspace. Naipakita na ang algorithm sa malaking sukat sa umiiral na quantum hardware [2], at maaaring suriin ang performance nito na may mapapatunayang garantiya [1].

Requirements

Bago simulan ang tutorial na ito, siguraduhing mayroon kayo ng sumusunod na naka-install:

• Qiskit SDK v2.3 o mas bago na may suporta sa visualization

• Qiskit Runtime v0.22 o mas bago ( pip install qiskit-ibm-runtime )

• SciPy (pip install scipy)

• Matplotlib (pip install matplotlib)

• Pandas (pip install pandas)

Ang hardware execution ay nangangailangan ng qiskit-ibm-runtime at access sa isang IBM Quantum® account.

Setup

Ini-import ng setup cell ang mga kinakailangang module at nagdedepine ng mga helper function para sa workflow:

1. bumuo ng Heisenberg Hamiltonian;

2. lutasin ang thresholded GEVP;

3. suriin ang learned Krylov filter;

4. i-convert ang mga filter value sa spectral weights;

5. i-plot ang reference at filtered energy distributions.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "pandas": "pandas", "qiskit": "qiskit", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

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))

Sa unang bahagi ng tutorial na ito, ipinapakita natin ang KQD method sa pamamagitan ng paggamit ng local statevector simulator. Pagkatapos, gagamit tayo ng tunay na quantum backend para tugunan ang isang utility-scale na problema.

Nagdedepine rin tayo ng fake backend upang ipakita ang backend-specific na transpilation at suriin ang resultang circuit.

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

backend = FakeMarrakesh()

Small-scale na halimbawa gamit ang simulator

Hakbang 1: I-map ang mga classical input sa isang quantum na problema

Hamiltonian at reference state

Ginagamit ng halimbawang ito ang isang 12-qubit open-boundary Heisenberg chain ($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*}$$

na may single-excitation product state

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

bilang reference state. Dahil ang Heisenberg Hamiltonian na ipinakita sa itaas ay nagpapanatili ng kabuuang bilang ng excitations, nananatili ang reference state sa single-excitation subspace, na ang dimension ay lumalaki lamang nang linear kasabay ng bilang ng qubits. Kaya naman, magagawa nating i-compute nang epektibo ang eksaktong ground-state energy, sa pamamagitan ng pag-diagonalize ng Hamiltonian na nilimitahan sa subspace na iyon, at gamitin ito bilang isang diagnostic benchmark lamang para sa KQD estimate. Tinatantiya mismo ng KQD workflow ang mga projected matrix element gamit ang Qiskit primitives at lulutasin nang classical ang resultang projected problem.

# 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)

Itakda ang mga parameter para sa algorithm

Batay sa mga upper-bound sa Hamiltonian norm, iminumungkahi nang heuristic ng Ref. [1] ang time step $\Delta t$ bilang $\pi/\|H\|$. Dahil mahirap i-compute ang spectral norm $\|H\|$, ginagamit natin sa halip ang upper bound nito:

$$\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*}$$

Itinatakda natin ang Krylov dimension sa $r=10$ at ang bilang ng Trotter steps kada time step sa $5$: isang Krylov space na malaki nang sapat upang malutas ang low-lying spectrum habang pinapanatiling abot-kaya ang pinakamalalim na circuit ($t_{\max}=(r-1)\Delta t$), at sapat na bilang ng Trotter steps upang panatilihing maliit ang discretization error sa pinakamalalim na circuit na iyon.

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

Bumuo ng circuit

Dito, bubuo tayo ng mga circuit upang tantiyahin ang mga matrix element na $\mathcal{H}_{k\ell}$ at $\mathcal{S}_{k\ell}$. Dahil nag-commute ang lahat ng powers ng $H$, mayroon tayong

$$\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*}$$

Such matrices where the elements depend on the index differences in this manner are called Toeplitz and can be reconstructed from first-row elements indexed by $d=\ell-k$.

Dito, ipinapakita natin ang circuit, tinatawag na extended-swap-test, na naghahanda ng

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

kung saan ang $|\psi_d\rangle=e^{-id\Delta tH}|\psi_{0}\rangle$.

Reference state

Inihahanda natin ang reference state $|\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))

Time evolution

Isakatuparan ang time-evolution operator na nabuo ng Hamiltonian, na tinatantiya sa pamamagitan ng simpleng 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
    )
)

Pinalawak na swap test circuit [3]

Inihahanda muna ng circuit ang reference state sa system register habang nananatili ang ancilla sa $|0\rangle$:

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

Pagkatapos, ang paglalagay ng Hadamard gate sa ancilla ay lumilikha ng isang coherent superposition ng dalawang 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*}$$

Sa huli, inilalapat ang controlled time-evolution gate:

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

tanging lamang kapag nasa $|1\rangle$ branch ang ancilla. Samakatuwid,

$$\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*}$$

Sa susunod na code block, ipinapatupad natin ang:

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

na itatalaga bilang $t=d\Delta t$ para sa $d=1,\cdots r-1$ sa execution step.

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))

Observables

Para sa anumang Hermitian system observable $O$, itinatakda natin dito ang mga observable na kakalkulahin:

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

Dahil ito ay dahil ang $O=I$ ay nagbibigay ng overlap element $\mathcal{S}_{0d}=\langle\psi_0|\psi_d\rangle$, habang ang $O=H$ ay nagbibigay ng Hamiltonian element $\mathcal{H}_{0d}=\langle\psi_0|H|\psi_d\rangle$.

Gamit ang $X=|0\rangle\langle 1|+|1\rangle\langle 0|$, mayroon tayong

$$\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*}$$

Katulad nito, gamit ang $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*}$$

Samakatuwid, mayroon tayong

$$\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*}$$

Sa huli, para sa bawat state $|\Psi_{0d}\rangle$, kailangan nating sukatin:

$$\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])

Para sa pagtantiya ng mga matrix element ng $\mathcal{H}$, mas malaki nang husto ang bilang ng Pauli terms kaysa sa mga matrix element ng $\mathcal{S}$.

Maaari na nating bawasan ang bilang ng mga sinusukat na Hamiltonian terms gamit ang isang shifting technique [4]. Hinahati natin ang Hamiltonian bilang:

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

kung saan pinili ang $T$ upang ang reference state ay maging eigenstate nito,

$$\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*}$$

ay ang shifted Hamiltonian matrix element. Samakatuwid, kailangan lang nating sukatin ang $X\otimes (H-T)$ at $Y\otimes (H-T)$. Ang kontribusyon mula sa $T$ ay muling ibinubuo nang classical gamit ang nasukat na nang overlap matrix element na $\mathcal{S}_{0d}$.

Sa halimbawang ito, natural na pagpipilian ang diagonal na bahagi ng Heisenberg Hamiltonian,

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

Dahil ang reference state ay isang computational-basis state, ito ay isang eigenstate ng bawat $Z_iZ_{i+1}$ term.

Gayunpaman, mas mainam na pagpipilian ang isama hindi lamang ang diagonal na mga $ZZ$ terms, kundi pati na rin ang mga $XX+YY$ terms na nag-a-annihilate sa reference state.

Para sa bawat katabing pares, natutugunan ng operator na $XX+YY$ ang

$$\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*}$$

Samakatuwid, nagbibigay ng kontribusyon ang $XX+YY$ term lamang kapag magkaiba ang occupation ng dalawang katabing qubit sa reference bitstring. Kung parehong $0$ o parehong $1$ ang dalawang qubit, na-a-annihilate ng term ang reference state at maaari ring alisin sa pamamagitan ng shifting.

Hayaang $|\psi_{\rm ref}\rangle = |z_0\cdots z_{n-1}\rangle$, kung saan ang $z_i \in \{0,1\}$. Sa gayon, maaari nating piliin ang

$$\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*}$$

Natutugunan pa rin ng operator na ito ang

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

dahil ang mga $ZZ$ terms ay kumikilos nang diagonal sa $|\psi_0\rangle$, habang nagbibigay ng zero ang mga shifted $XX+YY$ terms. Ang katumbas na eigenvalue ay natutukoy lamang, kung gayon, ng mga $ZZ$ terms,

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

Sa pagpiling ito, ang shifted Hamiltonian ay nagiging:

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

Bilang resulta, kailangan lamang sukatin ang mga edge na may magkaibang occupation sa reference state. Lahat ng $ZZ$ terms at lahat ng inactive na $XX+YY$ terms ay muling ibinubuo sa pamamagitan ng overlap contribution na $\tau \mathcal{S}_{0d}$, o nagbibigay ng zero na kontribusyon dahil sa disenyo nito.

Nagbibigay ito ng mas maliit na observable kaysa sa pag-shift lamang ng diagonal na bahagi. Partikular na, para sa isang computational-basis reference state na may localized excitation, ang mga edge lamang na katabi ng excitation ang nananatili sa $H-T$. Kaya naman, ang bilang ng mga Pauli terms sa $X\otimes(H-T)$ at $Y\otimes(H-T)$ ay maaaring bawasan nang malaki, habang ang muling ibinubuong matrix element

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

ay nananatiling eksaktong pareho.

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

Hakbang 2: I-optimize ang problema para sa pagpapatupad sa quantum hardware

Ngayon ay ginagawa nating hardware-oriented template ang abstract na extended-swap-test circuit. Bago iyon, pinapabuti pa natin ang circuit sa abstract level muna.

Ihambing ang pagkakasunod-sunod ng termino ng Hamiltonian

Una, ihahambing natin ang iba't ibang ordering ng mga Pauli terms sa Heisenberg Hamiltonian para sa Hamiltonian simulation. Hindi nagbabago ang mismong Hamiltonian, ngunit apektado ng ordering kung paano nabubuo ang product-formula circuit at kung gaano ka-parallel ang istruktura sa circuit. Halimbawa, inilalayo ng naive na ordering ang lahat ng nearest-neighbor $XX$ terms, pagkatapos lahat ng $YY$ terms, at pagkatapos lahat ng $ZZ$ terms. Inilalagay nito ang mga katabing edge tulad ng $(0,1)$ at $(1,2)$ sa tabi ng isa't isa, kaya't hindi sila maaaring i-execute nang parallel. Binibisita muna ng even-then-odd na ordering ang disjoint na mga even edge, sinusundan ng odd edges, na naglalantad ng parallel two-qubit layers. Lampas pa rito ang even-odd edge-grouped na ordering: para sa bawat edge, pinananatili nitong magkasama ang local na $XX$, $YY$, at $ZZ$ terms, habang binibisita pa rin ang even edges bago ang odd edges. Inaasahan nating babawasan ng even-then-odd at even-odd edge-grouped orderings ang circuit depth sa pamamagitan ng paglalantad ng parallel two-qubit layers, at babawasan pa ng edge-grouped ordering ang Trotter error dahil ang local two-qubit interaction sa parehong edge ay hinahawakan bilang compact block.

Apektado rin ng Hamiltonian ordering ang Trotter error. Kung inilalagay natin ang mga noncommuting terms sa tabi ng isa't isa, mas madalas magaganap ang basis transitions, na nagdudulot ng mas maraming Trotter error. Sa pamamagitan ng pagpapangkat ng mga terms na nangangailangan ng parehong Pauli-basis transformation, maiiwasan ang mga redundant na basis-changes.

Dito, ginagamit ng paghahambing ang pinakamalaking evolution time na lumilitaw sa first-row Krylov estimates, $t_{\max} = (r-1)\Delta t,$ gamit ang parehong transpilation condition.

Upang sukatin ang Trotter error, ginagamit natin ang process infidelity sa pagitan ng Trotterized circuit at ng eksaktong Hamiltonian evolution,

$$\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.

Gumagamit ang diagnostic na ito ng dense matrices, kaya't angkop ito para sa maliit na 12-qubit na halimbawang ito ngunit hindi ito nilalayong maging scalable na subroutine.

# 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

Pinananatili muna nating fixed ang synthesis rule sa isang first-order Trotter step at pinapabago lamang ang Pauli-term ordering. Ang layunin ng paghahambing na ito ay pangunahing makita kung gaano kalaki ang mababawasan sa circuit depth at two-qubit cost sa pamamagitan ng paglalantad ng disjoint na nearest-neighbor edges sa transpiler.

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

Sa reps=1, napapansin natin na parehong pinapababa ng even-then-odd at even-odd edge-grouped orderings ang depth mula 15 patungong 6 sa pamamagitan ng paglalantad ng parallel two-qubit layers, habang nananatiling pareho ang two-qubit gate count sa lahat ng tatlong ordering. Gayunpaman, malapit sa 1 ang infidelity ng lahat ng tatlong ordering, kaya't sinusuri natin ang bilang ng Trotter repetitions upang mas malinaw na ihiwalay ang mga ordering. Habang lumalaki ang bilang ng repetitions, mas mabilis bumaba ang infidelity ng even-odd edge-grouped na ordering kaysa sa dalawang iba pa, umaabot sa 0.045 sa reps=8 kumpara sa 0.075 para sa naive at even-then-odd na ordering.

Ihambing ang product-formula synthesis

Susunod, sinusuri natin ang iba't ibang advanced na setting ng Trotterization, na may fixed na Hamiltonian ordering sa even-odd edge-grouped ordering. Isinasaalang-alang natin ang first-order Lie-Trotter, second-order Suzuki-Trotter, at fourth-order 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

Nagbibigay ang Lie-Trotter ng pinakamababaw na circuit ngunit may pinakamalaking error, habang mas tumpak ang fourth-order Suzuki-Trotter ngunit nagpapataas ito ng circuit depth. Para sa natitirang bahagi ng tutorial, pinipili natin ang second-order Suzuki-Trotter dahil nagbibigay ito ng maliit na depth na circuit habang malaking-lakas na binabawasan ang Trotter error kumpara sa first-order formula.

Alisin ang controlled time-evolution gate

Sa extended swap test, ang pagkontrol sa time evolution gamit ang iisang ancilla qubit ay nangangailangan ng ancilla na kumontrol sa maraming gates sa buong system. Maaari itong magdulot ng malaking routing overhead at, sa pinakamasamang kaso, epektibong nangangailangan ng all-to-one connectivity. Upang maiwasan ito, posible ang karagdagang optimization sa pamamagitan ng pagpapalit ng controlled time-evolution gate ng control-free na bersyon, na sinasamantala ang symmetry ng Hamiltonian. Masdan natin ang sumusunod na circuit.

Dito, inihahanda ng $B_{\rm ref}$ ang reference state, $B_{\rm ref}|0^n\rangle = |\psi_{0}\rangle$.

Sa halip na unang ihanda ang $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)|\psi_0\rangle$ at pagkatapos ilapat ang $U(t)$ sa $|1\rangle$ branch lamang, direktang inihahanda ng circuit ang dalawang branch bilang

$$\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*}$$

Unang inilalapat ng circuit ang Hadamard gate sa ancilla at inihahanda ang reference state sa $|1\rangle$ branch lamang:

$$\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*}$$

Pagkatapos, inilalapat ang uncontrolled time-evolution operator sa parehong 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*}$$

Dahil pinapanatili ng Hamiltonian ang excitation number, makikita natin na ang $|0^n\rangle$ ay eigenstate nito, at kaya lamang naipupundo ng evolution operator ay isang phase sa ilalim ng Hamiltonian:

$$\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*}$$

Susunod, inilalapat ang $B_{\rm ref}$ sa $|0\rangle$ branch lamang.

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

Sa puntong ito, may karagdagang relative phase ang dalawang branch. Upang alisin ito, inilalapat natin ang ancilla phase gate

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

Binabago nito ang state bilang

$$\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*}$$

Kaya naman, maliban sa isang hindi mahalagang global phase, naihanda na natin sa huli ang

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

Sa susunod na code block, ipinapatupad natin ang control-free circuit na ito.

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

Transpilation

Ngayon, ini-transpile natin ang controlled- at control-free na circuits upang maging executable sa hardware. Ihambing natin ang resulta ng mga na-transpile na circuits.

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>
})

Hakbang 3: Isagawa gamit ang Qiskit primitives

Ang susunod na hakbang ay ang pagsumite ng parehong parameterized circuit para sa ilang halaga ng $d$. Para sa bawat $d$, tinatantiya natin ang apat na expectation value: $X\otimes I$, $Y\otimes I$, $X\otimes (H-T)$, at $Y\otimes (H-T)$. Ang apat na numerong ito ay pagkatapos pinagsasama upang bumuo ng complex first-row elements na $\mathcal{S}_{0d}$ at $\mathcal{H}_{0d}$.

Dito, $\mathcal{S}_{00}=1$ at $\mathcal{H}_{00}=\langle\psi_{0}|(H-T)|\psi_{0}\rangle=\langle\psi_{0}|H|\psi_{0}\rangle-\tau$ ay maaaring kalkulahin nang classical dahil sparse ang $\lvert\psi_0\rangle$, kaya't nilalaktawan natin ang kaso ng $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']

Hakbang 4: I-post-process at ibalik ang resulta sa nais na classical format

Pagkatapos tantiyahin ang mga projected matrices, nire-regularize natin at nilulutas ang GEVP

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

Ang pinakamaliit na generalized eigenvalue ay nagbibigay ng KQD estimate ng ground-state energy.

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

Ngayon nilulutas natin ang generalized eigenvalue problem gamit ang mga matrices na muling ibinuo mula sa mga estimate ng circuit. Sa isang ideal na statevector calculation na may eksaktong real-time evolution, dapat itong maireproduce ang eksaktong projected result. Sa praktika, maaaring magmula ang mga deviation sa Trotterization, sampling error, at numerical instability ng overlap matrix.

Sinusuri natin kung paano nag-converge ang energy habang pinapataas natin ang dimension ng Krylov subspace.

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)

Large-scale na halimbawa sa hardware

Gumamit ang nakaraang seksyon ng 12-qubit na modelo upang magamit ang statevector simulation bilang diagnostic. Ngayon ay pinapalaki natin ang parehong KQD workflow sa isang 30-qubit na Heisenberg chain at inihahanda ang workload para sa execution sa IBM Quantum hardware.

Mga Hakbang 1-4 na pinagsama sa iisang code block

Dito, pinagsasama-sama na natin ang lahat ng detalyeng ito sa iisang workflow sa mas malaking sukat, na pagkatapos ay pinapatakbo sa ating tunay na quantum hardware. Sa seksyong ito, inilalapat natin ang makatotohanang mga setting ng error mitigation upang pahusayin ang reliability ng mga resulta. Dahil ang mga matrix element na katumbas ng iba't ibang halaga ng $d$ ay maaaring suriin nang parallel, ginagamit natin ang Batch mode upang epektibong i-execute ang mga ito.

# -------------------------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

Apendise: Ang pananaw ng Hamiltonian-function (spectral-filter)

Ipinakita ng pangunahing workflow ang KQD nang operational: bumuo ng Krylov basis mula sa real-time-evolved states, tantiyahin ang projected matrices na $\mathcal{H}$ at $\mathcal{S}$, at lutasin ang GEVP. Binibisita muli ng appendix na ito ang parehong computation mula sa isang complementary na anggulo na nagpapaliwanag kung bakit gumagana ang KQD: ang pananaw na Hamiltonian-function, o spectral-filter [3], [5]. Muling ginagamit nito ang 12-qubit na modelo, time step $\Delta t$, at ang solusyon sa Krylov na nakuha na sa itaas; walang kinakailangang bagong execution ng circuit.

Ang reference state bilang isang energy distribution

Hayaang ang Hamiltonian ay may eigendecomposition na

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

na may energy eigenstates na $|E_m\rangle$. Anumang reference state ay maaaring i-expand sa eigenbasis na ito,

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

kaya't nagdadala ito ng spectral weight na $p_m=|a_m|^2$ sa bawat energy na $E_m$. Ang reference energy ay ang mean ng distribution na ito, $\langle H\rangle_{0}=\sum_m p_m E_m$.

Ang eigendecomposition ng isang pangkalahatang $n$-qubit Hamiltonian ay exponentially mahal, kaya't ang larawang ito ay diagnostic lamang, na hindi bahagi ng algorithm. Dito, gayunpaman, magagawa nating i-compute ito nang mura para sa parehong $n=12$ problema: pinapanatili ng Heisenberg Hamiltonian ang kabuuang bilang ng excitation, at ang reference na $|000001000000\rangle$ ay nagdadala ng iisang excitation, kaya't buong spectral content nito ay nasa single-excitation subspace na ang dimension ay lumalaki lamang nang linear kasama ng $n$. Muli nating ginagamit ang eksaktong single-excitation block (nagamit na bilang benchmark sa itaas) at kinukuha ang reference distribution na $\{(E_m, p_m)\}$ sa loob ng subspace na iyon.

# 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

Natututo ang KQD ng isang filter na humuhubog muli sa distribusyong ito

Ang Hamiltonian function na $f(H)$ ay tinutukoy sa pamamagitan ng spectral calculus,

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

o sa madaling salita, isang weighted sum ng eigenprojectors. Ang paglalapat nito sa reference ay nagbabago ng anyo ng bawat spectral amplitude, $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*}$$

Kung matalim na peaked ang $f$ sa pinakamababang energy ($f(E_0)=1$ at $f(E_m)\approx 0$ sa iba pang kaso) kung gayon ang $f(H)$ ay kikilos bilang ground-state projector, at ang normalized na output ay magiging (malapit sa) ground state. Samakatuwid, ang isang magandang low-pass spectral filter sa energy ang eksaktong kailangan natin.

Hindi tinutukoy ng KQD ang $f$ nang maaga. Sa halip, ini-expand nito ang filter sa real-time-evolution basis,

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

isang trigonometric na function ng enerhiya kung saan ang mga coefficient nitong $\{c_\ell\}$ ay eksaktong ang GEVP eigenvector na na-solve para sa itaas. Ang pag-minimize ng Rayleigh quotient $\mathbf{c}^\dagger\mathcal{H}\mathbf{c}/\mathbf{c}^\dagger\mathcal{S}\mathbf{c}$ ay katumbas din ng pag-aaral ng filter na pinakamahusay na sumusupress sa excited-state weight ng reference. Mas malaking Krylov dimension $r$ ay nagbibigay sa filter ng mas maraming degrees of freedom at mas matalim na peak sa ground-state energy.

Sinusuri ng helper sa ibaba ang natutunang filter na ito sa isang energy axis; ilalapat natin ito pagkatapos sa reference distribution na nakuha sa itaas.

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

I-visualize ang filter at ang flexibility nito

Ipapakita muna natin ang natutunang filter sa buong Krylov dimension na ginamit sa itaas, pagkatapos ay susubaybayan kung paano ito humahasa habang lumalaki ang dimension $r$.

Ipinapakita ng mga bar ang reference spectral weights $p_m$ (bago) at ang filtered weights $p_m|f_{\rm KQD}(E_m)|^2$ (pagkatapos), kasama ang natutunang filter intensity $|f_{\rm KQD}(E)|^2$ sa isang tuloy-tuloy na energy axis. Inikonsentra ng filter ang weight sa pinakamababang enerhiya ng single-excitation subspace — ang parehong enerhiya na kinombergehan ng KQD estimate sa small-scale example. Tandaan na ito ay ang ground state sa loob ng single-excitation sector, na siyang kaugnay na target para sa excitation-conserving reference state na ito, hindi ang global ground state.

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()

Pagtaas ng Krylov dimension: flexibility ng natutunang function

Tandaan na ang natutunang filter ay isang trigonometric polynomial sa enerhiya na may $r$ coefficients,

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

Ang Krylov dimension $r$ ay eksaktong ang bilang ng free coefficients, kaya kontrolado nito ang flexibility ng function. Ang maliit na $r$ ay makakagawa lamang ng malawak, mahinahong nagbabagong filter na tumutulo ang weight papunta sa mababang excited states; kapag tumaas ang $r$, maaaring bumuo ang filter ng mas makitid na peak sa target energy at supresahin nang mas agresibo ang natitirang excited-state weight. Ito ang spectral-filter counterpart ng energy convergence na naobserbahan sa small-scale example: habang lumalaki ang $r$, ang filtered distribution ay kumukolapsa papunta sa subspace ground state at ang estimated energy ay bumababa patungo dito.

Muli nating gagamitin ang mga projected matrix na na-estimate na sa itaas at simpleng lulutasin ang GEVP sa bawat leading $r\times r$ block, pagkatapos ay susuriin at ipa-plot ang katumbas na filter.

# 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

Mga susunod na hakbang

Kung nakita mong interesante ang gawaing ito, maaaring interesado ka sa mga sumusunod na materyal:

Mga Rekomendasyon

• Sample-based Krylov quantum diagonalization of a fermionic lattice model
• Ground-state energy estimation of the Heisenberg chain with VQE
• Quantum diagonalization algorithms course

References

[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).