{ "cells": [ { "cell_type": "markdown", "id": "afaf178f", "metadata": {}, "source": [ "---\n", "title: Krylov quantum diagonalization of lattice Hamiltonians\n", "description: Learn Krylov quantum diagonalization for lattice Hamiltonians with the concept of Hamiltonian function.\n", "---\n", "\n", "{/* cspell:ignore Fomichev GEVP Izmaylov KQDOLH Kang Nakatsukasa Toeplitz eigenbasis eigendecomposition eigvals eigvecs evals noncommuting orthogonalizer sharex thresholded */}\n", "\n", "# Krylov quantum diagonalization of lattice Hamiltonians\n", "\n", "*Usage estimate: 70 minutes on a Heron or Nighthawk processor (NOTE: This is an estimate only. Your runtime might vary.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "aa4ba954", "metadata": {}, "source": [ "## Learning outcomes\n", "\n", "* How to interpret Krylov quantum diagonalization (KQD) as learning a finite Hamiltonian function that acts as a spectral filter.\n", "* How to build the projected Hamiltonian and overlap matrices with extended swap-test measurements.\n", "* How to solve the resulting generalized eigenvalue problem (GEVP) and recover a ground-state energy estimate for a lattice Hamiltonian.\n", "\n", "## Prerequisites\n", "\n", "* [Qiskit primitives](/docs/guides/primitives)\n", "* [Ground-state energy estimation of the Heisenberg chain with VQE](/docs/tutorials/spin-chain-vqe)\n", "* [Quantum diagonalization algorithms course](/learning/courses/quantum-diagonalization-algorithms/krylov)\n", "\n", "## Background\n", "\n", "This tutorial demonstrates how to implement the Krylov quantum diagonalization (KQD) algorithm within the context of Qiskit patterns.\n", "You will first learn the theory behind the algorithm and then see a demonstration of its execution on a QPU.\n", "\n", "Estimating the low-energy properties of many-body Hamiltonians is a central task in quantum simulation.\n", "For example, ground-state energies and low-lying excitations are directly related to chemical stability, magnetic ordering, quantum phase transitions, and material response.\n", "On a classical computer, the Hilbert-space dimension grows exponentially with the number of orbitals or spins, so direct diagonalization quickly becomes impractical.\n", "\n", "There are several quantum-computing approaches to this problem.\n", "Near-term variational methods, such as the [variational quantum eigensolver (VQE)](/docs/tutorials/spin-chain-vqe), use relatively shallow parameterized circuits, but they require a nonlinear classical optimization loop with many quantum-circuit evaluations.\n", "At the other end, [quantum phase estimation (QPE)](https://arxiv.org/abs/quant-ph/0604193) provides a more direct route to eigenvalue estimation with rigorous guarantees, but standard QPE requires long coherent circuits and is mainly suited to fault-tolerant quantum computers.\n", "KQD sits between these two approaches: it uses real-time Hamiltonian evolution, as in phase-estimation-based algorithms, but replaces full phase estimation with a compact projected eigenvalue problem that can be solved classically.\n", "\n", "Consider an $n$-qubit Hamiltonian $H$ and a reference state $\\lvert \\psi_{0}\\rangle$.\n", "The KQD method builds a Krylov subspace from real-time evolved states,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\lvert \\psi_\\ell\\rangle = e^{-i\\ell\\Delta t H}\\lvert \\psi_{0}\\rangle,\n", "\\qquad\n", "\\ell = 0,1,\\ldots,r-1,\n", "\\end{equation*}\n", "$$\n", "\n", "where $r$ is the Krylov dimension and $\\Delta t$ is the time step.\n", "Any state in the Krylov subspace is then represented as a linear combination of these basis states,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\lvert \\psi(\\mathbf{c})\\rangle\n", "=\n", "\\frac{\\sum_{\\ell=0}^{r-1} c_\\ell \\lvert \\psi_\\ell\\rangle}\n", "{\\left\\|\\sum_{\\ell=0}^{r-1} c_\\ell \\lvert \\psi_\\ell\\rangle\\right\\|},\n", "\\end{equation*}\n", "$$\n", "\n", "where the denominator normalizes the state.\n", "\n", "With simple algebra, we can see that the corresponding energy is written as the Rayleigh quotient,\n", "\n", "$$\n", "\\begin{equation*}\n", "E(\\mathbf{c})\n", "=\\langle \\psi(\\mathbf{c})|H|\\psi(\\mathbf{c}) \\rangle=\n", "\\frac{\n", "\\sum_{k,\\ell} c_k^* c_\\ell \\langle \\psi_k\\vert H\\vert\\psi_\\ell\\rangle\n", "}{\n", "\\sum_{k,\\ell} c_k^* c_\\ell \\langle \\psi_k\\vert\\psi_\\ell\\rangle\n", "}\n", "=\n", "\\frac{\n", "\\mathbf{c}^{\\dagger}\\mathcal{H}\\mathbf{c}\n", "}{\n", "\\mathbf{c}^{\\dagger}\\mathcal{S}\\mathbf{c}\n", "}.\n", "\\end{equation*}\n", "$$\n", "\n", "Here, the matrices $\\mathcal{S}$ and $\\mathcal{H}$,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{S}_{k\\ell}=\\langle \\psi_k\\vert\\psi_\\ell\\rangle,\n", "\\qquad\n", "\\mathcal{H}_{k\\ell}=\\langle \\psi_k\\vert H\\vert\\psi_\\ell\\rangle\n", "\\end{equation*}\n", "$$\n", "\n", "define the projected overlap and Hamiltonian matrices.\n", "Their entries are estimated using quantum-circuit measurements.\n", "\n", "We aim to find the coefficient $\\mathbf{c}$ that gives the minimum $E(\\mathbf{c})$:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\min_{\\mathbf{c}\\neq \\mathbf{0}} E(\\mathbf{c}).\n", "\\end{equation*}\n", "$$\n", "\n", "By the Rayleigh-Ritz theorem, this minimization is equivalent to solving the generalized eigenvalue problem (GEVP),\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{H}\\mathbf{c}=E \\mathcal{S}\\mathbf{c}.\n", "\\end{equation*}\n", "$$\n", "\n", "Note that the dimension $r$ can be small enough for a classical computer to solve the GEVP.\n", "\n", "This is the same variational principle used in classical subspace diagonalization, but here the basis states are generated by quantum time evolution.\n", "Compared with VQE, KQD usually requires deeper circuits because it relies on real-time evolution.\n", "In return, KQD avoids nonlinear parameter optimization and iterative quantum-hardware execution, and it systematically improves as the projected subspace is enlarged.\n", "The algorithm has been demonstrated at large scale on existing quantum hardware [\\[2\\]](#references), and its performance can be analyzed with provable guarantees [\\[1\\]](#references).\n", "\n" ] }, { "cell_type": "markdown", "id": "3285eedc", "metadata": {}, "source": [ "## Requirements\n", "\n", "Before starting this tutorial, ensure that you have the following installed:\n", "\n", "* Qiskit SDK v2.3 or later with [visualization](/docs/api/qiskit/visualization) support\n", "* Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n", "* SciPy (`pip install scipy`)\n", "* Matplotlib (`pip install matplotlib`)\n", "* Pandas (`pip install pandas`)\n", "\n", "Hardware execution requires `qiskit-ibm-runtime` and access to an IBM Quantum® account.\n", "\n" ] }, { "cell_type": "markdown", "id": "047d8dcd", "metadata": {}, "source": [ "## Setup\n", "\n", "The setup cell imports the required modules and defines helper functions for the workflow:\n", "\n", "1. build the Heisenberg Hamiltonian;\n", "2. solve the thresholded GEVP;\n", "3. evaluate the learned Krylov filter;\n", "4. convert filter values into spectral weights;\n", "5. plot the reference and filtered energy distributions.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "7bb044cf", "metadata": {}, "outputs": [], "source": [ "from __future__ import annotations\n", "\n", "import warnings\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import scipy.linalg as la\n", "import matplotlib.pyplot as plt\n", "\n", "from qiskit import QuantumCircuit, transpile\n", "from qiskit.circuit import Parameter\n", "from qiskit.circuit.library import PauliEvolutionGate\n", "from qiskit.primitives import StatevectorEstimator\n", "from qiskit.quantum_info import Operator, SparsePauliOp\n", "from qiskit.synthesis import LieTrotter, SuzukiTrotter\n", "from qiskit.transpiler import PassManager, Layout\n", "from qiskit.transpiler.passes import CommutativeOptimization\n", "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "from qiskit_ibm_runtime import QiskitRuntimeService, EstimatorV2, Batch\n", "from qiskit_ibm_runtime.fake_provider import FakeMarrakesh\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "\n", "def make_heisenberg_hamiltonian(\n", " num_qubits: int,\n", " coupling: float = 1.0,\n", ") -> SparsePauliOp:\n", " \"\"\"Make a Heisenberg Hamiltonian for a 1D chain of qubits with nearest-neighbor interactions.\"\"\"\n", " terms: list[tuple[str, complex]] = []\n", "\n", " def append_term(q0: int, q1: int, pauli: str):\n", " label = [\"I\"] * num_qubits\n", " label[num_qubits - 1 - q0] = pauli[0]\n", " label[num_qubits - 1 - q1] = pauli[1]\n", " terms.append((\"\".join(label), coupling))\n", "\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " for q in range(num_qubits - 1):\n", " append_term(q, q + 1, pauli)\n", " return SparsePauliOp.from_list(terms).simplify()\n", "\n", "\n", "def _basis_state_transition_amplitude_sparse(\n", " hamiltonian: SparsePauliOp,\n", " bra_state: int,\n", " ket_state: int,\n", ") -> complex:\n", " \"\"\"Evaluate for computational-basis states.\"\"\"\n", " num_qubits = hamiltonian.num_qubits\n", " amplitude = 0.0 + 0.0j\n", " for pauli, coeff in zip(hamiltonian.paulis, hamiltonian.coeffs):\n", " new_state = ket_state\n", " phase = 1.0 + 0.0j\n", " for q in range(num_qubits):\n", " x = bool(pauli.x[q])\n", " z = bool(pauli.z[q])\n", " if not x and not z:\n", " continue\n", " bit = (new_state >> q) & 1\n", " if x and z:\n", " # Y|0> = i|1>, Y|1> = -i|0>\n", " phase *= 1j if bit == 0 else -1j\n", " new_state ^= 1 << q\n", " elif x:\n", " new_state ^= 1 << q\n", " else:\n", " # Z|0> = |0>, Z|1> = -|1>\n", " if bit:\n", " phase *= -1\n", " if new_state == bra_state:\n", " amplitude += coeff * phase\n", " return amplitude\n", "\n", "\n", "def basis_state_expectation_sparse(\n", " hamiltonian: SparsePauliOp,\n", " bitstring: str,\n", ") -> complex:\n", " \"\"\"Evaluate .\"\"\"\n", " state = int(bitstring, 2)\n", " return _basis_state_transition_amplitude_sparse(\n", " hamiltonian,\n", " bra_state=state,\n", " ket_state=state,\n", " )\n", "\n", "\n", "def diagonalize_single_1_subspace(\n", " hamiltonian: SparsePauliOp,\n", ") -> np.ndarray:\n", " \"\"\"Diagonalize the Hamiltonian projected onto the single-excitation subspace.\"\"\"\n", " num_qubits = hamiltonian.num_qubits\n", "\n", " # Integer basis states |...010...>, with the excitation at qubit k.\n", " basis = [1 << k for k in range(num_qubits)]\n", "\n", " h_single = np.empty((num_qubits, num_qubits), dtype=complex)\n", "\n", " for row, bra_state in enumerate(basis):\n", " for col, ket_state in enumerate(basis):\n", " h_single[row, col] = _basis_state_transition_amplitude_sparse(\n", " hamiltonian,\n", " bra_state=bra_state,\n", " ket_state=ket_state,\n", " )\n", " # Remove floating-point-level asymmetry.\n", " h_single = 0.5 * (h_single + h_single.conj().T)\n", " evals, _ = np.linalg.eigh(h_single)\n", " return np.real(evals)\n", "\n", "\n", "def simple_transpilation(circuit: QuantumCircuit) -> QuantumCircuit:\n", " \"\"\"Transpilation to simplify the circuit\"\"\"\n", " pm = PassManager(\n", " [\n", " CommutativeOptimization(),\n", " ]\n", " )\n", " circuit = transpile(circuit, optimization_level=3)\n", " circuit = pm.run(circuit)\n", " return circuit\n", "\n", "\n", "def summarize_circuit(circuit: QuantumCircuit) -> dict[str, int | str]:\n", " \"\"\"Summarize the circuit with depth, size, and 2-qubit gate information.\"\"\"\n", " two_qubit_total = sum(\n", " inst.operation.num_qubits == 2 for inst in circuit.data\n", " )\n", " two_qubit_depth = circuit.depth(lambda x: x[0].num_qubits == 2)\n", " return {\n", " \"depth\": circuit.depth(),\n", " \"size\": circuit.size(),\n", " \"2q gates\": two_qubit_total,\n", " \"2q depth\": two_qubit_depth,\n", " }\n", "\n", "\n", "def solve_thresholded_gevp(\n", " h_matrix: np.ndarray,\n", " s_matrix: np.ndarray,\n", " threshold: float = 1e-10,\n", ") -> tuple[float, np.ndarray, int]:\n", " \"\"\"Solve H c = E S c using canonical orthogonalization of S.\"\"\"\n", " s_vals, s_vecs = la.eigh(s_matrix)\n", "\n", " valid = s_vals > threshold\n", " if not np.any(valid):\n", " raise ValueError(\n", " \"All overlap eigenvalues were removed by thresholding.\"\n", " )\n", "\n", " keep = valid\n", "\n", " orthogonalizer = s_vecs[:, keep] @ np.diag(1.0 / np.sqrt(s_vals[keep]))\n", " h_orth = orthogonalizer.conj().T @ h_matrix @ orthogonalizer\n", " h_orth = 0.5 * (h_orth + h_orth.conj().T)\n", "\n", " eigvals, eigvecs = la.eigh(h_orth)\n", "\n", " coeffs = orthogonalizer @ eigvecs[:, 0]\n", " normalization = np.sqrt(np.real(coeffs.conj().T @ s_matrix @ coeffs))\n", " coeffs /= normalization\n", "\n", " return float(np.real(eigvals[0])), coeffs, int(np.sum(keep))" ] }, { "cell_type": "markdown", "id": "7a03a282", "metadata": {}, "source": [ "In the first part of this tutorial, we demonstrate the KQD method by using a local statevector simulator. Later, we use a real quantum backend to address a utility-scale problem.\n", "\n", "We also define a fake backend to demonstrate backend-specific transpilation and inspect the resulting circuit.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8026556a", "metadata": {}, "outputs": [], "source": [ "try:\n", " service = QiskitRuntimeService()\n", "except Exception:\n", " QiskitRuntimeService.save_account(\n", " token=\"\", instance=\"\", overwrite=True\n", " )\n", " service = QiskitRuntimeService()\n", "\n", "backend = FakeMarrakesh()" ] }, { "cell_type": "markdown", "id": "c23d1b80", "metadata": {}, "source": [ "## Small-scale simulator example\n", "\n" ] }, { "cell_type": "markdown", "id": "af6d069a", "metadata": {}, "source": [ "### Step 1: Map classical inputs to a quantum problem\n", "\n", "#### Hamiltonian and reference state\n", "\n", "This example uses a 12-qubit open-boundary Heisenberg chain ($n=12$),\n", "\n", "$$\n", "\\begin{equation*}\n", "H=\\sum_{i=0}^{n-2}\\left(X_iX_{i+1}+Y_iY_{i+1}+Z_iZ_{i+1}\\right),\n", "\\end{equation*}\n", "$$\n", "\n", "with a single-excitation product state\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\psi_{0}\\rangle=|000001000000\\rangle\n", "\\end{equation*}\n", "$$\n", "\n", "as the reference state.\n", "Because the Heisenberg Hamiltonian defined above conserves the total number of excitations, the reference state stays in the single-excitation subspace, whose dimension grows only linearly with the number of qubits. We can therefore compute the exact ground-state energy efficiently, by diagonalizing the Hamiltonian restricted to that subspace, and use it purely as a diagnostic benchmark for the KQD estimate. The KQD workflow itself estimates projected matrix elements by using Qiskit primitives and solves the resulting projected problem classically.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "e156e4c7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hamiltonian:\n", "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'],\n", " 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,\n", " 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,\n", " 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,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "Reference state: |000001000000>\n", "Reference energy: (7+0j)\n" ] } ], "source": [ "# Problem definition for the simulator example\n", "num_qubits = 12\n", "hamiltonian = make_heisenberg_hamiltonian(num_qubits=num_qubits, coupling=1.0)\n", "ref_bitstring = \"000001000000\"\n", "ref_energy = basis_state_expectation_sparse(hamiltonian, ref_bitstring)\n", "\n", "print(\"Hamiltonian:\")\n", "print(hamiltonian)\n", "print(f\"Reference state: |{ref_bitstring}>\")\n", "print(\"Reference energy: \", ref_energy)" ] }, { "cell_type": "markdown", "id": "7f874d2a", "metadata": {}, "source": [ "#### Set the parameters for the algorithm\n", "\n", "Based on the upper-bounds on the Hamiltonian norm, Ref. [\\[1\\]](#references) heuristically suggests the time step $\\Delta t$ as $\\pi/\\|H\\|$.\n", "Because the spectral norm $\\|H\\|$ is hard to compute, we instead use its upper bound:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\|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).\n", "\\end{equation*}\n", "$$\n", "\n", "We set the Krylov dimension to $r=10$ and the number of Trotter steps per time step to $5$: a large enough Krylov space to resolve the low-lying spectrum while keeping the deepest circuit ($t_{\\max}=(r-1)\\Delta t$) affordable, and enough Trotter steps to keep the discretization error small at that deepest circuit.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "af9a86df", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dt in Krylov basis: 0.09519977738150888\n" ] } ], "source": [ "dt = np.pi / (3 * (num_qubits - 1))\n", "print(\"dt in Krylov basis: \", dt)\n", "\n", "krylov_dim = 10\n", "num_trotter_steps = 5" ] }, { "cell_type": "markdown", "id": "5f02e4f3", "metadata": {}, "source": [ "#### Build circuit\n", "\n", "Here, we build the circuits to estimate the matrix elements $\\mathcal{H}_{k\\ell}$ and $\\mathcal{S}_{k\\ell}$.\n", "Because all powers of $H$ commute, we have\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{H}_{k\\ell}\n", "=\n", "\\langle\\psi_{0}|H e^{-i(\\ell-k)\\Delta tH}|\\psi_{0}\\rangle=\\mathcal{H}_{0,\\ell-k},\n", "\\qquad\n", "\\mathcal{S}_{k\\ell}\n", "=\n", "\\langle\\psi_{0}|e^{-i(\\ell-k)\\Delta tH}|\\psi_{0}\\rangle=\\mathcal{S}_{0,\\ell-k}.\n", "\\end{equation*}\n", "$$\n", "\n", "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$.\n", "\n", "Here, we present the circuit, called **extended-swap-test**, which prepares\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\Phi_{0d}\\rangle=\n", "\\frac{|0\\rangle|\\psi_0\\rangle+|1\\rangle|\\psi_d\\rangle}{\\sqrt{2}},\n", "\\end{equation*}\n", "$$\n", "\n", "where $|\\psi_d\\rangle=e^{-id\\Delta tH}|\\psi_{0}\\rangle$.\n", "\n", "#### Reference state\n", "\n", "We prepare the reference state $|\\psi_0\\rangle$.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "fb2879c5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "qc_ref = QuantumCircuit(num_qubits)\n", "for i, b in enumerate(reversed(ref_bitstring)):\n", " if b == \"1\":\n", " qc_ref.x(i)\n", "display(qc_ref.draw(\"mpl\", scale=0.5))" ] }, { "cell_type": "markdown", "id": "4a96258c", "metadata": {}, "source": [ "#### Time evolution\n", "\n", "Realize the time-evolution operator generated by the Hamiltonian, approximated by simple [Lie-Trotterization](/docs/api/qiskit/qiskit.synthesis.LieTrotter).\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "017b8548", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "t = Parameter(\"t\")\n", "\n", "evol_gate = PauliEvolutionGate(\n", " hamiltonian,\n", " time=t,\n", " synthesis=LieTrotter(reps=num_trotter_steps),\n", " label=\"U(t)\",\n", ")\n", "\n", "# Synthesize U(t) first and then control the synthesized circuit.\n", "# This makes the controlled structure visible in the circuit drawer.\n", "evolution_circuit = QuantumCircuit(num_qubits, name=\"U(t)\")\n", "evolution_circuit.append(evol_gate, range(num_qubits))\n", "evolution_circuit = simple_transpilation(evolution_circuit)\n", "\n", "# Make a controlled version of the evolution circuit.\n", "controlled_evolution_gate = evolution_circuit.to_gate(label=\"U(t)\").control(\n", " 1, label=\"C-U(t)\"\n", ")\n", "display(\n", " evolution_circuit.assign_parameters({t: 1.5}).draw(\n", " \"mpl\", scale=0.5, fold=-1\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "efa6076d", "metadata": {}, "source": [ "#### Extended swap test circuit [\\[3\\]](#references)\n", "\n", "The circuit first prepares the reference state on the system register while the ancilla remains in $|0\\rangle$:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle |0\\rangle^{\\otimes n}\n", "\\longrightarrow\n", "|0\\rangle |\\psi_{0}\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Then, applying a Hadamard gate to the ancilla creates a coherent superposition of two branches:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle |\\psi_0\\rangle\n", "\\longrightarrow\n", "\\frac{|0\\rangle + |1\\rangle}{\\sqrt{2}} |\\psi_0\\rangle\n", "=\n", "\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Finally, the controlled time-evolution gate applies:\n", "\n", "$$\n", "\\begin{equation*}\n", "U(t=d\\Delta t) = e^{-id\\Delta t H}\n", "\\end{equation*}\n", "$$\n", "\n", "only when the ancilla is in the $|1\\rangle$ branch. Therefore,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle|\\psi_0\\rangle}{\\sqrt{2}}\n", "\\longrightarrow\n", "\\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle U(d\\Delta t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "In the next code block, we implement:\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\Psi(t)\\rangle = \\frac{|0\\rangle|\\psi_0\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}},\n", "\\end{equation*}\n", "$$\n", "\n", "which will be assigned as $t=d\\Delta t$ for $d=1,\\cdots r-1$ in the execution step.\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "1260f66a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ancilla = 0\n", "system_qubits = list(range(1, num_qubits + 1))\n", "\n", "extended_swap_test = QuantumCircuit(num_qubits + 1)\n", "\n", "# Append state preparation part\n", "extended_swap_test = extended_swap_test.compose(qc_ref, system_qubits)\n", "\n", "# Prepare the coherent branch label, (|0> + |1>) / sqrt(2).\n", "extended_swap_test.h(ancilla)\n", "\n", "# Apply U(t) only to the |1> branch of the ancilla.\n", "extended_swap_test.append(\n", " controlled_evolution_gate, [ancilla] + system_qubits\n", ")\n", "\n", "# Decompose once more for visualization so that control bullets are visible.\n", "display(extended_swap_test.draw(\"mpl\", fold=-1))" ] }, { "cell_type": "markdown", "id": "b05e2e7d", "metadata": {}, "source": [ "#### Observables\n", "\n", "For any Hermitian system observable $O$, here we set the observables to calculate:\n", "\n", "$$\n", "\\begin{equation*}\n", "z=\\langle \\psi_0|O|\\psi_d\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "This is because $O=I$ gives the overlap element $\\mathcal{S}_{0d}=\\langle\\psi_0|\\psi_d\\rangle$, while $O=H$ gives the Hamiltonian element $\\mathcal{H}_{0d}=\\langle\\psi_0|H|\\psi_d\\rangle$.\n", "\n", "Using $X=|0\\rangle\\langle 1|+|1\\rangle\\langle 0|$, we have\n", "\n", "$$\n", "\\begin{align*}\n", "\\langle \\Psi_{0d}|X\\otimes O| \\Psi_{0d} \\rangle\n", "&=\n", "\\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) \\\\\n", "&=\n", "\\frac{1}{2}\n", "\\left(\n", "\\langle \\psi_0|O|\\psi_d\\rangle\n", "+\n", "\\langle \\psi_d|O|\\psi_0\\rangle\n", "\\right) \\\\\n", "&=\n", "\\frac{z+z^*}{2}\n", "=\n", "\\operatorname{Re} z.\n", "\\end{align*}\n", "$$\n", "\n", "Similarly, using $Y=-i|0\\rangle\\langle 1|+i|1\\rangle\\langle 0|$,\n", "\n", "$$\n", "\\begin{align*}\n", "\\langle \\Psi_{0d}|Y\\otimes O| \\Psi_{0d} \\rangle\n", "&=\n", "\\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) \\\\\n", "&=\n", "\\frac{1}{2}\n", "\\left(\n", "-i\\langle \\psi_0|O|\\psi_d\\rangle\n", "+\n", "i\\langle \\psi_d|O|\\psi_0\\rangle\n", "\\right) \\\\\n", "&=\n", "\\frac{-iz+iz^*}{2}\n", "=\n", "\\operatorname{Im} z.\n", "\\end{align*}\n", "$$\n", "\n", "Therefore, we have\n", "\n", "$$\n", "\\begin{equation*}\n", "\\langle \\Phi_{0d}|X\\otimes O|\\Phi_{0d}\\rangle\n", "=\n", "\\operatorname{Re}\\langle\\psi_0|O|\\psi_d\\rangle ,\n", "\\quad\n", "\\langle \\Phi_{0d}|Y\\otimes O|\\Phi_{0d}\\rangle\n", "=\n", "\\operatorname{Im}\\langle\\psi_0|O|\\psi_d\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Finally, for each state $|\\Psi_{0d}\\rangle$, we need to measure:\n", "\n", "$$\n", "\\begin{align*}\n", "X \\otimes I \\text{ for } \\operatorname{Re}\\mathcal{S}_{0d},\\\\\n", "Y \\otimes I \\text{ for } \\operatorname{Im}\\mathcal{S}_{0d},\\\\\n", "X \\otimes H \\text{ for } \\operatorname{Re}\\mathcal{H}_{0d},\\\\\n", "Y \\otimes H \\text{ for } \\operatorname{Im}\\mathcal{H}_{0d}.\\\\\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "4c563534", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Observable: Re S_0d\n", "SparsePauliOp(['IIIIIIIIIIIIX'],\n", " coeffs=[1.+0.j])\n", "\n", "Observable: Im S_0d\n", "SparsePauliOp(['IIIIIIIIIIIIY'],\n", " coeffs=[1.+0.j])\n", "\n", "Observable: Re H_0d\n", "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'],\n", " 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,\n", " 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,\n", " 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,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n", "Observable: Im H_0d\n", "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'],\n", " 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,\n", " 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,\n", " 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,\n", " 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n" ] } ], "source": [ "n_qubits = hamiltonian.num_qubits\n", "\n", "observable_labels = [\n", " \"Re S_0d\",\n", " \"Im S_0d\",\n", " \"Re H_0d\",\n", " \"Im H_0d\",\n", "]\n", "\n", "# X ⊗ I and Y ⊗ I.\n", "# Qiskit's Pauli-label convention places qubit 0 on the rightmost character,\n", "# so the ancilla Pauli is appended to the right.\n", "obs_x_identity = SparsePauliOp(\"I\" * n_qubits + \"X\")\n", "obs_y_identity = SparsePauliOp(\"I\" * n_qubits + \"Y\")\n", "\n", "# X ⊗ H and Y ⊗ H.\n", "obs_x_hamiltonian = SparsePauliOp.from_list(\n", " [\n", " (label + \"X\", coeff)\n", " for label, coeff in zip(\n", " hamiltonian.paulis.to_labels(),\n", " hamiltonian.coeffs,\n", " )\n", " ]\n", ")\n", "\n", "obs_y_hamiltonian = SparsePauliOp.from_list(\n", " [\n", " (label + \"Y\", coeff)\n", " for label, coeff in zip(\n", " hamiltonian.paulis.to_labels(),\n", " hamiltonian.coeffs,\n", " )\n", " ]\n", ")\n", "\n", "observables = [\n", " obs_x_identity,\n", " obs_y_identity,\n", " obs_x_hamiltonian,\n", " obs_y_hamiltonian,\n", "]\n", "\n", "for obs, label in zip(observables, observable_labels):\n", " print(f\"Observable: {label}\")\n", " print(obs)\n", " print()" ] }, { "cell_type": "markdown", "id": "780486da", "metadata": {}, "source": [ "For the estimation of the $\\mathcal{H}$ matrix elements, the number of Pauli terms is much larger than that for the $\\mathcal{S}$ matrix elements.\n", "\n", "We can now reduce the number of measured Hamiltonian terms by using a shifting technique [\\[4\\]](#references).\n", "We split the Hamiltonian as:\n", "\n", "$$\n", "\\begin{equation*}\n", "H = (H-T) + T,\n", "\\end{equation*}\n", "$$\n", "\n", "where $T$ is chosen such that the reference state is its eigenstate,\n", "\n", "$$\n", "\\begin{equation*}\n", "T|\\psi_0\\rangle = \\tau |\\psi_0\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Then,\n", "\n", "$$\n", "\\begin{align*}\n", "\\mathcal{H}_{0d}\n", "&=\n", "\\langle\\psi_0|H|\\psi_d\\rangle \\\\\n", "&=\n", "\\langle\\psi_0|(H-T)|\\psi_d\\rangle\n", "+\n", "\\langle\\psi_0|T|\\psi_d\\rangle \\\\\n", "&=\n", "\\langle\\psi_0|(H-T)|\\psi_d\\rangle\n", "+\n", "\\tau \\langle\\psi_0|\\psi_d\\rangle \\\\\n", "&=\n", "\\widetilde{\\mathcal{H}}_{0d}\n", "+\n", "\\tau \\mathcal{S}_{0d}.\n", "\\end{align*}\n", "$$\n", "\n", "Here,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\widetilde{\\mathcal{H}}_{0d}\n", "=\n", "\\langle\\psi_0|(H-T)|\\psi_d\\rangle\n", "\\end{equation*}\n", "$$\n", "\n", "is the shifted Hamiltonian matrix element.\n", "Therefore, we only need to measure $X\\otimes (H-T)$ and $Y\\otimes (H-T)$.\n", "The contribution from $T$ is reconstructed classically using the already measured overlap matrix element $\\mathcal{S}_{0d}$.\n", "\n", "In this example, a natural choice is the diagonal part of the Heisenberg Hamiltonian,\n", "\n", "$$\n", "\\begin{equation*}\n", "T\n", "=\n", "\\sum_{i=0}^{n-2} Z_iZ_{i+1}.\n", "\\end{equation*}\n", "$$\n", "\n", "Because the reference state is a computational-basis state, it is an eigenstate of every $Z_iZ_{i+1}$ term.\n", "\n", "However, a more advantageous choice is to include not only the diagonal $ZZ$ terms, but also the $XX+YY$ terms that annihilate the reference state.\n", "\n", "For each neighboring pair, the operator $XX+YY$ satisfies\n", "\n", "$$\n", "\\begin{equation*}\n", "(XX+YY)|00\\rangle = 0,\n", "\\qquad\n", "(XX+YY)|11\\rangle = 0,\n", "\\end{equation*}\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\begin{equation*}\n", "(XX+YY)|01\\rangle = 2|10\\rangle,\n", "\\qquad\n", "(XX+YY)|10\\rangle = 2|01\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Therefore, the $XX+YY$ term contributes only when the two neighboring qubits have different occupations in the reference bitstring.\n", "If the two qubits are both $0$ or both $1$, the term annihilates the reference state and can also be shifted out.\n", "\n", "Let $|\\psi_{\\rm ref}\\rangle = |z_0\\cdots z_{n-1}\\rangle$, where $z_i \\in \\{0,1\\}$.\n", "Thus, we can choose\n", "\n", "$$\n", "\\begin{equation*}\n", "T=\\sum_{i=0}^{n-2} Z_iZ_{i+1}\n", "+\n", "\\sum_{z_i={z_{i+1}}}\n", "\\left(\n", "X_iX_{i+1}\n", "+\n", "Y_iY_{i+1}\n", "\\right).\n", "\\end{equation*}\n", "$$\n", "\n", "This operator still satisfies\n", "\n", "$$\n", "\\begin{equation*}\n", "T|\\psi_0\\rangle = \\tau |\\psi_0\\rangle ,\n", "\\end{equation*}\n", "$$\n", "\n", "because the $ZZ$ terms act diagonally on $|\\psi_0\\rangle$, while the shifted $XX+YY$ terms give zero.\n", "The corresponding eigenvalue is therefore determined only by the $ZZ$ terms,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\tau = \\sum_{i=0}^{n-2}\n", "(-1)^{z_i}\n", "(-1)^{z_{i+1}} .\n", "\\end{equation*}\n", "$$\n", "\n", "With this choice, the shifted Hamiltonian becomes:\n", "\n", "$$\n", "\\begin{equation*}\n", "H-T = \\sum_{z_i\\ne z_{i+1}}\n", "\\left(\n", "X_iX_{i+1}\n", "+\n", "Y_iY_{i+1}\n", "\\right).\n", "\\end{equation*}\n", "$$\n", "\n", "As a result, only the edges with different occupations in the reference state need to be measured.\n", "All $ZZ$ terms and all inactive $XX+YY$ terms are reconstructed through the overlap contribution $\\tau \\mathcal{S}_{0d}$, or give zero contribution by construction.\n", "\n", "This gives a smaller observable than shifting only the diagonal part.\n", "In particular, for a computational-basis reference state with a localized excitation, only the edges adjacent to the excitation remain in $H-T$.\n", "Therefore, the number of Pauli terms in $X\\otimes(H-T)$ and $Y\\otimes(H-T)$ can be substantially reduced, while the reconstructed matrix element\n", "\n", "$$\n", "\\begin{equation*}\n", "\\widetilde{\\mathcal{H}}_{0d}\n", "+\n", "\\tau \\mathcal{S}_{0d}\n", "\\end{equation*}\n", "$$\n", "\n", "remains exactly the same.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "52a9b116", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Observable: Re shifted H_0d\n", "SparsePauliOp(['IIIIIXXIIIIIX', 'IIIIIYYIIIIIX', 'IIIIXXIIIIIIX', 'IIIIYYIIIIIIX'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n", "Observable: Im shifted H_0d\n", "SparsePauliOp(['IIIIIXXIIIIIY', 'IIIIIYYIIIIIY', 'IIIIXXIIIIIIY', 'IIIIYYIIIIIIY'],\n", " coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])\n", "\n", "tau = 7.0\n" ] } ], "source": [ "def make_reduced_heisenberg_observables(\n", " ref_bitstring: str,\n", " coupling: float = 1.0,\n", ") -> tuple[SparsePauliOp, SparsePauliOp, float]:\n", " \"\"\"Build X⊗(H-T), Y⊗(H-T), and tau.\"\"\"\n", " n_qubits = len(ref_bitstring)\n", "\n", " shifted_terms: list[tuple[str, complex]] = []\n", " tau = 0.0\n", "\n", " def bit(q: int) -> str:\n", " return ref_bitstring[n_qubits - 1 - q]\n", "\n", " def append_term(q0: int, q1: int, pauli: str):\n", " label = [\"I\"] * n_qubits\n", " label[n_qubits - 1 - q0] = pauli[0]\n", " label[n_qubits - 1 - q1] = pauli[1]\n", " shifted_terms.append((\"\".join(label), coupling))\n", "\n", " for q in range(n_qubits - 1):\n", " same_occupation = bit(q) == bit(q + 1)\n", "\n", " # ZZ contribution to tau\n", " tau += coupling * (1.0 if same_occupation else -1.0)\n", "\n", " # XX + YY survives only for opposite occupations.\n", " if not same_occupation:\n", " append_term(q, q + 1, \"XX\")\n", " append_term(q, q + 1, \"YY\")\n", "\n", " if shifted_terms:\n", " obs_x_shifted_hamiltonian = SparsePauliOp.from_list(\n", " [(label + \"X\", coeff) for label, coeff in shifted_terms]\n", " )\n", " obs_y_shifted_hamiltonian = SparsePauliOp.from_list(\n", " [(label + \"Y\", coeff) for label, coeff in shifted_terms]\n", " )\n", " else:\n", " obs_x_shifted_hamiltonian = SparsePauliOp(\n", " \"I\" * n_qubits + \"X\", coeffs=[0.0]\n", " )\n", " obs_y_shifted_hamiltonian = SparsePauliOp(\n", " \"I\" * n_qubits + \"Y\", coeffs=[0.0]\n", " )\n", "\n", " return obs_x_shifted_hamiltonian, obs_y_shifted_hamiltonian, tau\n", "\n", "\n", "obs_x_shifted_hamiltonian, obs_y_shifted_hamiltonian, shift_tau = (\n", " make_reduced_heisenberg_observables(ref_bitstring)\n", ")\n", "\n", "print(\"Observable: Re shifted H_0d\")\n", "print(obs_x_shifted_hamiltonian)\n", "print()\n", "\n", "print(\"Observable: Im shifted H_0d\")\n", "print(obs_y_shifted_hamiltonian)\n", "print()\n", "\n", "print(\"tau =\", shift_tau)\n", "\n", "observables = [\n", " obs_x_identity,\n", " obs_y_identity,\n", " obs_x_shifted_hamiltonian,\n", " obs_y_shifted_hamiltonian,\n", "]\n", "observable_labels = [\n", " \"Re S_0d\",\n", " \"Im S_0d\",\n", " \"Re shifted H_0d\",\n", " \"Im shifted H_0d\",\n", "]" ] }, { "cell_type": "markdown", "id": "7c7daa55", "metadata": {}, "source": [ "### Step 2: Optimize problem for quantum hardware execution\n", "\n", "Now we turn the abstract extended-swap-test circuit into a hardware-oriented template.\n", "Before that, we further optimize the circuit in the abstract level first.\n", "\n", "#### Compare Hamiltonian term ordering\n", "\n", "First, we compare different orderings of the Pauli terms in the Heisenberg Hamiltonian for the Hamiltonian simulation.\n", "The Hamiltonian itself is unchanged, but the ordering affects how the product-formula circuit is generated and how much the structure can be parallelized in the circuit.\n", "For example, the naive ordering lists all nearest-neighbor $XX$ terms, then all $YY$ terms, then all $ZZ$ terms.\n", "This places adjacent edges such as $(0,1)$ and $(1,2)$ next to each other, so they cannot be executed in parallel.\n", "The even-then-odd ordering visits disjoint even edges first, followed by odd edges, which exposes parallel two-qubit layers.\n", "The even-odd edge-grouped ordering goes one step further: for each edge, it keeps the local $XX$, $YY$, and $ZZ$ terms together, while still visiting even edges before odd edges.\n", "We expect the even-then-odd and even-odd edge-grouped orderings to reduce circuit depth by exposing parallel two-qubit layers, and the edge-grouped ordering to additionally reduce the Trotter error because the local two-qubit interaction on the same edge is handled as a compact block.\n", "\n", "![Pauli evolutions with different term ordering](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-quantum-diagonalization/ham_ordering.avif)\n", "\n", "The Hamiltonian ordering also affects the Trotter error.\n", "If we place noncommuting terms next to each other, the basis transitions occur more often, which induces more Trotter error.\n", "By grouping terms that require the same Pauli-basis transformation, redundant basis-changes can be avoided.\n", "\n", "Here, the comparison uses the largest evolution time that appears in the first-row Krylov estimates,\n", "$t_{\\max} = (r-1)\\Delta t,$\n", "using the same transpilation condition.\n", "\n", "In order to measure the Trotter error, we use **process infidelity** between the Trotterized circuit and the exact Hamiltonian evolution,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\text{Infidelity}(U,V)\n", "=\n", "1-F(U,V)\n", "=\n", "1 - \\frac{|\\operatorname{Tr}(U^\\dagger V)|^2}{d^2},\n", "\\end{equation*}\n", "$$\n", "\n", "where $d=2^n$ is the dimension of the Hilbert space.\n", "\n", "This diagnostic uses dense matrices, so it is suitable for this small 12-qubit example but not intended as a scalable subroutine.\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "d21eae64", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Comparison time: 0.8567979964335799\n", "\n", "naive:\n", "['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']\n", "\n", "even-then-odd:\n", "['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']\n", "\n", "even-odd edge-grouped:\n", "['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']\n", "\n", "Precomputing the exact evolution operator... Done\n" ] } ], "source": [ "# The comparison uses the largest time that appears in the first-row Krylov estimates.\n", "# Circuit depth does not depend on this numeric value, but the Trotter error does.\n", "comparison_time = (krylov_dim - 1) * dt\n", "\n", "\n", "def make_heisenberg_hamiltonian_ordered(\n", " num_qubits: int,\n", " ordering: str,\n", " coupling: float = 1.0,\n", ") -> SparsePauliOp:\n", " \"\"\"Return the same Heisenberg Hamiltonian with a specified term ordering.\"\"\"\n", " terms: list[tuple[str, complex]] = []\n", " even_edges = [(q, q + 1) for q in range(0, num_qubits - 1, 2)]\n", " odd_edges = [(q, q + 1) for q in range(1, num_qubits - 1, 2)]\n", "\n", " def append_term(q0: int, q1: int, pauli: str):\n", " label = [\"I\"] * num_qubits\n", " label[num_qubits - 1 - q0] = pauli[0]\n", " label[num_qubits - 1 - q1] = pauli[1]\n", " terms.append((\"\".join(label), coupling))\n", "\n", " if ordering == \"naive\":\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " for q in range(num_qubits - 1):\n", " append_term(q, q + 1, pauli)\n", " elif ordering == \"even-then-odd\":\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " for q0, q1 in even_edges + odd_edges:\n", " append_term(q0, q1, pauli)\n", " elif ordering == \"even-odd edge-grouped\":\n", " for q0, q1 in even_edges + odd_edges:\n", " for pauli in (\"XX\", \"YY\", \"ZZ\"):\n", " append_term(q0, q1, pauli)\n", " else:\n", " raise ValueError(f\"Unknown ordering: {ordering}\")\n", "\n", " return SparsePauliOp.from_list(terms).simplify()\n", "\n", "\n", "def build_numeric_evolution_circuit(\n", " hamiltonian: SparsePauliOp,\n", " synthesis,\n", " time_value: float,\n", " **synthesis_kwargs,\n", ") -> QuantumCircuit:\n", " \"\"\"Build a numeric circuit for exp(-i H t) with a chosen synthesis rule.\"\"\"\n", " evolution_gate = PauliEvolutionGate(\n", " hamiltonian,\n", " time=time_value,\n", " synthesis=synthesis(**synthesis_kwargs),\n", " )\n", " circuit = QuantumCircuit(hamiltonian.num_qubits)\n", " circuit.append(evolution_gate, range(hamiltonian.num_qubits))\n", " return circuit\n", "\n", "\n", "def process_infidelity(\n", " circuit: QuantumCircuit,\n", " exact_matrix: np.ndarray,\n", ") -> float:\n", " \"\"\"Return 1 - |Tr(U_circuit† U_exact) / d|².\"\"\"\n", " circuit_matrix = np.asarray(Operator(circuit).data)\n", " dim = circuit_matrix.shape[0]\n", " normalized_trace = np.vdot(circuit_matrix, exact_matrix) / dim\n", " fidelity = np.abs(normalized_trace) ** 2\n", " return float(np.clip(1.0 - fidelity, 0.0, 1.0))\n", "\n", "\n", "hamiltonians_by_ordering = {\n", " ordering: make_heisenberg_hamiltonian_ordered(\n", " num_qubits, ordering, coupling=1.0\n", " )\n", " for ordering in [\"naive\", \"even-then-odd\", \"even-odd edge-grouped\"]\n", "}\n", "\n", "print(f\"Comparison time: {comparison_time}\\n\")\n", "\n", "for order_name, ham_ordered in hamiltonians_by_ordering.items():\n", " print(f\"{order_name}:\")\n", " print([op for op, _ in ham_ordered.to_list()])\n", " print()\n", "\n", "print(\"Precomputing the exact evolution operator... \", end=\"\")\n", "exact_matrix = la.expm(-1j * comparison_time * hamiltonian.to_matrix())\n", "print(\"Done\")" ] }, { "cell_type": "markdown", "id": "d3fcb5e1", "metadata": {}, "source": [ "We first keep the synthesis rule fixed to one first-order Trotter step and vary only the Pauli-term ordering.\n", "The goal of this comparison is mainly to see how much circuit depth and two-qubit cost can be reduced by exposing disjoint nearest-neighbor edges to the transpiler.\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "21eaab38", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Circuit for naive ordering:\n", "['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']\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Circuit for even-then-odd ordering:\n", "['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']\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Circuit for even-odd edge-grouped ordering:\n", "['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']\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
orderingsynthesisinfidelitydepthsize2q gates2q depth
0naiveLieTrotter(reps=1)0.99991715333315
1naiveLieTrotter(reps=2)0.80599821666621
2naiveLieTrotter(reps=4)0.2713893313213233
3naiveLieTrotter(reps=8)0.0745905726426457
4even-then-oddLieTrotter(reps=1)0.999917633336
5even-then-oddLieTrotter(reps=2)0.80599812666612
6even-then-oddLieTrotter(reps=4)0.2713892413213224
7even-then-oddLieTrotter(reps=8)0.0745904826426448
8even-odd edge-groupedLieTrotter(reps=1)0.998432633336
9even-odd edge-groupedLieTrotter(reps=2)0.65384312666612
10even-odd edge-groupedLieTrotter(reps=4)0.1815292413213224
11even-odd edge-groupedLieTrotter(reps=8)0.0452444826426448
\n", "
" ], "text/plain": [ " ordering synthesis infidelity depth size \\\n", "0 naive LieTrotter(reps=1) 0.999917 15 33 \n", "1 naive LieTrotter(reps=2) 0.805998 21 66 \n", "2 naive LieTrotter(reps=4) 0.271389 33 132 \n", "3 naive LieTrotter(reps=8) 0.074590 57 264 \n", "4 even-then-odd LieTrotter(reps=1) 0.999917 6 33 \n", "5 even-then-odd LieTrotter(reps=2) 0.805998 12 66 \n", "6 even-then-odd LieTrotter(reps=4) 0.271389 24 132 \n", "7 even-then-odd LieTrotter(reps=8) 0.074590 48 264 \n", "8 even-odd edge-grouped LieTrotter(reps=1) 0.998432 6 33 \n", "9 even-odd edge-grouped LieTrotter(reps=2) 0.653843 12 66 \n", "10 even-odd edge-grouped LieTrotter(reps=4) 0.181529 24 132 \n", "11 even-odd edge-grouped LieTrotter(reps=8) 0.045244 48 264 \n", "\n", " 2q gates 2q depth \n", "0 33 15 \n", "1 66 21 \n", "2 132 33 \n", "3 264 57 \n", "4 33 6 \n", "5 66 12 \n", "6 132 24 \n", "7 264 48 \n", "8 33 6 \n", "9 66 12 \n", "10 132 24 \n", "11 264 48 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ordering_comparison_rows = []\n", "infidelity_reps = [1, 2, 4, 8]\n", "\n", "for ordering, ham_ordered in hamiltonians_by_ordering.items():\n", " for reps in infidelity_reps:\n", " circuit = build_numeric_evolution_circuit(\n", " ham_ordered,\n", " LieTrotter,\n", " comparison_time,\n", " reps=reps,\n", " )\n", " decomposed_circuit = simple_transpilation(circuit)\n", " infidelity = process_infidelity(decomposed_circuit, exact_matrix)\n", " ordering_comparison_rows.append(\n", " {\n", " \"ordering\": ordering,\n", " \"synthesis\": f\"LieTrotter(reps={reps})\",\n", " \"infidelity\": infidelity,\n", " **summarize_circuit(decomposed_circuit),\n", " }\n", " )\n", " if reps == 1:\n", " print(f\"Circuit for {ordering} ordering:\")\n", " print([op for op, _ in ham_ordered.to_list()])\n", " display(decomposed_circuit.draw(\"mpl\", fold=-1, scale=0.6))\n", "\n", "ordering_comparison_df = pd.DataFrame(ordering_comparison_rows)\n", "display(ordering_comparison_df)\n", "\n", "hamiltonian_for_synthesis = hamiltonians_by_ordering[\"even-odd edge-grouped\"]" ] }, { "cell_type": "markdown", "id": "1c15299c", "metadata": {}, "source": [ "At `reps=1`, we observe that even-then-odd and even-odd edge-grouped orderings both collapse the depth from 15 to 6 by exposing parallel two-qubit layers, while the two-qubit gate count stays the same across all three orderings.\n", "However, all three orderings have infidelity close to 1, so we sweep the number of Trotter repetitions to separate the orderings more clearly.\n", "As the number of repetitions grows, the infidelity of `even-odd edge-grouped` ordering drops faster than the other two, reaching 0.045 at `reps=8` versus 0.075 for naive and even-then-odd ordering.\n", "\n", "#### Compare product-formula synthesis\n", "\n", "Next, we explore different advanced settings of Trotterization, with the fixed Hamiltonian ordering to the even-odd edge-grouped ordering.\n", "We consider first-order Lie-Trotter, second-order Suzuki-Trotter, and fourth-order Suzuki-Trotter.\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "7b6b7f2f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
synthesisrepsinfidelitydepthsize2q gates2q depth
0LieTrotter19.984324e-01633336
1SuzukiTrotter(order=2)19.733399e-01951519
3LieTrotter26.538427e-0112666612
4SuzukiTrotter(order=2)22.522533e-0115848415
6LieTrotter33.197242e-0118999918
7SuzukiTrotter(order=2)34.804050e-022111711721
9LieTrotter41.815291e-012413213224
10SuzukiTrotter(order=2)41.453103e-022715015027
12LieTrotter51.161770e-013016516530
2SuzukiTrotter(order=4)12.884402e-013318318333
13SuzukiTrotter(order=2)55.803455e-033318318333
5SuzukiTrotter(order=4)21.641162e-036334834863
8SuzukiTrotter(order=4)32.907076e-059351351393
11SuzukiTrotter(order=4)42.791061e-06123678678123
14SuzukiTrotter(order=4)54.736685e-07153843843153
\n", "
" ], "text/plain": [ " synthesis reps infidelity depth size 2q gates \\\n", "0 LieTrotter 1 9.984324e-01 6 33 33 \n", "1 SuzukiTrotter(order=2) 1 9.733399e-01 9 51 51 \n", "3 LieTrotter 2 6.538427e-01 12 66 66 \n", "4 SuzukiTrotter(order=2) 2 2.522533e-01 15 84 84 \n", "6 LieTrotter 3 3.197242e-01 18 99 99 \n", "7 SuzukiTrotter(order=2) 3 4.804050e-02 21 117 117 \n", "9 LieTrotter 4 1.815291e-01 24 132 132 \n", "10 SuzukiTrotter(order=2) 4 1.453103e-02 27 150 150 \n", "12 LieTrotter 5 1.161770e-01 30 165 165 \n", "2 SuzukiTrotter(order=4) 1 2.884402e-01 33 183 183 \n", "13 SuzukiTrotter(order=2) 5 5.803455e-03 33 183 183 \n", "5 SuzukiTrotter(order=4) 2 1.641162e-03 63 348 348 \n", "8 SuzukiTrotter(order=4) 3 2.907076e-05 93 513 513 \n", "11 SuzukiTrotter(order=4) 4 2.791061e-06 123 678 678 \n", "14 SuzukiTrotter(order=4) 5 4.736685e-07 153 843 843 \n", "\n", " 2q depth \n", "0 6 \n", "1 9 \n", "3 12 \n", "4 15 \n", "6 18 \n", "7 21 \n", "9 24 \n", "10 27 \n", "12 30 \n", "2 33 \n", "13 33 \n", "5 63 \n", "8 93 \n", "11 123 \n", "14 153 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "synthesis_comparison_rows = []\n", "\n", "for num_trotter_steps in [1, 2, 3, 4, 5]:\n", " synthesis_cases = [\n", " (\"LieTrotter\", LieTrotter, {\"reps\": num_trotter_steps}),\n", " (\n", " \"SuzukiTrotter(order=2)\",\n", " SuzukiTrotter,\n", " {\"order\": 2, \"reps\": num_trotter_steps},\n", " ),\n", " (\n", " \"SuzukiTrotter(order=4)\",\n", " SuzukiTrotter,\n", " {\"order\": 4, \"reps\": num_trotter_steps},\n", " ),\n", " ]\n", " for label, synthesis, kwargs in synthesis_cases:\n", " circuit = build_numeric_evolution_circuit(\n", " hamiltonian_for_synthesis,\n", " synthesis,\n", " comparison_time,\n", " **kwargs,\n", " )\n", " decomposed_circuit = simple_transpilation(circuit)\n", " synthesis_comparison_rows.append(\n", " {\n", " \"synthesis\": label,\n", " \"reps\": kwargs[\"reps\"],\n", " \"infidelity\": process_infidelity(\n", " decomposed_circuit, exact_matrix\n", " ),\n", " **summarize_circuit(decomposed_circuit),\n", " }\n", " )\n", "\n", "synthesis_comparison_df = pd.DataFrame(synthesis_comparison_rows)\n", "display(synthesis_comparison_df.sort_values([\"2q gates\", \"2q depth\"]))\n", "\n", "# For memory free\n", "exact_matrix = None" ] }, { "cell_type": "markdown", "id": "304f3357", "metadata": {}, "source": [ "Lie-Trotter gives the shallowest circuit but has the largest error, while the fourth-order Suzuki-Trotter is more accurate but increases the circuit depth.\n", "For the rest of the tutorial, we choose second-order Suzuki-Trotter because it gives a small-depth circuit while substantially reducing the Trotter error relative to the first-order formula.\n", "\n" ] }, { "cell_type": "markdown", "id": "2415fa62", "metadata": {}, "source": [ "#### Remove the controlled time-evolution gate\n", "\n", "In the extended swap test, controlling time evolution with a single ancilla qubit requires the ancilla to control many gates across the system.\n", "This can introduce substantial routing overhead and, in the worst case, effectively requires all-to-one connectivity.\n", "To avoid this, further optimization is possible by replacing the controlled time-evolution gate with a control-free version, exploiting the symmetry of the Hamiltonian.\n", "Let's observe the following circuit.\n", "\n", "![circuit optimization](https://quantum.cloud.ibm.com/docs/images/tutorials/krylov-quantum-diagonalization/circuit_opt.avif)\n", "\n", "Here, $B_{\\rm ref}$ prepares the reference state, $B_{\\rm ref}|0^n\\rangle = |\\psi_{0}\\rangle$.\n", "\n", "Instead of first preparing $\\frac{1}{\\sqrt{2}}(|0\\rangle+|1\\rangle)|\\psi_0\\rangle$ and then applying $U(t)$ only on the $|1\\rangle$ branch, the circuit directly prepares the two branches as\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle|\\psi_0\\rangle\n", "\\quad \\text{and} \\quad\n", "|1\\rangle|\\psi_d\\rangle ,\n", "\\end{equation*}\n", "$$\n", "\n", "where\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\psi_d\\rangle = U(d\\Delta t)|\\psi_0\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "The circuit first applies a Hadamard gate to the ancilla and prepares the reference state only on the $|1\\rangle$ branch:\n", "\n", "$$\n", "\\begin{equation*}\n", "|0\\rangle |0^n\\rangle\n", "\\longrightarrow\n", "\\frac{|0\\rangle |0^n\\rangle + |1\\rangle B_{\\rm ref}|0^n\\rangle}{\\sqrt{2}}\n", "=\n", "\\frac{|0\\rangle |0^n\\rangle + |1\\rangle |\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Then the uncontrolled time-evolution operator is applied to both branches:\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{|0\\rangle |0^n\\rangle + |1\\rangle |\\psi_0\\rangle}{\\sqrt{2}}\n", "\\longrightarrow\n", "\\frac{|0\\rangle U(t)|0^n\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Since the Hamiltonian preserves the excitation number, we can see that $|0^n\\rangle$ is its eigenstate, and thus\n", "the evolution operator only accumulates a phase under the Hamiltonian:\n", "\n", "$$\n", "\\begin{equation*}\n", "U(t)|0^n\\rangle\n", "=\n", "e^{-iE_{\\rm vac}t}|0^n\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "Therefore,\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{|0\\rangle U(t)|0^n\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}}\n", "=\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |0^n\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Next, $B_{\\rm ref}$ is applied only on the $|0\\rangle$ branch.\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |\\psi_0\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "At this point, the two branches have an extra relative phase. To remove it, we apply the ancilla phase gate\n", "\n", "$$\n", "\\begin{equation*}\n", "P(-E_{\\rm vac}t)\n", "=\n", "\\begin{pmatrix}\n", "1 & 0 \\\\\n", "0 & e^{-iE_{\\rm vac}t}\n", "\\end{pmatrix}.\n", "\\end{equation*}\n", "$$\n", "\n", "This transforms the state as\n", "\n", "$$\n", "\\begin{equation*}\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |\\psi_0\\rangle + |1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}}\n", "\\longrightarrow\n", "\\frac{e^{-iE_{\\rm vac}t}|0\\rangle |\\psi_0\\rangle\n", "+\n", "e^{-iE_{\\rm vac}t}|1\\rangle U(t)|\\psi_0\\rangle}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "Thus, up to an irrelevant global phase, we finally prepared\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\Phi_{0d}\\rangle\n", "=\n", "\\frac{\n", "|0\\rangle|\\psi_0\\rangle\n", "+\n", "|1\\rangle U(t)|\\psi_0\\rangle\n", "}{\\sqrt{2}} .\n", "\\end{equation*}\n", "$$\n", "\n", "In the next code block, we implement this control-free circuit.\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "790e5e0c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Vacuum energy E_vac = 11.0+0.0j\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "controlled_extended_swap_test = extended_swap_test\n", "\n", "# Reuse the ordering and synthesis rule selected by the comparison above.\n", "evol_gate_optimized = PauliEvolutionGate(\n", " hamiltonian_for_synthesis,\n", " time=t,\n", " synthesis=SuzukiTrotter(order=2, reps=num_trotter_steps),\n", ")\n", "\n", "uncontrolled_evolution = QuantumCircuit(num_qubits, name=\"U_ST2(t)\")\n", "uncontrolled_evolution.append(evol_gate_optimized, range(num_qubits))\n", "\n", "controlled_state_prep = QuantumCircuit(num_qubits + 1, name=\"C-Prep\")\n", "for q, bit in enumerate(reversed(ref_bitstring)):\n", " if bit == \"1\":\n", " controlled_state_prep.cx(ancilla, system_qubits[q])\n", "\n", "vacuum_bitstring = \"0\" * num_qubits\n", "vacuum_energy = basis_state_expectation_sparse(hamiltonian, \"0\" * num_qubits)\n", "\n", "optimized_extended_swap_test = QuantumCircuit(num_qubits + 1)\n", "optimized_extended_swap_test.h(ancilla)\n", "\n", "# Prepare |psi_ref> only on the |1> branch.\n", "optimized_extended_swap_test.compose(controlled_state_prep, inplace=True)\n", "optimized_extended_swap_test.barrier()\n", "\n", "# Apply the Trotterized time evolution without control.\n", "optimized_extended_swap_test.compose(\n", " uncontrolled_evolution,\n", " qubits=system_qubits,\n", " inplace=True,\n", ")\n", "optimized_extended_swap_test.barrier()\n", "\n", "# Map |0>|0...0> to |0>|psi_ref>, leaving the |1> branch unchanged.\n", "optimized_extended_swap_test.x(ancilla)\n", "optimized_extended_swap_test.compose(controlled_state_prep, inplace=True)\n", "optimized_extended_swap_test.x(ancilla)\n", "\n", "# Cancel the known vacuum phase so that the same X/Y observables can be used.\n", "optimized_extended_swap_test.p(-vacuum_energy * t, ancilla)\n", "optimized_extended_swap_test = simple_transpilation(\n", " optimized_extended_swap_test\n", ")\n", "\n", "print(f\"Vacuum energy E_vac = {vacuum_energy:.1f}\")\n", "display(\n", " optimized_extended_swap_test.assign_parameters({t: 1.0}).draw(\n", " \"mpl\", scale=0.5, fold=26\n", " )\n", ")" ] }, { "cell_type": "markdown", "id": "e223d1a1", "metadata": {}, "source": [ "#### Transpilation\n", "\n", "Now, we transpile the controlled- and control-free circuits to become executable on the hardware.\n", "Let's compare the result of the transpiled circuits.\n", "\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "6750197b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
labeldepthsize2q gates2q depth
0abstract controlled U(t)154572378646864580
1optimized non-controlled U(t)261171630757
\n", "
" ], "text/plain": [ " label depth size 2q gates 2q depth\n", "0 abstract controlled U(t) 15457 23786 4686 4580\n", "1 optimized non-controlled U(t) 261 1716 307 57" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Layout({\n", "18: ,\n", "5: ,\n", "6: ,\n", "7: ,\n", "8: ,\n", "9: ,\n", "10: ,\n", "11: ,\n", "12: ,\n", "13: ,\n", "14: ,\n", "15: ,\n", "19: \n", "})\n" ] } ], "source": [ "pass_manager = generate_preset_pass_manager(\n", " backend=backend,\n", " optimization_level=3,\n", ")\n", "\n", "isa_controlled_extended_swap_test = pass_manager.run(\n", " controlled_extended_swap_test\n", ")\n", "\n", "pass_manager = generate_preset_pass_manager(\n", " backend=backend, optimization_level=3, routing_method=\"none\"\n", ")\n", "isa_optimized_extended_swap_test = pass_manager.run(\n", " optimized_extended_swap_test\n", ")\n", "\n", "transpilation_result = [\n", " {\n", " \"label\": \"abstract controlled U(t)\",\n", " **summarize_circuit(isa_controlled_extended_swap_test),\n", " },\n", " {\n", " \"label\": \"optimized non-controlled U(t)\",\n", " **summarize_circuit(isa_optimized_extended_swap_test),\n", " },\n", "]\n", "\n", "display(pd.DataFrame(transpilation_result))\n", "\n", "\n", "def filter_qubits_from_layout(layout):\n", " q_layout = Layout(\n", " {\n", " physical: virtual\n", " for physical, virtual in layout.get_physical_bits().items()\n", " if virtual._register.name == \"q\"\n", " }\n", " )\n", " return q_layout\n", "\n", "\n", "print(\n", " filter_qubits_from_layout(\n", " isa_optimized_extended_swap_test.layout.initial_layout\n", " )\n", ")\n", "\n", "isa_observables = [\n", " op.apply_layout(isa_optimized_extended_swap_test.layout)\n", " for op in observables\n", "]" ] }, { "cell_type": "markdown", "id": "39aef2af", "metadata": {}, "source": [ "### Step 3: Execute using Qiskit primitives\n", "\n", "The next step is to submit the same parameterized circuit for several values of $d$. For each $d$, we estimate four expectation values: $X\\otimes I$, $Y\\otimes I$, $X\\otimes (H-T)$, and $Y\\otimes (H-T)$. These four numbers are then combined into the complex first-row elements $\\mathcal{S}_{0d}$ and $\\mathcal{H}_{0d}$.\n", "\n", "Here, $\\mathcal{S}_{00}=1$ and $\\mathcal{H}_{00}=\\langle\\psi_{0}|(H-T)|\\psi_{0}\\rangle=\\langle\\psi_{0}|H|\\psi_{0}\\rangle-\\tau$ can be classically calculated since $\\lvert\\psi_0\\rangle$ is sparse, so we skip the $d=0$ case.\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "83a49027", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of Krylov basis states: r = 9\n", "Number of PUBs: 36\n", "Each d uses observables: ['Re S_0d', 'Im S_0d', 'Re shifted H_0d', 'Im shifted H_0d']\n" ] } ], "source": [ "pub_list = []\n", "d_values = list(range(1, krylov_dim))\n", "\n", "# Exact local statevector estimator.\n", "estimator = StatevectorEstimator()\n", "\n", "# We use the circuit before the transpilation for the local simulator,\n", "# but we will use the transpiled circuit for the real backend.\n", "for d in d_values:\n", " parameter_values = [d * dt]\n", "\n", " for ob in observables:\n", " pub_list.append(\n", " (\n", " optimized_extended_swap_test,\n", " ob,\n", " parameter_values,\n", " )\n", " )\n", "\n", "job = estimator.run(pub_list)\n", "\n", "# Local PrimitiveJob does not provide Runtime-style job inputs,\n", "# so preserve the inputs directly.\n", "inputs = pub_list\n", "result = job.result()\n", "\n", "print(f\"Number of Krylov basis states: r = {len(d_values)}\")\n", "print(f\"Number of PUBs: {len(pub_list)}\")\n", "print(f\"Each d uses observables: {observable_labels}\")" ] }, { "cell_type": "markdown", "id": "a0fd2bc9", "metadata": {}, "source": [ "### Step 4: Post-process and return result in desired classical format\n", "\n", "After estimating the projected matrices, we regularize and solve the GEVP\n", "\n", "$$\n", "\\begin{equation*}\n", "\\mathcal{H}\\mathbf{c} = E\\mathcal{S}\\mathbf{c}.\n", "\\end{equation*}\n", "$$\n", "\n", "The smallest generalized eigenvalue gives the KQD estimate of the ground-state energy.\n", "\n" ] }, { "cell_type": "code", "execution_count": 16, "id": "1019edb3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated first row of S:\n", "[ 1. +0.j 0.758-0.596j 0.203-0.836j -0.291-0.636j -0.444-0.229j\n", " -0.276+0.053j -0.044+0.051j 0.006-0.121j -0.155-0.218j -0.345-0.101j]\n", "\n", "Estimated first row of H:\n", "[ 7. +0.j 4.842-4.76j 0.044-6.185j -3.791-3.653j -4.137+0.39j\n", " -1.495+2.654j 1.331+1.777j 1.85 -0.763j -0.032-2.278j -2.222-1.372j]\n", "\n", "Eigenvalues of the estimated overlap matrix S:\n", "[-0. 0. 0. 0. 0. 0. 0.01 0.355 3.526 6.109]\n", "Condition number above 1e-12: 1.904e+12\n", "\n" ] } ], "source": [ "h_shifted_row_est = np.zeros(krylov_dim, dtype=complex)\n", "s_row_est = np.zeros(krylov_dim, dtype=complex)\n", "\n", "h_shifted_row_est[0] = ref_energy - shift_tau\n", "s_row_est[0] = 1.0\n", "\n", "for idx, (pub_input, pub_result) in enumerate(zip(inputs, result)):\n", " d_index, obs_index = divmod(idx, len(observables))\n", " ev = np.asarray(pub_result.data.evs).reshape(-1)[0]\n", " std = np.asarray(pub_result.data.stds).reshape(-1)[0]\n", "\n", " if obs_index == 0:\n", " s_row_est[d_index + 1] = ev\n", " elif obs_index == 1:\n", " s_row_est[d_index + 1] += 1j * ev\n", " elif obs_index == 2:\n", " h_shifted_row_est[d_index + 1] = ev\n", " elif obs_index == 3:\n", " h_shifted_row_est[d_index + 1] += 1j * ev\n", "\n", "# H_0d = shifted_H_0d + tau * S_0d.\n", "h_row_est = h_shifted_row_est + shift_tau * s_row_est\n", "\n", "h_matrix_est = la.toeplitz(h_row_est.conj(), h_row_est)\n", "s_matrix_est = la.toeplitz(s_row_est.conj(), s_row_est)\n", "\n", "s_eigvals = la.eigvalsh(0.5 * (s_matrix_est + s_matrix_est.conj().T))\n", "positive_s_eigvals = s_eigvals[s_eigvals > 1e-12]\n", "s_condition_number = (\n", " positive_s_eigvals[-1] / positive_s_eigvals[0]\n", " if len(positive_s_eigvals) > 0\n", " else np.inf\n", ")\n", "\n", "\n", "with np.printoptions(precision=3, suppress=True):\n", " print(\"Estimated first row of S:\")\n", " print(s_row_est)\n", " print()\n", " print(\"Estimated first row of H:\")\n", " print(h_row_est)\n", " print()\n", " print(\"Eigenvalues of the estimated overlap matrix S:\")\n", " print(s_eigvals)\n", " print(f\"Condition number above 1e-12: {s_condition_number:.3e}\")\n", " print()" ] }, { "cell_type": "markdown", "id": "176d4283", "metadata": {}, "source": [ "Now we solve the generalized eigenvalue problem using the matrices reconstructed from the circuit estimates.\n", "In an ideal statevector calculation with exact real-time evolution, this should reproduce the exact projected result.\n", "In practice, deviations can come from Trotterization, sampling error, and numerical instability of the overlap matrix.\n", "\n", "We observe how the energy converges as we increase the dimension of the Krylov subspace.\n", "\n" ] }, { "cell_type": "code", "execution_count": 17, "id": "e789b2e0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "exact ground state energy: 3.136296694843727\n", "Krylov ground state energy (dim=1, retained=1): 7.0\n", "Krylov ground state energy (dim=2, retained=2): 4.184510657551266\n", "Krylov ground state energy (dim=3, retained=3): 3.5539074630394136\n", "Krylov ground state energy (dim=4, retained=4): 3.3366270761341044\n", "Krylov ground state energy (dim=5, retained=5): 3.252017453225087\n", "Krylov ground state energy (dim=6, retained=6): 3.2300275138879186\n", "Krylov ground state energy (dim=7, retained=7): 3.2299154099085685\n", "Krylov ground state energy (dim=8, retained=7): 3.2298063744216776\n", "Krylov ground state energy (dim=9, retained=7): 3.2296778282872456\n", "Krylov ground state energy (dim=10, retained=8): 3.223647515867734\n" ] } ], "source": [ "exact_evals = diagonalize_single_1_subspace(hamiltonian)\n", "exact_ground = min(exact_evals)\n", "print(\"exact ground state energy: \", exact_ground)\n", "\n", "threshold = 1e-12\n", "energy_convergence = []\n", "\n", "for r in range(1, krylov_dim + 1):\n", " energy_est_kqd, coeffs_est, retained_est = solve_thresholded_gevp(\n", " h_matrix_est[:r, :r],\n", " s_matrix_est[:r, :r],\n", " threshold=threshold,\n", " )\n", " energy_convergence.append(energy_est_kqd)\n", " print(\n", " f\"Krylov ground state energy (dim={r}, retained={retained_est}): \",\n", " energy_est_kqd,\n", " )" ] }, { "cell_type": "code", "execution_count": 18, "id": "08d57f33", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def plot_energy_convergence(energy_convergence, exact_ground, krylov_dim):\n", " fig, ax = plt.subplots(figsize=(7, 4.5))\n", "\n", " ax.plot(\n", " range(1, krylov_dim + 1),\n", " energy_convergence,\n", " marker=\"o\",\n", " label=\"KQD estimate\",\n", " )\n", "\n", " ax.axhline(\n", " exact_ground,\n", " linestyle=\"--\",\n", " label=f\"Exact ground energy = {exact_ground:.6f}\",\n", " )\n", "\n", " ax.set_xlabel(\"Krylov dimension\")\n", " ax.set_ylabel(\"Ground-state energy\")\n", " ax.grid(True, alpha=0.3)\n", " ax.legend()\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "plot_energy_convergence(energy_convergence, exact_ground, krylov_dim)" ] }, { "cell_type": "markdown", "id": "23c4eb47", "metadata": {}, "source": [ "## Large-scale hardware example\n", "\n", "The previous section used a 12-qubit model so that the statevector simulation could be used as a diagnostic. We now scale the same KQD workflow to a 30-qubit Heisenberg chain and prepare the workload for execution on IBM Quantum hardware.\n", "\n", "### Steps 1-4 compressed into a single code block\n", "\n", "Here we now put all of these details together into a singular workflow at a larger scale, which is then run on our real quantum hardware. In this section, we apply realistic error-mitigation settings to improve the reliability of the results. Since the matrix elements corresponding to different values of $d$ can be evaluated in parallel, we use `Batch` mode to execute them efficiently.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e3e0d1f3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Backend: ibm_boston\n", "System qubits: 30\n", "Total circuit qubits: 156\n", "Krylov dimension: 7\n", "Time step: 0.036110\n", "Shift tau: 25.0\n", "ISA circuit depth: 219\n", "ISA two-qubit gates: 728\n", "ISA two-qubit depth: 52\n", "Job IDs by d: {1: 'd9fl4p4jeosc73fk4dg0', 2: 'd9fl4pineu4c739poecg', 3: 'd9fl4q2neu4c739poedg', 4: 'd9fl4qhhtsac739fjgn0', 5: 'd9fl4r4jeosc73fk4dhg', 6: 'd9fl4rkjeosc73fk4dj0'}\n", "Estimated first row of S:\n", "[ 1. +0.j 0.42055-0.65466j -0.17883-0.73125j -0.64523-0.31601j\n", " -0.70781+0.55574j -0.02467+0.70525j 0.59301+0.55602j]\n", "Estimated first row of H:\n", "[ 25. +0.j 10.3419 -16.5586j -5.017 -18.03933j\n", " -16.44797 -6.98173j -17.17081+14.84502j 0.62722+17.81196j\n", " 15.68886+12.74875j]\n", "exact ground state energy: 21.021912418526902\n", "Krylov ground state energy (dim=1, retained=1): 25.0\n", "Krylov ground state energy (dim=2, retained=2): 24.432409110686205\n", "Krylov ground state energy (dim=3, retained=3): 24.256856893278556\n", "Krylov ground state energy (dim=4, retained=4): 23.727409826799715\n", "Krylov ground state energy (dim=5, retained=4): 23.324720470780864\n", "Krylov ground state energy (dim=6, retained=5): 21.91957579005085\n", "Krylov ground state energy (dim=7, retained=5): 21.548331214122644\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# -------------------------Step 1-------------------------\n", "\n", "# Map the classical problem to quantum circuits and observables.\n", "# Problem and KQD parameters.\n", "large_num_qubits = 30\n", "large_krylov_dim = 7\n", "large_num_trotter_steps = 3\n", "large_dt = np.pi / (3 * (large_num_qubits - 1))\n", "large_t = Parameter(\"t_large\")\n", "\n", "# Use a single excitation near the center of the chain.\n", "large_excitation_qubit = large_num_qubits // 2\n", "large_ref_label = [\"0\"] * large_num_qubits\n", "large_ref_label[large_num_qubits - 1 - large_excitation_qubit] = \"1\"\n", "large_ref_bitstring = \"\".join(large_ref_label)\n", "\n", "# Use the ordering and product formula selected in the preceding section.\n", "large_hamiltonian = make_heisenberg_hamiltonian_ordered(\n", " large_num_qubits,\n", " ordering=\"even-odd edge-grouped\",\n", " coupling=1.0,\n", ")\n", "large_ref_energy = float(\n", " np.real(\n", " basis_state_expectation_sparse(large_hamiltonian, large_ref_bitstring)\n", " )\n", ")\n", "large_vacuum_energy = float(\n", " np.real(\n", " basis_state_expectation_sparse(\n", " large_hamiltonian, \"0\" * large_num_qubits\n", " )\n", " )\n", ")\n", "\n", "large_evolution_gate = PauliEvolutionGate(\n", " large_hamiltonian,\n", " time=large_t,\n", " synthesis=SuzukiTrotter(order=2, reps=large_num_trotter_steps),\n", ")\n", "large_uncontrolled_evolution = QuantumCircuit(\n", " large_num_qubits,\n", " name=\"U_ST2_large(t)\",\n", ")\n", "large_uncontrolled_evolution.append(\n", " large_evolution_gate,\n", " range(large_num_qubits),\n", ")\n", "\n", "# Build the control-free extended-swap-test circuit.\n", "large_ancilla = 0\n", "large_system_qubits = list(range(1, large_num_qubits + 1))\n", "\n", "large_controlled_state_prep = QuantumCircuit(\n", " large_num_qubits + 1,\n", " name=\"C-Prep-large\",\n", ")\n", "large_controlled_state_prep.cx(\n", " large_ancilla,\n", " large_system_qubits[large_excitation_qubit],\n", ")\n", "\n", "large_extended_swap_test = QuantumCircuit(large_num_qubits + 1)\n", "large_extended_swap_test.h(large_ancilla)\n", "large_extended_swap_test.compose(large_controlled_state_prep, inplace=True)\n", "large_extended_swap_test.compose(\n", " large_uncontrolled_evolution,\n", " qubits=large_system_qubits,\n", " inplace=True,\n", ")\n", "large_extended_swap_test.x(large_ancilla)\n", "large_extended_swap_test.compose(\n", " large_controlled_state_prep.inverse(),\n", " inplace=True,\n", ")\n", "large_extended_swap_test.x(large_ancilla)\n", "large_extended_swap_test.p(\n", " -large_vacuum_energy * large_t,\n", " large_ancilla,\n", ")\n", "\n", "# Reuse the Hamiltonian-shifting construction from the preceding section.\n", "(\n", " large_obs_x_shifted_hamiltonian,\n", " large_obs_y_shifted_hamiltonian,\n", " large_shift_tau,\n", ") = make_reduced_heisenberg_observables(large_ref_bitstring)\n", "\n", "large_observables = [\n", " SparsePauliOp(\"I\" * large_num_qubits + \"X\"),\n", " SparsePauliOp(\"I\" * large_num_qubits + \"Y\"),\n", " large_obs_x_shifted_hamiltonian,\n", " large_obs_y_shifted_hamiltonian,\n", "]\n", "large_observable_labels = [\n", " \"Re S_0d\",\n", " \"Im S_0d\",\n", " \"Re shifted H_0d\",\n", " \"Im shifted H_0d\",\n", "]\n", "\n", "# -------------------------Step 2-------------------------\n", "\n", "# Optimize the problem for quantum execution.\n", "# Select a real backend and transpile the parameterized circuit to ISA form.\n", "large_backend = service.backend(\"ibm_boston\")\n", "large_pass_manager = generate_preset_pass_manager(\n", " backend=large_backend, optimization_level=3, routing_method=\"none\"\n", ")\n", "large_isa_circuit = large_pass_manager.run(large_extended_swap_test)\n", "large_isa_observables = [\n", " observable.apply_layout(large_isa_circuit.layout)\n", " for observable in large_observables\n", "]\n", "\n", "large_two_qubit_gate_count = sum(\n", " instruction.operation.num_qubits == 2\n", " for instruction in large_isa_circuit.data\n", ")\n", "\n", "print(f\"Backend: {large_backend.name}\")\n", "print(f\"System qubits: {large_num_qubits}\")\n", "print(f\"Total circuit qubits: {large_isa_circuit.num_qubits}\")\n", "print(f\"Krylov dimension: {large_krylov_dim}\")\n", "print(f\"Time step: {large_dt:.6f}\")\n", "print(f\"Shift tau: {large_shift_tau}\")\n", "print(f\"ISA circuit depth: {large_isa_circuit.depth()}\")\n", "print(f\"ISA two-qubit gates: {large_two_qubit_gate_count}\")\n", "print(\n", " f\"ISA two-qubit depth: {large_isa_circuit.depth(lambda x: x[0].num_qubits == 2)}\"\n", ")\n", "\n", "# -------------------------Step 3-------------------------\n", "\n", "# Execute on quantum hardware with Qiskit Runtime primitives.\n", "# Submit one job per d, with all four observables in that job.\n", "large_d_values = list(range(1, large_krylov_dim))\n", "retrieve_batch_id = None\n", "large_jobs = []\n", "\n", "if retrieve_batch_id is None:\n", " large_estimator_options = {\n", " \"default_shots\": 8192,\n", " \"dynamical_decoupling\": {\n", " \"enable\": True,\n", " \"sequence_type\": \"XpXm\",\n", " },\n", " \"resilience\": {\n", " \"measure_mitigation\": True,\n", " \"measure_noise_learning\": {\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 256,\n", " },\n", " \"layer_noise_learning\": {\n", " \"max_layers_to_learn\": 4,\n", " \"layer_pair_depths\": [0, 1, 2, 4, 16, 32],\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 128,\n", " },\n", " \"zne_mitigation\": True,\n", " \"zne\": {\n", " \"amplifier\": \"pea\",\n", " \"noise_factors\": [1.0, 1.5, 2.0],\n", " \"extrapolator\": (\"exponential\", \"linear\"),\n", " },\n", " },\n", " \"twirling\": {\n", " \"enable_gates\": True,\n", " \"enable_measure\": True,\n", " \"num_randomizations\": 32,\n", " \"shots_per_randomization\": 256,\n", " \"strategy\": \"active-accum\",\n", " },\n", " }\n", "\n", " with Batch(backend=large_backend) as large_batch:\n", " large_batch_id = large_batch.session_id\n", " large_estimator = EstimatorV2(\n", " mode=large_batch,\n", " options=large_estimator_options,\n", " )\n", " # Krylov quantum diagonalization of lattice Hamiltonians -> TUT_KQDOLH.\n", " large_estimator.options.environment.job_tags = [\"TUT_KQDOLH\"]\n", "\n", " for d in large_d_values:\n", " parameter_values = [d * large_dt]\n", " pubs_for_d = [\n", " (large_isa_circuit, observable, parameter_values)\n", " for observable in large_isa_observables\n", " ]\n", " job = large_estimator.run(pubs_for_d)\n", " large_jobs.append(job)\n", " print(\n", " f\"Submitted d={d}: job_id={job.job_id()}, \"\n", " f\"PUBs={len(pubs_for_d)}\"\n", " )\n", "\n", " print(f\"Batch ID: {large_batch_id}\")\n", "else:\n", " large_batch_id = retrieve_batch_id\n", " large_jobs = service.jobs(\n", " session_id=large_batch_id,\n", " limit=None,\n", " descending=False,\n", " )\n", "\n", "large_job_ids_by_d = {\n", " d: job.job_id() for d, job in zip(large_d_values, large_jobs)\n", "}\n", "print(f\"Job IDs by d: {large_job_ids_by_d}\")\n", "\n", "# -------------------------Step 4-------------------------\n", "\n", "# Post-process the quantum results and solve the classical GEVP.\n", "# Reconstruct the first rows of S and the shifted Hamiltonian matrix.\n", "large_s_row_est = np.zeros(large_krylov_dim, dtype=complex)\n", "large_h_shifted_row_est = np.zeros(large_krylov_dim, dtype=complex)\n", "large_s_row_est[0] = 1.0\n", "large_h_shifted_row_est[0] = large_ref_energy - large_shift_tau\n", "\n", "for d, job in zip(large_d_values, large_jobs):\n", " job_result = job.result()\n", "\n", " if len(job_result) != len(large_observables):\n", " raise RuntimeError(\n", " f\"Expected {len(large_observables)} PUB results for d={d}, \"\n", " f\"but received {len(job_result)}.\"\n", " )\n", "\n", " expectation_values = [\n", " np.asarray(pub_result.data.evs).reshape(-1)[0]\n", " for pub_result in job_result\n", " ]\n", "\n", " large_s_row_est[d] = expectation_values[0] + 1j * expectation_values[1]\n", " large_h_shifted_row_est[d] = (\n", " expectation_values[2] + 1j * expectation_values[3]\n", " )\n", "\n", "# H_0d = shifted_H_0d + tau * S_0d.\n", "large_h_row_est = large_h_shifted_row_est + large_shift_tau * large_s_row_est\n", "\n", "large_s_matrix_est = la.toeplitz(\n", " large_s_row_est.conj(),\n", " large_s_row_est,\n", ")\n", "large_h_matrix_est = la.toeplitz(\n", " large_h_row_est.conj(),\n", " large_h_row_est,\n", ")\n", "\n", "large_exact_gnd = min(diagonalize_single_1_subspace(large_hamiltonian))\n", "\n", "large_energy_convergence = []\n", "\n", "with np.printoptions(precision=5, suppress=True):\n", " print(\"Estimated first row of S:\")\n", " print(large_s_row_est)\n", " print(\"Estimated first row of H:\")\n", " print(large_h_row_est)\n", " print(\"exact ground state energy: \", large_exact_gnd)\n", "\n", " for r in range(1, large_krylov_dim + 1):\n", " energy_est_kqd, coeffs_est, retained_est = solve_thresholded_gevp(\n", " large_h_matrix_est[:r, :r],\n", " large_s_matrix_est[:r, :r],\n", " threshold=5e-2,\n", " )\n", " large_energy_convergence.append(energy_est_kqd)\n", " print(\n", " f\"Krylov ground state energy (dim={r}, retained={retained_est}): \",\n", " energy_est_kqd,\n", " )\n", "\n", "plot_energy_convergence(\n", " large_energy_convergence, large_exact_gnd, large_krylov_dim\n", ")" ] }, { "cell_type": "markdown", "id": "app-md-intro", "metadata": {}, "source": [ "## Appendix: The Hamiltonian-function (spectral-filter) viewpoint\n", "\n", "The main workflow presented KQD operationally: build a Krylov basis from real-time-evolved states, estimate the projected matrices $\\mathcal{H}$ and $\\mathcal{S}$, and solve the GEVP.\n", "This appendix revisits the *same* computation from a complementary angle that explains **why** KQD works: the *Hamiltonian-function*, or *spectral-filter*, viewpoint [\\[3\\]](#references), [\\[5\\]](#references).\n", "It reuses the 12-qubit model, time step $\\Delta t$, and the Krylov solution already obtained above; no new circuit execution is required.\n", "\n", "### The reference state as an energy distribution\n", "\n", "Let the Hamiltonian have the eigendecomposition\n", "\n", "$$\n", "\\begin{equation*}\n", "H=\\sum_{m} E_m\\,|E_m\\rangle\\!\\langle E_m|,\n", "\\qquad E_0\\le E_1\\le\\cdots,\n", "\\end{equation*}\n", "$$\n", "\n", "with energy eigenstates $|E_m\\rangle$.\n", "Any reference state can be expanded in this eigenbasis,\n", "\n", "$$\n", "\\begin{equation*}\n", "|\\psi_{0}\\rangle=\\sum_m a_m\\,|E_m\\rangle,\n", "\\qquad a_m=\\langle E_m|\\psi_{0}\\rangle,\n", "\\end{equation*}\n", "$$\n", "\n", "so it carries a **spectral weight** $p_m=|a_m|^2$ at each energy $E_m$.\n", "The reference energy is the mean of this distribution, $\\langle H\\rangle_{0}=\\sum_m p_m E_m$.\n", "\n", "The eigendecomposition of a general $n$-qubit Hamiltonian is exponentially costly, so this picture is only a diagnostic, which is *not* part of the algorithm.\n", "Here, however, we can compute it cheaply for the *same* $n=12$ problem: the Heisenberg Hamiltonian conserves the total excitation number, and the reference $|000001000000\\rangle$ carries a single excitation, so its entire spectral content lives in the single-excitation subspace whose dimension grows only linearly with $n$.\n", "We therefore reuse the exact single-excitation block (already used above as a benchmark) and read off the reference distribution $\\{(E_m, p_m)\\}$ inside that subspace.\n", "\n" ] }, { "cell_type": "code", "execution_count": 20, "id": "app-code-decomp", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Single-excitation subspace dimension: 12\n", "Subspace ground-state energy: 3.136297\n", "Reference energy (sum p_m E_m): 7.000000\n", "Reference weight on subspace ground: 0.163827\n" ] } ], "source": [ "# Exact single-excitation-subspace decomposition of the reference state.\n", "# This reuses make_heisenberg_hamiltonian / _basis_state_transition_amplitude_sparse\n", "# and the n=12 `hamiltonian` and `ref_bitstring` defined in the small-scale example.\n", "single_excitation_states = [1 << k for k in range(num_qubits)]\n", "\n", "h_single = np.array(\n", " [\n", " [\n", " _basis_state_transition_amplitude_sparse(hamiltonian, bra, ket)\n", " for ket in single_excitation_states\n", " ]\n", " for bra in single_excitation_states\n", " ]\n", ")\n", "h_single = 0.5 * (h_single + h_single.conj().T)\n", "\n", "subspace_evals, subspace_evecs = la.eigh(h_single)\n", "subspace_evals = np.real(subspace_evals)\n", "\n", "# Reference-state coordinates inside the single-excitation subspace.\n", "ref_position = single_excitation_states.index(int(ref_bitstring, 2))\n", "ref_in_subspace = np.zeros(num_qubits, dtype=complex)\n", "ref_in_subspace[ref_position] = 1.0\n", "\n", "# Amplitudes and spectral weights of the reference in the energy eigenbasis.\n", "ref_eigen_amplitudes = subspace_evecs.conj().T @ ref_in_subspace\n", "ref_spectral_weights = np.abs(ref_eigen_amplitudes) ** 2\n", "\n", "print(f\"Single-excitation subspace dimension: {num_qubits}\")\n", "print(f\"Subspace ground-state energy: {subspace_evals[0]:.6f}\")\n", "print(\n", " f\"Reference energy (sum p_m E_m): {np.sum(ref_spectral_weights * subspace_evals):.6f}\"\n", ")\n", "print(f\"Reference weight on subspace ground: {ref_spectral_weights[0]:.6f}\")" ] }, { "cell_type": "markdown", "id": "app-md-filter", "metadata": {}, "source": [ "### KQD learns a filter that reshapes this distribution\n", "\n", "A **Hamiltonian function** $f(H)$ is defined through spectral calculus,\n", "\n", "$$\n", "\\begin{equation*}\n", "f(H)=\\sum_m f(E_m)\\,|E_m\\rangle\\!\\langle E_m|,\n", "\\end{equation*}\n", "$$\n", "\n", "or in other words, a weighted sum of eigenprojectors.\n", "Applying it to the reference reshapes each spectral amplitude, $a_m \\to a_m f(E_m)$:\n", "\n", "$$\n", "\\begin{equation*}\n", "f(H)\\,|\\psi_{0}\\rangle=\\sum_m a_m\\,f(E_m)\\,|E_m\\rangle .\n", "\\end{equation*}\n", "$$\n", "\n", "If $f$ were sharply peaked at the lowest energy ($f(E_0)=1$ and $f(E_m)\\approx 0$ otherwise) then $f(H)$ would act as a ground-state projector, and the normalized output would be (close to) the ground state.\n", "Therefore, a good low-pass **spectral filter** in energy is exactly what we want.\n", "\n", "KQD does not prescribe $f$ in advance. Instead it expands the filter in the real-time-evolution basis,\n", "\n", "$$\n", "\\begin{equation*}\n", "f_{\\rm KQD}(E)=\\sum_{\\ell=0}^{r-1} c_\\ell\\, e^{-i\\ell\\Delta t E},\n", "\\end{equation*}\n", "$$\n", "\n", "a *trigonometric* function of the energy whose coefficients $\\{c_\\ell\\}$ are precisely the GEVP eigenvector solved for above.\n", "Minimizing the Rayleigh quotient $\\mathbf{c}^\\dagger\\mathcal{H}\\mathbf{c}/\\mathbf{c}^\\dagger\\mathcal{S}\\mathbf{c}$ is therefore the same as *learning* the filter that best suppresses the excited-state weight of the reference.\n", "Larger Krylov dimension $r$ gives the filter more degrees of freedom and a sharper peak at the ground-state energy.\n", "\n", "The helper below evaluates this learned filter on an energy axis; we then apply it to the reference distribution obtained above.\n", "\n" ] }, { "cell_type": "code", "execution_count": 21, "id": "app-code-filter", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Ground-state overlap (reference): 0.1638\n", "Ground-state overlap (filtered): 0.9638\n", "Mean energy (reference): 7.000000\n", "Mean energy (filtered): 3.164503\n" ] } ], "source": [ "def trigonometric_krylov_filter(\n", " coeffs: np.ndarray,\n", " energies: np.ndarray,\n", " time_step: float,\n", ") -> np.ndarray:\n", " \"\"\"Evaluate the learned Krylov filter f(E) = sum_l c_l exp(-i l dt E).\"\"\"\n", " values = np.zeros_like(energies, dtype=complex)\n", " for ell, coeff in enumerate(coeffs):\n", " values += coeff * np.exp(-1j * ell * time_step * energies)\n", " return values\n", "\n", "\n", "def filtered_spectral_weights(\n", " weights: np.ndarray,\n", " filter_values: np.ndarray,\n", ") -> np.ndarray:\n", " \"\"\"Reshape spectral weights by |f(E)|^2 and renormalize.\"\"\"\n", " reshaped = weights * np.abs(filter_values) ** 2\n", " return reshaped / np.sum(reshaped)\n", "\n", "\n", "# Recover the KQD coefficients from the already-estimated projected matrices.\n", "# The shift only moves H by tau * S, so it does not change the GEVP eigenvector;\n", "# we solve at the full Krylov dimension used in the small-scale example.\n", "_, kqd_coeffs, _ = solve_thresholded_gevp(\n", " h_matrix_est,\n", " s_matrix_est,\n", " threshold=1e-12,\n", ")\n", "\n", "filter_on_spectrum = trigonometric_krylov_filter(\n", " kqd_coeffs, subspace_evals, dt\n", ")\n", "filtered_weights = filtered_spectral_weights(\n", " ref_spectral_weights, filter_on_spectrum\n", ")\n", "\n", "print(f\"Ground-state overlap (reference): {ref_spectral_weights[0]:.4f}\")\n", "print(f\"Ground-state overlap (filtered): {filtered_weights[0]:.4f}\")\n", "print(\n", " f\"Mean energy (reference): {np.sum(ref_spectral_weights * subspace_evals):.6f}\"\n", ")\n", "print(\n", " f\"Mean energy (filtered): {np.sum(filtered_weights * subspace_evals):.6f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "app-md-plot", "metadata": {}, "source": [ "### Visualize the filter and its flexibility\n", "\n", "We first show the learned filter at the full Krylov dimension used above, then track how it sharpens as the dimension $r$ grows.\n", "\n", "The bars show the reference spectral weights $p_m$ (before) and the filtered weights $p_m|f_{\\rm KQD}(E_m)|^2$ (after), together with the learned filter intensity $|f_{\\rm KQD}(E)|^2$ on a continuous energy axis.\n", "The filter concentrates the weight onto the lowest energy of the single-excitation subspace — the same energy that the KQD estimate converged to in the small-scale example.\n", "Note that this is the ground state *within the single-excitation sector*, which is the relevant target for this excitation-conserving reference state, not the global ground state.\n", "\n" ] }, { "cell_type": "code", "execution_count": 22, "id": "app-code-plot", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(8, 4))\n", "\n", "visible = (ref_spectral_weights > 1e-4) | (filtered_weights > 1e-4)\n", "ax.bar(\n", " subspace_evals[visible],\n", " ref_spectral_weights[visible],\n", " width=0.18,\n", " alpha=0.45,\n", " label=\"reference $p_m$\",\n", ")\n", "ax.bar(\n", " subspace_evals[visible],\n", " filtered_weights[visible],\n", " width=0.14,\n", " alpha=0.9,\n", " label=r\"filtered $p_m\\,|f_{\\rm KQD}(E_m)|^2$\",\n", ")\n", "\n", "energy_grid = np.linspace(\n", " subspace_evals.min() - 0.5, subspace_evals.max() + 0.5, 800\n", ")\n", "filter_intensity = (\n", " np.abs(trigonometric_krylov_filter(kqd_coeffs, energy_grid, dt)) ** 2\n", ")\n", "filter_intensity /= filter_intensity.max()\n", "ax.plot(\n", " energy_grid,\n", " filter_intensity,\n", " color=\"k\",\n", " linewidth=2,\n", " label=r\"$|f_{\\rm KQD}(E)|^2$ (normalized)\",\n", ")\n", "\n", "ax.axvline(\n", " subspace_evals[0],\n", " color=\"C3\",\n", " linestyle=\"--\",\n", " linewidth=1,\n", " label=\"subspace ground energy\",\n", ")\n", "ax.set_xlabel(\"Energy eigenvalue $E_m$\")\n", "ax.set_ylabel(\"Spectral weight\")\n", "ax.set_ylim(0, 1)\n", "ax.legend(loc=\"upper right\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "app-md-sweep", "metadata": {}, "source": [ "### Increasing the Krylov dimension: flexibility of the learned function\n", "\n", "Recall that the learned filter is a trigonometric polynomial in the energy with $r$ coefficients,\n", "\n", "$$\n", "\\begin{equation*}\n", "f_{\\rm KQD}(E)=\\sum_{\\ell=0}^{r-1} c_\\ell\\, e^{-i\\ell\\Delta t E}.\n", "\\end{equation*}\n", "$$\n", "\n", "The Krylov dimension $r$ is exactly the number of free coefficients, so it controls the **flexibility** of the function.\n", "A small $r$ can only produce a broad, gently varying filter that leaks weight into low-lying excited states; when $r$ increases, the filter can form a narrower peak at the target energy and suppress the remaining excited-state weight more aggressively.\n", "This is the spectral-filter counterpart of the energy convergence observed in the small-scale example: as $r$ grows, the filtered distribution collapses onto the subspace ground state and the estimated energy decreases toward it.\n", "\n", "We reuse the projected matrices already estimated above and simply solve the GEVP at each leading $r\\times r$ block, then evaluate and plot the corresponding filter.\n", "\n" ] }, { "cell_type": "code", "execution_count": 23, "id": "app-code-sweep", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " r retained ground overlap filtered energy\n", "-- -------- -------------- ---------------\n", " 2 2 0.4549 4.184510\n", " 4 4 0.8373 3.310168\n", " 6 6 0.9706 3.158100\n", " 8 7 0.9701 3.158725\n", "10 8 0.9638 3.164503\n", "\n", "Subspace ground-state energy (target): 3.136297\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Sweep the Krylov dimension using the leading r x r blocks of the estimated matrices.\n", "sweep_dims = [r for r in (2, 4, 6, 8, krylov_dim) if r <= krylov_dim]\n", "sweep_dims = sorted(set(sweep_dims))\n", "\n", "energy_grid = np.linspace(\n", " subspace_evals.min() - 0.5, subspace_evals.max() + 0.5, 800\n", ")\n", "\n", "sweep_cases = []\n", "print(\" r retained ground overlap filtered energy\")\n", "print(\"-- -------- -------------- ---------------\")\n", "for r in sweep_dims:\n", " _, coeffs_r, retained_r = solve_thresholded_gevp(\n", " h_matrix_est[:r, :r],\n", " s_matrix_est[:r, :r],\n", " threshold=1e-12,\n", " )\n", " filter_on_spectrum_r = trigonometric_krylov_filter(\n", " coeffs_r, subspace_evals, dt\n", " )\n", " filtered_weights_r = filtered_spectral_weights(\n", " ref_spectral_weights, filter_on_spectrum_r\n", " )\n", " filtered_energy_r = float(np.sum(filtered_weights_r * subspace_evals))\n", "\n", " sweep_cases.append((r, coeffs_r, filtered_weights_r))\n", " print(\n", " f\"{r:2d} {retained_r:8d} {filtered_weights_r[0]:14.4f} {filtered_energy_r:15.6f}\"\n", " )\n", "\n", "print(f\"\\nSubspace ground-state energy (target): {subspace_evals[0]:.6f}\")\n", "\n", "# One panel per Krylov dimension: filtered spectrum (bars) + filter intensity (curve).\n", "fig, axes = plt.subplots(\n", " len(sweep_cases),\n", " 1,\n", " figsize=(8, 2.1 * len(sweep_cases)),\n", " sharex=True,\n", ")\n", "axes = np.atleast_1d(axes)\n", "\n", "for idx, (ax, (r, coeffs_r, filtered_weights_r)) in enumerate(\n", " zip(axes, sweep_cases)\n", "):\n", " visible = (ref_spectral_weights > 1e-4) | (filtered_weights_r > 1e-4)\n", " ax.bar(\n", " subspace_evals[visible],\n", " ref_spectral_weights[visible],\n", " width=0.18,\n", " alpha=0.35,\n", " color=\"C0\",\n", " label=\"reference $p_m$\" if idx == 0 else None,\n", " )\n", " ax.bar(\n", " subspace_evals[visible],\n", " filtered_weights_r[visible],\n", " width=0.14,\n", " alpha=0.9,\n", " color=\"C1\",\n", " label=\"filtered weights\" if idx == 0 else None,\n", " )\n", "\n", " filter_intensity_r = (\n", " np.abs(trigonometric_krylov_filter(coeffs_r, energy_grid, dt)) ** 2\n", " )\n", " filter_intensity_r /= filter_intensity_r.max()\n", " ax.plot(energy_grid, filter_intensity_r, color=\"k\", linewidth=2)\n", "\n", " ax.axvline(subspace_evals[0], color=\"C3\", linestyle=\"--\", linewidth=1)\n", " ax.set_ylim(0, 1)\n", " ax.set_ylabel(\"weight\")\n", " ax.legend(loc=\"upper right\", title=f\"$r={r}$\")\n", "\n", "axes[-1].set_xlabel(\"Energy eigenvalue $E_m$\")\n", "fig.suptitle(\n", " r\"KQD-learned filter $|f_{\\rm KQD}(E)|^2$ sharpening with Krylov dimension $r$\",\n", " y=1.0,\n", ")\n", "fig.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "4a6b7ea8", "metadata": {}, "source": [ "## Next steps\n", "\n", "If you found this work interesting, you might be interested in the following material:\n", "\n", "\n", " * [Sample-based Krylov quantum diagonalization of a fermionic lattice model](/docs/tutorials/sample-based-krylov-quantum-diagonalization)\n", " * [Ground-state energy estimation of the Heisenberg chain with VQE](/docs/tutorials/spin-chain-vqe)\n", " * [Quantum diagonalization algorithms course](/learning/courses/quantum-diagonalization-algorithms/krylov)\n", "\n", "\n", "## References\n", "\n", "\\[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).\n", "\n", "\\[2] N. Yoshioka, M. Amico, W. Kirby, et al., *Diagonalization of large many-body Hamiltonians on a quantum processor*, [arXiv:2407.14431](https://arxiv.org/abs/2407.14431) (2024).\n", "\n", "\\[3] R. M. Parrish and P. L. McMahon, *Quantum filter diagonalization: quantum eigendecomposition without full quantum phase estimation*, *Physical Review Letters* **122**, 230401 (2019).\n", "\n", "\\[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).\n", "\n", "\\[5] G. Lee, M. Kang, J. Hong, S. Fomichev and J. Huh, *Filtered Quantum Phase Estimation*, [arXiv:2510.04294](https://arxiv.org/abs/2510.04294) (2025).\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "hours": 2, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" }, "qpuSeconds": 4200 }, "nbformat": 4, "nbformat_minor": 5 }