格子ハミルトニアンのクリロフ量子対角化

使用時間の目安: HeronまたはNighthawkプロセッサで70分（注: これは目安です。実行時間は異なる場合があります。）

学習成果

• クリロフ量子対角化（KQD）を、スペクトル・フィルターとして機能する有限のハミルトニアン関数を学習することとして解釈する方法。

• 拡張スワップ・テスト測定を用いて射影ハミルトニアン行列と重なり行列を構築する方法。

• 得られた一般化固有値問題（GEVP）を解き、格子ハミルトニアンの基底状態エネルギー推定値を求める方法。

前提条件

• Qiskitプリミティブ

• VQEを用いたハイゼンベルク鎖の基底状態エネルギー推定

• 量子対角化アルゴリズムコース

背景

このチュートリアルでは、Qiskitパターンの文脈でクリロフ量子対角化アルゴリズム（KQD）を実装する方法を示します。まずアルゴリズムの背後にある理論を学び、その後QPU上での実行のデモンストレーションを行います。

多体ハミルトニアンの低エネルギー特性を推定することは、量子シミュレーションにおける中心的な課題です。 例えば、基底状態エネルギーや低励起状態は、化学的安定性、磁気秩序、量子相転移、および物質の応答特性と直接関係しています。 古典コンピュータでは、ヒルベルト空間の次元が軌道数やスピン数に対して指数関数的に増大するため、直接対角化はすぐに現実的ではなくなります。

この問題に対しては、いくつかの量子コンピューティングのアプローチがあります。 変分量子固有値ソルバー（VQE）などの近い将来向けの変分法は、比較的浅いパラメータ化回路を使用しますが、多くの量子回路評価を伴う非線形な古典最適化ループが必要です。 一方、量子位相推定（QPE）は、厳密な保証を伴うより直接的な固有値推定の手段を提供しますが、標準的なQPEは長いコヒーレント回路を必要とし、主に誤り耐性量子コンピュータに適しています。 KQDはこれら2つのアプローチの中間に位置します。位相推定ベースのアルゴリズムと同様にリアルタイムのハミルトニアン発展を使用しますが、完全な位相推定の代わりに、古典的に解ける小規模な射影固有値問題を用います。

$n$量子ビットのハミルトニアン $H$ と参照状態 $\lvert \psi_{0}\rangle$ を考えます。 KQD法は、リアルタイムに発展させた状態からクリロフ部分空間を構築します。

$$\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$ はクリロフ次元、$\Delta t$ は時間刻みです。 クリロフ部分空間内の任意の状態は、これらの基底状態の線形結合として表されます。

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

射影された重なり行列とハミルトニアン行列を定義します。 これらの要素は、量子回路測定を用いて推定されます。

$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はリアルタイム発展に依存するため、通常より深い回路を必要とします。 その代わり、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. 学習されたクリロフ・フィルターを評価する。

4. フィルター値をスペクトル重みに変換する。

5. 参照エネルギー分布とフィルター後のエネルギー分布をプロットする。

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

このチュートリアルの前半では、ローカルの状態ベクトル・シミュレーターを使用してKQD法を実演します。後半では、実際の量子バックエンドを使用して、ユーティリティ・スケールの問題に取り組みます。

また、バックエンド固有のトランスパイルを実演し、得られた回路を検証するために、フェイク・バックエンドも定義します。

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

backend = FakeMarrakesh()

小規模シミュレーター例

ステップ1: 古典的な入力を量子問題にマッピングする

ハミルトニアンと参照状態

この例では、12量子ビットの開放境界ハイゼンベルク鎖（$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*}$$

参照状態とします。 上で定義したハイゼンベルク・ハミルトニアンは励起数の総和を保存するため、参照状態は単一励起部分空間内にとどまり、その次元は量子ビット数に対して線形にしか増加しません。したがって、その部分空間に制限したハミルトニアンを対角化することで、厳密な基底状態エネルギーを効率的に計算でき、これはKQD推定値を検証するための参照値としてのみ使用します。KQDワークフロー自体は、Qiskitプリミティブを用いて射影された行列要素を推定し、得られた射影問題を古典的に解きます。

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

クリロフ次元を $r=10$、1時間刻みあたりのトロッター・ステップ数を $5$ に設定します。これは、最も深い回路（$t_{\max}=(r-1)\Delta t$）を許容範囲に保ちながら低励起スペクトルを解像するのに十分な大きさのクリロフ空間であり、また、その最も深い回路において離散化誤差を小さく保つのに十分なトロッター・ステップ数です。

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

回路の構築

ここでは、行列要素 $\mathcal{H}_{k\ell}$ と $\mathcal{S}_{k\ell}$ を推定するための回路を構築します。 $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*}$$

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

ここでは、拡張スワップ・テストと呼ばれる回路を紹介します。この回路は、

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

時間発展

ハミルトニアンによって生成される時間発展演算子を、単純なリー・トロッター化によって近似し、実現します。

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

次に、アンシラにアダマール・ゲートを適用すると、2つの分岐のコヒーレントな重ね合わせが生成されます。

$$\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$ 分岐にある場合にのみ適用されます。したがって、

$$\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{ は } \operatorname{Re}\mathcal{S}_{0d},\\ Y \otimes I \text{ は } \operatorname{Im}\mathcal{S}_{0d},\\ X \otimes H \text{ は } \operatorname{Re}\mathcal{H}_{0d},\\ Y \otimes H \text{ は } \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}$ 行列要素の推定では、パウリ項の数が $\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$ 項は、参照ビット列において隣接する2つの量子ビットの占有数が異なる場合にのみ寄与します。 2つの量子ビットがともに $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$ 項はゼロを与えるからです。 したがって、対応する固有値は $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}$ を通じて再構成されるか、構成上ゼロの寄与を与えます。

これにより、対角部分のみをシフトする場合よりも小さなオブザーバブルが得られます。 特に、局在した励起を持つ計算基底の参照状態では、その励起に隣接するエッジのみが $H-T$ に残ります。 したがって、$X\otimes(H-T)$ と $Y\otimes(H-T)$ に含まれるパウリ項の数を大幅に削減できる一方で、再構成された行列要素

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

ステップ 2: 量子ハードウェア実行のための問題の最適化

ここでは、抽象的な拡張スワップ・テスト回路を、ハードウェア指向のテンプレートに変換します。 その前に、まず抽象レベルで回路をさらに最適化します。

ハミルトニアン項の順序の比較

まず、ハミルトニアン・シミュレーションのために、ハイゼンベルク・ハミルトニアン内のパウリ項の異なる順序を比較します。 ハミルトニアン自体は変わりませんが、順序は積公式回路がどのように生成されるか、また回路内の構造をどの程度並列化できるかに影響します。 例えば、素朴な順序では、まずすべての最近接 $XX$ 項をリストし、次にすべての $YY$ 項、最後にすべての $ZZ$ 項を並べます。 これにより、$(0,1)$ や $(1,2)$ のような隣接するエッジが互いに隣り合って配置されるため、並列実行できません。 偶数・奇数順序（even-then-odd）は、まず互いに素な偶数エッジを訪れ、続いて奇数エッジを訪れるため、並列な2量子ビット層が現れます。 偶数・奇数エッジ・グループ化順序（even-odd edge-grouped）はさらに一歩進み、各エッジについて局所的な $XX$、$YY$、$ZZ$ 項をまとめて扱いながら、依然として奇数エッジより先に偶数エッジを訪れます。 偶数・奇数順序と偶数・奇数エッジ・グループ化順序は、並列な2量子ビット層を露出させることで回路深さを削減すると予想され、エッジ・グループ化順序はさらに、同一エッジ上の局所的な2量子ビット相互作用がコンパクトなブロックとして処理されるため、トロッター誤差も削減すると予想されます。

ハミルトニアンの順序は、トロッター誤差にも影響します。 非可換な項を隣り合わせに配置すると、基底変換がより頻繁に発生し、より大きなトロッター誤差を引き起こします。 同じパウリ基底変換を必要とする項をグループ化することで、冗長な基底変更を回避できます。

ここでは、同じトランスパイル条件を用いて、第1行のクリロフ推定に現れる最大の発展時間 $t_{\max} = (r-1)\Delta t,$ を用いて比較を行います。

トロッター誤差を測定するために、トロッター化された回路と厳密なハミルトニアン発展との間のプロセス非忠実度を使用します。

$$\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次のトロッター・ステップ1回に固定し、パウリ項の順序のみを変化させます。 この比較の主な目的は、互いに素な最近接エッジをトランスパイラーに露出させることで、回路深さと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量子ビット・ゲート数は3つの順序すべてで変わらないことが観察されます。 しかし、3つの順序すべてで非忠実度が1に近いため、トロッター繰り返し回数を掃引して順序をより明確に区別します。 繰り返し回数が増えるにつれて、even-odd edge-grouped 順序の非忠実度は他の2つよりも速く低下し、reps=8 で0.045に達するのに対し、素朴な順序と偶数・奇数順序では0.075となります。

積公式合成の比較

次に、ハミルトニアンの順序を偶数・奇数エッジ・グループ化順序に固定した上で、トロッター化のさまざまな高度な設定を検討します。 1次のリー・トロッター、2次の鈴木・トロッター、4次の鈴木・トロッターを考えます。

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

リー・トロッターは最も浅い回路を与えますが、誤差が最も大きく、一方で4次の鈴木・トロッターはより正確ですが回路深さが増加します。 このチュートリアルの残りの部分では、2次の鈴木・トロッターを選択します。これは、1次の公式と比較してトロッター誤差を大幅に削減しつつ、小さな深さの回路を与えるためです。

制御時間発展ゲートの除去

拡張スワップ・テストにおいて、単一のアンシラ量子ビットで時間発展を制御するには、アンシラがシステム全体にわたる多くのゲートを制御する必要があります。 これにより、大きなルーティング・オーバーヘッドが生じる可能性があり、最悪の場合、実質的にオール・ツー・ワンの接続性が必要になります。 これを回避するため、ハミルトニアンの対称性を利用して、制御時間発展ゲートを制御なしバージョンに置き換えることで、さらなる最適化が可能です。 次の回路を見てみましょう。

ここで、$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$ 分岐にのみ $U(t)$ を適用する代わりに、この回路は2つの分岐を直接次のように準備します。

$$\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$ 分岐にのみ参照状態を準備します。

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

次に、制御なしの時間発展演算子が両方の分岐に適用されます。

$$\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$ 分岐にのみ適用されます。

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

この時点で、2つの分岐には余分な相対位相が生じています。それを取り除くために、アンシラ位相ゲートを適用します。

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

ステップ 3: Qiskit プリミティブを使用した実行

次のステップでは、同じパラメータ化回路を複数の $d$ の値に対して投入します。それぞれの $d$ に対して、$X\otimes I$、$Y\otimes I$、$X\otimes (H-T)$、$Y\otimes (H-T)$ の4つの期待値を推定します。これら4つの数値は、複素数の第1行要素 $\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

ここでは、回路推定値から再構成された行列を用いて、一般化固有値問題を解きます。 厳密なリアルタイム発展を用いた理想的な状態ベクトル計算では、これは厳密な射影結果を再現するはずです。 実際には、トロッター化、サンプリング誤差、重なり行列の数値的不安定性から偏差が生じる可能性があります。

クリロフ部分空間の次元を増やすにつれて、エネルギーがどのように収束するかを観察します。

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を1つのコード・ブロックにまとめる

ここでは、これらの詳細をすべて1つのワークフローにまとめ、より大きなスケールで実際の量子ハードウェア上で実行します。このセクションでは、結果の信頼性を向上させるために、現実的なエラー緩和設定を適用します。異なる $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を操作的に提示しました。すなわち、リアルタイムに発展させた状態からクリロフ基底を構築し、射影行列 $\mathcal{H}$ と $\mathcal{S}$ を推定し、GEVPを解くという流れです。 この付録では、KQDがなぜうまくいくのかを説明する補完的な視点、すなわちハミルトニアン関数、あるいはスペクトル・フィルターの観点 [3]、[5] から、同じ計算を再考します。 ここでは、上ですでに得た12量子ビットのモデル、時間刻み $\Delta t$、およびクリロフ解を再利用するため、新たな回路実行は必要ありません。

エネルギー分布としての参照状態

ハミルトニアンが次の固有分解を持つとします。

$$\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 の固有ベクトルそのものです。 したがって、レイリー商 $\mathbf{c}^\dagger\mathcal{H}\mathbf{c}/\mathbf{c}^\dagger\mathcal{S}\mathbf{c}$ を最小化することは、参照状態の励起状態成分を最も良く抑制するフィルターを学習することと同じです。 クリロフ次元 $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

フィルターとその柔軟性を可視化する

まず、上で使用した完全なクリロフ次元における学習済みフィルターを示し、その後、次元 $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()

クリロフ次元を増やす：学習された関数の柔軟性

学習されたフィルターは、$r$ 個の係数を持つエネルギーの三角多項式であったことを思い出してください。

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

クリロフ次元 $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

次のステップ

この内容に興味を持った方は、以下の資料にも興味があるかもしれません:

おすすめ

• フェルミオン格子モデルのサンプルベース・クリロフ量子対角化
• 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).