{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "ee388de5-2507-48ff-97a6-6777698d6256",
      "metadata": {},
      "source": [
        "---\n",
        "title: Sample-based Krylov quantum diagonalization of a fermionic lattice model\n",
        "description: Use the sample-based quantum diagonalization algorithm to simulate the single-impurity Anderson model using noisy quantum hardware.\n",
        "---\n",
        "\n",
        "{/* cspell:ignore fontdict fontsize nocc SQKD DMRG textrm varepsilon vecs pqrs ijkl */}\n",
        "\n",
        "# Sample-based Krylov quantum diagonalization of a fermionic lattice model\n",
        "\n",
        "*Usage estimate: Nine seconds on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a1b2c3d4-e5f6-7890-abcd-ef1234567890",
      "metadata": {},
      "source": [
        "## Learning outcomes\n",
        "\n",
        "After going through this tutorial, users should understand:\n",
        "\n",
        "* How to use the [SQD Qiskit addon](https://github.com/Qiskit/qiskit-addon-sqd) to approximate the ground state energy of a lattice model using bitstrings sampled from a quantum processing unit (QPU).\n",
        "* How to use [ffsim](https://github.com/qiskit-community/ffsim) to construct time evolution circuits for fermionic simulation.\n",
        "* How to combine samples from multiple circuits for post-processing with the sample-based Krylov diagonalization (SKQD) algorithm.\n",
        "\n",
        "## Prerequisites\n",
        "\n",
        "We suggest that users are familiar with the following topics before going through this tutorial:\n",
        "\n",
        "* [Sample-based quantum diagonalization of a chemistry Hamiltonian](/docs/tutorials/sample-based-quantum-diagonalization)\n",
        "* [Krylov quantum diagonalization of lattice Hamiltonians](/docs/tutorials/krylov-quantum-diagonalization)\n",
        "* [Qiskit primitives](/docs/guides/primitives)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "dc5cc74e-06bf-45ac-a69b-81778138e08f",
      "metadata": {},
      "source": [
        "## Background\n",
        "\n",
        "This tutorial shows how to use sample-based quantum diagonalization (SQD) to estimate the ground state energy of a fermionic lattice model. Specifically, we study the one-dimensional single-impurity Anderson model (SIAM), which is used to describe magnetic impurities embedded in metals.\n",
        "\n",
        "This tutorial follows a similar workflow to the related tutorial [Sample-based quantum diagonalization of a chemistry Hamiltonian](/docs/tutorials/sample-based-quantum-diagonalization). However, a key difference lies in how the quantum circuits are built. The other tutorial uses a heuristic variational ansatz, which is appealing for chemistry Hamiltonians with potentially millions of interaction terms. On the other hand, this tutorial uses circuits that approximate time evolution by the Hamiltonian. Such circuits can be deep, which makes this approach better for applications to lattice models. The state vectors prepared by these circuits form the basis for a [Krylov subspace](https://en.wikipedia.org/wiki/Krylov_subspace), and as a result, the algorithm provably and efficiently converges to the ground state, under suitable assumptions.\n",
        "\n",
        "The approach used in this tutorial can be viewed as a combination of the techniques used in SQD and [Krylov quantum diagonalization (KQD)](https://arxiv.org/abs/2407.14431). The combined approach is sometimes referred to as sample-based Krylov quantum diagonalization (SQKD). See [Krylov quantum diagonalization of lattice Hamiltonians](/docs/tutorials/krylov-quantum-diagonalization) for a tutorial on the KQD method.\n",
        "\n",
        "This tutorial is based on the work [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702), which can be referred to for more details.\n",
        "\n",
        "### Single-impurity Anderson model (SIAM)\n",
        "\n",
        "The one-dimensional SIAM Hamiltonian is a sum of three terms:\n",
        "\n",
        "$$\n",
        "H = H_{\\textrm{imp}}+ H_\\textrm{bath} + H_\\textrm{hyb},\n",
        "$$\n",
        "\n",
        "where\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "  H_\\textrm{imp} &= \\varepsilon \\left( \\hat{n}_{d\\uparrow} + \\hat{n}_{d\\downarrow} \\right) + U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}, \\\\\n",
        "  H_\\textrm{bath} &= -t \\sum_{\\substack{\\mathbf{j} = 0\\\\ \\sigma\\in \\{\\uparrow, \\downarrow\\}}}^{L-1} \\left(\\hat{c}^\\dagger_{\\mathbf{j}, \\sigma}\\hat{c}_{\\mathbf{j}+1, \\sigma} + \\hat{c}^\\dagger_{\\mathbf{j}+1, \\sigma}\\hat{c}_{\\mathbf{j}, \\sigma} \\right), \\\\\n",
        "  H_\\textrm{hyb} &= V\\sum_{\\sigma \\in \\{\\uparrow, \\downarrow \\}} \\left(\\hat{d}^\\dagger_\\sigma \\hat{c}_{0, \\sigma} + \\hat{c}^\\dagger_{0, \\sigma} \\hat{d}_{\\sigma} \\right).\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "Here, $c^\\dagger_{\\mathbf{j},\\sigma}/c_{\\mathbf{j},\\sigma}$ are the fermionic creation/annihilation operators for the $\\mathbf{j}^{\\textrm{th}}$ bath site with spin $\\sigma$, $\\hat{d}^\\dagger_{\\sigma}/\\hat{d}_{\\sigma}$ are creation/annihilation operators for the impurity mode, and $\\hat{n}_{d\\sigma} = \\hat{d}^\\dagger_{\\sigma} \\hat{d}_{\\sigma}$. $t$, $U$, and $V$ are real numbers describing the hopping, on-site, hybridization interactions, and $\\varepsilon$ is a real number specifying the chemical potential.\n",
        "\n",
        "Note that the Hamiltonian is a specific instance of the generic interaction-electron Hamiltonian,\n",
        "\n",
        "$$\n",
        "\\begin{align*}\n",
        "  H &= \\sum_{\\substack{p, q \\\\ \\sigma}} h_{pq} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}_{q\\sigma}  +  \\sum_{\\substack{p, q, r, s \\\\ \\sigma \\tau}} \\frac{h_{pqrs}}{2} \\hat{a}^\\dagger_{p\\sigma} \\hat{a}^\\dagger_{q\\tau} \\hat{a}_{s\\tau} \\hat{a}_{r\\sigma} \\\\\n",
        "  &= H_1 + H_2,\n",
        "\\end{align*}\n",
        "$$\n",
        "\n",
        "where $H_1$ consists of one-body terms, which are quadratic in the fermionic creation and annihilation operators, and $H_2$ consists of two-body terms, which are quartic. For the SIAM,\n",
        "\n",
        "$$\n",
        "H_2 = U \\hat{n}_{d\\uparrow}\\hat{n}_{d\\downarrow}\n",
        "$$\n",
        "\n",
        "and $H_1$ contains the rest of the terms in the Hamiltonian. In order to represent the Hamiltonian programmatically, we store the matrix $h_{pq}$ and the tensor $h_{pqrs}$.\n",
        "\n",
        "### Position and momentum bases\n",
        "\n",
        "Due to the approximate translational symmetry in $H_\\textrm{bath}$, we don't expect the ground state to be sparse in the position basis (the orbital basis in which the Hamiltonian is specified above). The performance of SQD is guaranteed only if the ground state is sparse, that is, it has significant weight on only a small number of computational basis states. To improve the sparsity of the ground state, we perform the simulation in the orbital basis in which $H_\\textrm{bath}$ is diagonal. We call this basis the *momentum basis*. Because $H_\\textrm{bath}$ is a quadratic fermionic Hamiltonian, it can be efficiently diagonalized by an orbital rotation.\n",
        "\n",
        "### Approximate time evolution by the Hamiltonian\n",
        "\n",
        "To approximate time evolution by the Hamiltonian, we use a second order Trotter-Suzuki decomposition,\n",
        "\n",
        "$$\n",
        "  e^{-i \\Delta t H} \\approx e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2}.\n",
        "$$\n",
        "\n",
        "Under the [Jordan-Wigner transformation](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation), time evolution by $H_2$ amounts to a single [CPhase](/docs/api/qiskit/qiskit.circuit.library.CPhaseGate) gate between the spin-up and spin-down orbitals at the impurity site. Because $H_1$ is a quadratic fermionic Hamiltonian, time evolution by $H_1$ amounts to an orbital rotation.\n",
        "\n",
        "The Krylov basis states $\\{ |\\psi_k\\rangle \\}_{k=0}^{D-1}$, where $D$ is the dimension of the Krylov subspace, are formed by repeated application of a single Trotter step, so\n",
        "\n",
        "$$\n",
        "  |\\psi_k\\rangle \\approx \\left[e^{-i\\frac{\\Delta t}{2} H_2} e^{-i\\Delta t H_1} e^{-i\\frac{\\Delta t}{2} H_2} \\right]^k\\ket{\\psi_0}.\n",
        "$$\n",
        "\n",
        "In the following SQD-based workflow, we will sample from this set of circuits and post-process the combined set of bitstrings with SQD. This approach contrasts with the one used in the related tutorial [Sample-based quantum diagonalization of a chemistry Hamiltonian](/docs/tutorials/sample-based-quantum-diagonalization), where samples were drawn from a single heuristic variational circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "df66d697-102c-4a2c-80f3-f67fdda05573",
      "metadata": {},
      "source": [
        "## Requirements\n",
        "\n",
        "Before starting this tutorial, ensure that you have the following installed:\n",
        "\n",
        "* Qiskit SDK v1.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
        "* Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n",
        "* SQD Qiskit addon v0.11 or later (`pip install qiskit-addon-sqd`)\n",
        "* ffsim v0.0.72 or later (`pip install ffsim`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c3d4e5f6-a7b8-9012-cdef-123456789012",
      "metadata": {},
      "source": [
        "## Small-scale simulator example\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8540487a-8033-49c2-9f30-022336105f64",
      "metadata": {},
      "source": [
        "### Step 1: Map problem to a quantum circuit\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4e6652bd-c97d-4f2f-98a7-d44857805cf1",
      "metadata": {},
      "source": [
        "First, we generate the SIAM Hamiltonian in the position basis. The Hamiltonian is represented by the matrix $h_{pq}$ and the tensor $h_{pqrs}$. Then, we rotate it to the momentum basis. In the position basis, we place the impurity at the first site. However, when we rotate to the momentum basis, we move the impurity to a central site to facilitate interactions with other orbitals.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "b5cb9c28-4721-4141-8665-96885038e210",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pyscf.fci\n",
        "\n",
        "\n",
        "def siam_hamiltonian(\n",
        "    norb: int,\n",
        "    hopping: float,\n",
        "    onsite: float,\n",
        "    hybridization: float,\n",
        "    chemical_potential: float,\n",
        ") -> tuple[np.ndarray, np.ndarray]:\n",
        "    \"\"\"Hamiltonian for the single-impurity Anderson model.\"\"\"\n",
        "    # Place the impurity on the first site\n",
        "    impurity_orb = 0\n",
        "\n",
        "    # One body matrix elements in the \"position\" basis\n",
        "    h1e = np.zeros((norb, norb))\n",
        "    np.fill_diagonal(h1e[:, 1:], -hopping)\n",
        "    np.fill_diagonal(h1e[1:, :], -hopping)\n",
        "    h1e[impurity_orb, impurity_orb + 1] = -hybridization\n",
        "    h1e[impurity_orb + 1, impurity_orb] = -hybridization\n",
        "    h1e[impurity_orb, impurity_orb] = chemical_potential\n",
        "\n",
        "    # Two body matrix elements in the \"position\" basis\n",
        "    h2e = np.zeros((norb, norb, norb, norb))\n",
        "    h2e[impurity_orb, impurity_orb, impurity_orb, impurity_orb] = onsite\n",
        "\n",
        "    return h1e, h2e\n",
        "\n",
        "\n",
        "def momentum_basis(norb: int) -> np.ndarray:\n",
        "    \"\"\"Get the orbital rotation to change from the position to the momentum basis.\"\"\"\n",
        "    n_bath = norb - 1\n",
        "\n",
        "    # Orbital rotation that diagonalizes the bath (non-interacting system)\n",
        "    hopping_matrix = np.zeros((n_bath, n_bath))\n",
        "    np.fill_diagonal(hopping_matrix[:, 1:], -1)\n",
        "    np.fill_diagonal(hopping_matrix[1:, :], -1)\n",
        "    _, vecs = np.linalg.eigh(hopping_matrix)\n",
        "\n",
        "    # Expand to include impurity\n",
        "    orbital_rotation = np.zeros((norb, norb))\n",
        "    # Impurity is on the first site\n",
        "    orbital_rotation[0, 0] = 1\n",
        "    orbital_rotation[1:, 1:] = vecs\n",
        "\n",
        "    # Move the impurity to the center\n",
        "    new_index = n_bath // 2\n",
        "    perm = np.r_[1 : (new_index + 1), 0, (new_index + 1) : norb]\n",
        "    orbital_rotation = orbital_rotation[:, perm]\n",
        "\n",
        "    return orbital_rotation\n",
        "\n",
        "\n",
        "def rotated(\n",
        "    h1e: np.ndarray, h2e: np.ndarray, orbital_rotation: np.ndarray\n",
        ") -> tuple[np.ndarray, np.ndarray]:\n",
        "    \"\"\"Rotate the orbital basis of a Hamiltonian.\"\"\"\n",
        "    h1e_rotated = np.einsum(\n",
        "        \"ab,Aa,Bb->AB\",\n",
        "        h1e,\n",
        "        orbital_rotation,\n",
        "        orbital_rotation.conj(),\n",
        "        optimize=\"greedy\",\n",
        "    )\n",
        "    h2e_rotated = np.einsum(\n",
        "        \"abcd,Aa,Bb,Cc,Dd->ABCD\",\n",
        "        h2e,\n",
        "        orbital_rotation,\n",
        "        orbital_rotation.conj(),\n",
        "        orbital_rotation,\n",
        "        orbital_rotation.conj(),\n",
        "        optimize=\"greedy\",\n",
        "    )\n",
        "    return h1e_rotated, h2e_rotated\n",
        "\n",
        "\n",
        "# Total number of spatial orbitals, including the bath sites and the impurity\n",
        "# This should be an even number\n",
        "norb = 8\n",
        "\n",
        "# System is half-filled\n",
        "nelec = (norb // 2, norb // 2)\n",
        "# One orbital is the impurity, the rest are bath sites\n",
        "n_bath = norb - 1\n",
        "\n",
        "# Hamiltonian parameters\n",
        "hybridization = 1.0\n",
        "hopping = 1.0\n",
        "onsite = 10.0\n",
        "chemical_potential = -0.5 * onsite\n",
        "\n",
        "# Generate Hamiltonian in position basis\n",
        "h1e, h2e = siam_hamiltonian(\n",
        "    norb=norb,\n",
        "    hopping=hopping,\n",
        "    onsite=onsite,\n",
        "    hybridization=hybridization,\n",
        "    chemical_potential=chemical_potential,\n",
        ")\n",
        "\n",
        "# Rotate to momentum basis\n",
        "orbital_rotation = momentum_basis(norb)\n",
        "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n",
        "# In the momentum basis, the impurity is placed in the center\n",
        "impurity_index = n_bath // 2\n",
        "\n",
        "# Use PySCF to compute the exact ground state energy\n",
        "reference_energy, _ = pyscf.fci.direct_spin1.kernel(h1e, h2e, norb, nelec)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e888edf2-7865-41a8-be57-6ddb72dd0cc7",
      "metadata": {},
      "source": [
        "Next, we generate the circuits to produce the Krylov basis states.\n",
        "For each spin species, the initial state $\\ket{\\psi_0}$ is given by the superposition of all possible excitations of the three electrons closest to the Fermi level into the 4 closest empty modes starting from the state $|00\\cdots 0011 \\cdots 11\\rangle$, and realized by the application of seven [XXPlusYYGate](/docs/api/qiskit/qiskit.circuit.library.XXPlusYYGate)s.\n",
        "The time-evolved states are produced by successive applications of a second-order Trotter step.\n",
        "\n",
        "For a more detailed description of this model and how the circuits are designed, refer to [\"Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization\"](https://arxiv.org/abs/2501.09702).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "0f729f86-1814-4d5a-ae65-3f9614e103b3",
      "metadata": {},
      "outputs": [],
      "source": [
        "from typing import Sequence\n",
        "\n",
        "import ffsim\n",
        "import scipy\n",
        "from qiskit import QuantumCircuit, QuantumRegister\n",
        "from qiskit.circuit import CircuitInstruction, Qubit\n",
        "from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate\n",
        "\n",
        "\n",
        "def prepare_initial_state(qubits: Sequence[Qubit], norb: int, nocc: int):\n",
        "    \"\"\"Prepare initial state.\"\"\"\n",
        "    assert norb >= 8\n",
        "    x_gate = XGate()\n",
        "    rot = XXPlusYYGate(0.5 * np.pi, -0.5 * np.pi)\n",
        "    for i in range(nocc):\n",
        "        yield CircuitInstruction(x_gate, [qubits[i]])\n",
        "        yield CircuitInstruction(x_gate, [qubits[norb + i]])\n",
        "    for i in range(3):\n",
        "        for j in range(nocc - i - 1, nocc + i, 2):\n",
        "            yield CircuitInstruction(rot, [qubits[j], qubits[j + 1]])\n",
        "            yield CircuitInstruction(\n",
        "                rot, [qubits[norb + j], qubits[norb + j + 1]]\n",
        "            )\n",
        "    yield CircuitInstruction(rot, [qubits[j + 1], qubits[j + 2]])\n",
        "    yield CircuitInstruction(\n",
        "        rot, [qubits[norb + j + 1], qubits[norb + j + 2]]\n",
        "    )\n",
        "\n",
        "\n",
        "def trotter_step(\n",
        "    qubits: Sequence[Qubit],\n",
        "    time_step: float,\n",
        "    one_body_evolution: np.ndarray,\n",
        "    h2e: np.ndarray,\n",
        "    impurity_index: int,\n",
        "    norb: int,\n",
        "):\n",
        "    \"\"\"A Trotter step.\"\"\"\n",
        "    # Assume the two-body interaction is just the on-site interaction of the impurity\n",
        "    onsite = h2e[\n",
        "        impurity_index, impurity_index, impurity_index, impurity_index\n",
        "    ]\n",
        "    # Two-body evolution for half the time\n",
        "    yield CircuitInstruction(\n",
        "        CPhaseGate(-0.5 * time_step * onsite),\n",
        "        [qubits[impurity_index], qubits[norb + impurity_index]],\n",
        "    )\n",
        "    # One-body evolution for the full time\n",
        "    yield CircuitInstruction(\n",
        "        ffsim.qiskit.OrbitalRotationJW(norb, one_body_evolution), qubits\n",
        "    )\n",
        "    # Two-body evolution for half the time\n",
        "    yield CircuitInstruction(\n",
        "        CPhaseGate(-0.5 * time_step * onsite),\n",
        "        [qubits[impurity_index], qubits[norb + impurity_index]],\n",
        "    )\n",
        "\n",
        "\n",
        "# Time step\n",
        "time_step = 0.2\n",
        "# Number of Krylov basis states\n",
        "krylov_dim = 8\n",
        "\n",
        "# Initialize circuit\n",
        "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
        "circuit = QuantumCircuit(qubits)\n",
        "\n",
        "# Generate initial state\n",
        "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n",
        "    circuit.append(instruction)\n",
        "circuit.measure_all()\n",
        "\n",
        "# Create list of circuits, starting with the initial state circuit\n",
        "circuits = [circuit.copy()]\n",
        "\n",
        "# Add time evolution circuits to the list\n",
        "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n",
        "for i in range(krylov_dim - 1):\n",
        "    # Remove measurements\n",
        "    circuit.remove_final_measurements()\n",
        "    # Append another Trotter step\n",
        "    for instruction in trotter_step(\n",
        "        qubits,\n",
        "        time_step,\n",
        "        one_body_evolution,\n",
        "        h2e_momentum,\n",
        "        impurity_index,\n",
        "        norb,\n",
        "    ):\n",
        "        circuit.append(instruction)\n",
        "    # Measure qubits\n",
        "    circuit.measure_all()\n",
        "    # Add a copy of the circuit to the list\n",
        "    circuits.append(circuit.copy())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "9f2cc4d4-ecac-457a-bcae-558319668e1f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-krylov-quantum-diagonalization/extracted-outputs/9f2cc4d4-ecac-457a-bcae-558319668e1f-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circuits[0].draw(\"mpl\", scale=0.4, fold=-1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "827976ec-4815-4707-80b1-e13fb2fef309",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-krylov-quantum-diagonalization/extracted-outputs/827976ec-4815-4707-80b1-e13fb2fef309-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circuits[-1].draw(\"mpl\", scale=0.4, fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b8ca6be4-61d9-47be-8099-8712c7ecc774",
      "metadata": {},
      "source": [
        "### Step 2: Optimize problem for quantum execution\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a3304e1b-9c7c-4212-8744-d1c62292eced",
      "metadata": {},
      "source": [
        "Next, we optimize the circuit for a target hardware. For now, we'll create a generic backend with a specified number of qubits and a gate set that the time evolution circuits naturally decompose to.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "2d2fdbff-1e22-45af-a2eb-c334e4328c59",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.providers.fake_provider import GenericBackendV2\n",
        "\n",
        "backend = GenericBackendV2(\n",
        "    2 * norb, basis_gates=[\"cp\", \"xx_plus_yy\", \"p\", \"x\"]\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ce3e6a52-7b99-49b6-8294-b005efa59cfc",
      "metadata": {},
      "source": [
        "Now, we use Qiskit to transpile the circuits to the target backend.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "c8643533-9fec-40bf-a307-da8839b1e444",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.transpiler import generate_preset_pass_manager\n",
        "\n",
        "pass_manager = generate_preset_pass_manager(\n",
        "    optimization_level=3, backend=backend\n",
        ")\n",
        "isa_circuits = pass_manager.run(circuits)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6cfd3eea-e2d9-40a5-a449-1d3d790a5f2d",
      "metadata": {},
      "source": [
        "### Step 3: Execute using Qiskit primitives\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ad48e3e6-1013-45a1-942c-63766fed5819",
      "metadata": {},
      "source": [
        "After optimizing the circuits for hardware execution, we are ready to run them on the target hardware and collect samples for ground state energy estimation. After using the Sampler primitive to sample bitstrings from each circuit, we combine all of the results into a single counts dictionary and plot the top 20 most commonly sampled bitstrings.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "80eee553-60d6-4258-88ab-d8d120418c36",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.visualization import plot_histogram\n",
        "from qiskit.primitives import StatevectorSampler\n",
        "\n",
        "# Sample from the circuits\n",
        "sampler = StatevectorSampler()\n",
        "job = sampler.run(isa_circuits, shots=500)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "10af4663-7375-4b50-bae6-9f3d5106457b",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-krylov-quantum-diagonalization/extracted-outputs/10af4663-7375-4b50-bae6-9f3d5106457b-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "from qiskit.primitives import BitArray\n",
        "\n",
        "# Combine the shots from the individual Trotter circuits\n",
        "bit_array = BitArray.concatenate_shots(\n",
        "    [result.data.meas for result in job.result()]\n",
        ")\n",
        "\n",
        "plot_histogram(bit_array.get_counts(), number_to_keep=20)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2aa74455-d16b-4ac3-a354-54a79d5c5759",
      "metadata": {},
      "source": [
        "### Step 4: Post-process and return result to desired classical format\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d3f7713c-a4e6-407f-94de-0e5cfeb1134c",
      "metadata": {},
      "source": [
        "Now, we run the SQD algorithm using the `diagonalize_fermionic_hamiltonian` function. See the [API documentation](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) for explanations of the arguments to this function.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "7609d1e1-e8ef-48e1-a965-97927f403163",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -13.4222953188441\n",
            "\t\tSubspace dimension: 529\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -13.42237556285828\n",
            "\t\tSubspace dimension: 784\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -13.422045397387413\n",
            "\t\tSubspace dimension: 529\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -13.422379583305478\n",
            "\t\tSubspace dimension: 900\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -13.422376197704326\n",
            "\t\tSubspace dimension: 841\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -13.422421162849295\n",
            "\t\tSubspace dimension: 1089\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -13.422421164670345\n",
            "\t\tSubspace dimension: 1156\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -13.422421492737689\n",
            "\t\tSubspace dimension: 1156\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -13.422421205869572\n",
            "\t\tSubspace dimension: 1156\n",
            "Iteration 4\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -13.422421494558726\n",
            "\t\tSubspace dimension: 1225\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -13.422421492737689\n",
            "\t\tSubspace dimension: 1156\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -13.422421492737689\n",
            "\t\tSubspace dimension: 1156\n"
          ]
        }
      ],
      "source": [
        "from qiskit_addon_sqd.fermion import (\n",
        "    SCIResult,\n",
        "    diagonalize_fermionic_hamiltonian,\n",
        ")\n",
        "\n",
        "# List to capture intermediate results\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        print(f\"\\tSubsample {i}\")\n",
        "        print(f\"\\t\\tEnergy: {result.energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "rng = np.random.default_rng(24)\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    h1e_momentum,\n",
        "    h2e_momentum,\n",
        "    bit_array,\n",
        "    samples_per_batch=100,\n",
        "    norb=norb,\n",
        "    nelec=nelec,\n",
        "    num_batches=3,\n",
        "    max_iterations=5,\n",
        "    symmetrize_spin=True,\n",
        "    callback=callback,\n",
        "    seed=rng,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3dee9c61-fc42-48e9-8888-af8fc831cd5c",
      "metadata": {},
      "source": [
        "The following code cell plots the results. The first plot shows the computed energy as a function of the number of configuration recovery iterations, and the second plot shows the average occupancy of each spatial orbital after the final iteration. Since this is such a small problem, the first iteration already brings us very close to the exact energy (note the scale of the y axis).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "b6879566-8bf5-4c28-bfb6-b2686692e3d3",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Reference energy: -13.42249\n",
            "SQD energy: -13.42242\n",
            "Absolute error: 0.00007\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-krylov-quantum-diagonalization/extracted-outputs/b6879566-8bf5-4c28-bfb6-b2686692e3d3-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "min_es = [\n",
        "    min(result, key=lambda res: res.energy).energy\n",
        "    for result in result_history\n",
        "]\n",
        "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n",
        "\n",
        "# Data for energies plot\n",
        "x1 = range(len(result_history))\n",
        "\n",
        "# Data for avg spatial orbital occupancy\n",
        "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
        "x2 = range(len(y2))\n",
        "\n",
        "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
        "\n",
        "# Plot energies\n",
        "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n",
        "axs[0].set_xticks(x1)\n",
        "axs[0].set_xticklabels(x1)\n",
        "axs[0].axhline(\n",
        "    y=reference_energy,\n",
        "    color=\"#BF5700\",\n",
        "    linestyle=\"--\",\n",
        "    label=\"reference energy\",\n",
        ")\n",
        "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n",
        "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
        "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n",
        "axs[0].legend()\n",
        "\n",
        "# Plot orbital occupancy\n",
        "axs[1].bar(x2, y2, width=0.8)\n",
        "axs[1].set_xticks(x2)\n",
        "axs[1].set_xticklabels(x2)\n",
        "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
        "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
        "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
        "\n",
        "print(f\"Reference energy: {reference_energy:.5f}\")\n",
        "print(f\"SQD energy: {min_e:.5f}\")\n",
        "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0c9f4976-d770-426a-822e-e6756c2cfbe7",
      "metadata": {},
      "source": [
        "### Verify the energy\n",
        "\n",
        "The energy returned by SQD is guaranteed to be an upper bound to the true ground state energy. The value of the energy can be verified because SQD also returns the coefficients of the state vector approximating the ground state. You can compute the energy from the state vector using its one- and two-particle reduced density matrices, as demonstrated in the following code cell.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "e2b9de72-61cf-49d3-a1b5-f043e4b16956",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Recomputed energy: -13.42242\n"
          ]
        }
      ],
      "source": [
        "rdm1 = result.sci_state.rdm(rank=1, spin_summed=True)\n",
        "rdm2 = result.sci_state.rdm(rank=2, spin_summed=True)\n",
        "\n",
        "energy = np.sum(h1e_momentum * rdm1) + 0.5 * np.sum(h2e_momentum * rdm2)\n",
        "\n",
        "print(f\"Recomputed energy: {energy:.5f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5221928a-79ff-4f54-90b4-1fe4d7739aae",
      "metadata": {},
      "source": [
        "## Large-scale hardware example\n",
        "\n",
        "Now, we run a larger example on a real QPU.\n",
        "For the reference energy, we use the results of a [DMRG](https://en.wikipedia.org/wiki/Density_matrix_renormalization_group) calculation that was performed separately.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "933037d8-847e-4986-80da-5ac8d677b2ff",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Using backend ibm_boston\n",
            "Iteration 1\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -28.63965951544449\n",
            "\t\tSubspace dimension: 9801\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -28.625588929202006\n",
            "\t\tSubspace dimension: 9409\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -28.647371834135498\n",
            "\t\tSubspace dimension: 8281\n",
            "Iteration 2\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -28.67213260849567\n",
            "\t\tSubspace dimension: 29584\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -28.670340686158816\n",
            "\t\tSubspace dimension: 27225\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -28.669976379525988\n",
            "\t\tSubspace dimension: 31329\n",
            "Iteration 3\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -28.68622875601382\n",
            "\t\tSubspace dimension: 36100\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -28.698569623143126\n",
            "\t\tSubspace dimension: 34225\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -28.694848533971882\n",
            "\t\tSubspace dimension: 33856\n",
            "Iteration 4\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -28.69883392844593\n",
            "\t\tSubspace dimension: 42025\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -28.701289495200996\n",
            "\t\tSubspace dimension: 38025\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -28.699319594978245\n",
            "\t\tSubspace dimension: 45369\n",
            "Iteration 5\n",
            "\tSubsample 0\n",
            "\t\tEnergy: -28.701936886834154\n",
            "\t\tSubspace dimension: 51076\n",
            "\tSubsample 1\n",
            "\t\tEnergy: -28.702468711812013\n",
            "\t\tSubspace dimension: 53824\n",
            "\tSubsample 2\n",
            "\t\tEnergy: -28.702298147575938\n",
            "\t\tSubspace dimension: 52900\n",
            "Reference energy: -28.70660\n",
            "SQD energy: -28.70247\n",
            "Absolute error: 0.00413\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/sample-based-krylov-quantum-diagonalization/extracted-outputs/933037d8-847e-4986-80da-5ac8d677b2ff-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "\n",
        "# Model parameters\n",
        "norb = 20\n",
        "nelec = (norb // 2, norb // 2)\n",
        "n_bath = norb - 1\n",
        "hybridization = 1.0\n",
        "hopping = 1.0\n",
        "onsite = 10.0\n",
        "chemical_potential = -0.5 * onsite\n",
        "\n",
        "# Generate Hamiltonian and orbital rotation\n",
        "h1e, h2e = siam_hamiltonian(\n",
        "    norb=norb,\n",
        "    hopping=hopping,\n",
        "    onsite=onsite,\n",
        "    hybridization=hybridization,\n",
        "    chemical_potential=chemical_potential,\n",
        ")\n",
        "orbital_rotation = momentum_basis(norb)\n",
        "h1e_momentum, h2e_momentum = rotated(h1e, h2e, orbital_rotation.T.conj())\n",
        "impurity_index = n_bath // 2\n",
        "\n",
        "# Set reference energy to DMRG value computed separately\n",
        "reference_energy = -28.70659686\n",
        "\n",
        "# Algorithm parameters\n",
        "time_step = 0.2\n",
        "krylov_dim = 8\n",
        "\n",
        "# Construct circuits\n",
        "qubits = QuantumRegister(2 * norb, name=\"q\")\n",
        "circuit = QuantumCircuit(qubits)\n",
        "for instruction in prepare_initial_state(qubits, norb=norb, nocc=norb // 2):\n",
        "    circuit.append(instruction)\n",
        "circuit.measure_all()\n",
        "circuits = [circuit.copy()]\n",
        "one_body_evolution = scipy.linalg.expm(-1j * time_step * h1e_momentum)\n",
        "for i in range(krylov_dim - 1):\n",
        "    circuit.remove_final_measurements()\n",
        "    for instruction in trotter_step(\n",
        "        qubits,\n",
        "        time_step,\n",
        "        one_body_evolution,\n",
        "        h2e_momentum,\n",
        "        impurity_index,\n",
        "        norb,\n",
        "    ):\n",
        "        circuit.append(instruction)\n",
        "    circuit.measure_all()\n",
        "    circuits.append(circuit.copy())\n",
        "\n",
        "# Initialize hardware backend\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(\n",
        "    operational=True, simulator=False, min_num_qubits=127\n",
        ")\n",
        "print(f\"Using backend {backend.name}\")\n",
        "\n",
        "# Transpile to backend\n",
        "pass_manager = generate_preset_pass_manager(\n",
        "    optimization_level=3, backend=backend\n",
        ")\n",
        "isa_circuits = pass_manager.run(circuits)\n",
        "\n",
        "# Sample from the circuits\n",
        "sampler = Sampler(backend)\n",
        "sampler.options.environment.job_tags = [\"TUT_SKQD\"]\n",
        "job = sampler.run(isa_circuits, shots=500)\n",
        "\n",
        "# Combine the shots from the individual Trotter circuits\n",
        "bit_array = BitArray.concatenate_shots(\n",
        "    [result.data.meas for result in job.result()]\n",
        ")\n",
        "\n",
        "# Run configuration recovery and diagonalization\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        print(f\"\\tSubsample {i}\")\n",
        "        print(f\"\\t\\tEnergy: {result.energy}\")\n",
        "        print(\n",
        "            f\"\\t\\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "\n",
        "\n",
        "rng = np.random.default_rng(24)\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    h1e_momentum,\n",
        "    h2e_momentum,\n",
        "    bit_array,\n",
        "    samples_per_batch=100,\n",
        "    norb=norb,\n",
        "    nelec=nelec,\n",
        "    num_batches=3,\n",
        "    max_iterations=5,\n",
        "    symmetrize_spin=True,\n",
        "    callback=callback,\n",
        "    seed=rng,\n",
        ")\n",
        "\n",
        "\n",
        "# Plot results\n",
        "min_es = [\n",
        "    min(result, key=lambda res: res.energy).energy\n",
        "    for result in result_history\n",
        "]\n",
        "min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])\n",
        "x1 = range(len(result_history))\n",
        "y2 = np.sum(result.orbital_occupancies, axis=0)\n",
        "x2 = range(len(y2))\n",
        "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n",
        "axs[0].plot(x1, min_es, label=\"energy\", marker=\"o\")\n",
        "axs[0].set_xticks(x1)\n",
        "axs[0].set_xticklabels(x1)\n",
        "axs[0].axhline(\n",
        "    y=reference_energy,\n",
        "    color=\"#BF5700\",\n",
        "    linestyle=\"--\",\n",
        "    label=\"reference energy\",\n",
        ")\n",
        "axs[0].set_title(\"Approximated Ground State Energy vs SQD Iterations\")\n",
        "axs[0].set_xlabel(\"Iteration Index\", fontdict={\"fontsize\": 12})\n",
        "axs[0].set_ylabel(\"Energy\", fontdict={\"fontsize\": 12})\n",
        "axs[0].legend()\n",
        "axs[1].bar(x2, y2, width=0.8)\n",
        "axs[1].set_xticks(x2)\n",
        "axs[1].set_xticklabels(x2)\n",
        "axs[1].set_title(\"Avg Occupancy per Spatial Orbital\")\n",
        "axs[1].set_xlabel(\"Orbital Index\", fontdict={\"fontsize\": 12})\n",
        "axs[1].set_ylabel(\"Avg Occupancy\", fontdict={\"fontsize\": 12})\n",
        "print(f\"Reference energy: {reference_energy:.5f}\")\n",
        "print(f\"SQD energy: {min_e:.5f}\")\n",
        "print(f\"Absolute error: {abs(min_e - reference_energy):.5f}\")\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "482ebea3-84b8-471b-bddc-b282e23162ad",
      "metadata": {},
      "source": [
        "## Next steps\n",
        "\n",
        "<Admonition type=\"tip\" title=\"Recommendations\">\n",
        "  If you found this work interesting, you might be interested in the following material:\n",
        "\n",
        "  * [Sample-based quantum diagonalization of a chemistry Hamiltonian](/docs/tutorials/sample-based-quantum-diagonalization) - a related tutorial using a heuristic variational ansatz instead of Trotter circuits\n",
        "  * [Krylov quantum diagonalization of lattice Hamiltonians](/docs/tutorials/krylov-quantum-diagonalization) - a tutorial on the KQD method\n",
        "  * [SQD addon API documentation](https://qiskit.github.io/qiskit-addon-sqd/apidocs/qiskit_addon_sqd.fermion.html#qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian) - reference for the `diagonalize_fermionic_hamiltonian` function\n",
        "  * [*Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization*](https://arxiv.org/abs/2501.09702) - the paper this tutorial is based on\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "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"
    },
    "hours": 1.5,
    "qpuSeconds": 9
  },
  "nbformat": 4,
  "nbformat_minor": 4
}