diff --git a/docs/tutorials/multi-product-formula.ipynb b/docs/tutorials/multi-product-formula.ipynb index d2628acee87..97595a2f800 100644 --- a/docs/tutorials/multi-product-formula.ipynb +++ b/docs/tutorials/multi-product-formula.ipynb @@ -2,153 +2,262 @@ "cells": [ { "cell_type": "markdown", - "id": "9857bace", + "id": "454a9dfd", "metadata": {}, "source": [ "---\n", "title: Multi-product formulas to reduce Trotter error\n", - "description: Multi-product formulas can be used in observable estimation to reduce Trotter error or to implement time-evolution with fixed Trotter error in lower depth.\n", + "description: Use multi-product formulas in observable estimation to reduce Trotter error or implement time-evolution with fixed Trotter error in lower depth.\n", "---\n", "\n", + "{/* cspell:ignore ncol circo Layerwise markersize unbiasedness infty ndash lesssim propto tenpy unfused Néel correlator Neel exponentiating gtrsim */}\n", "\n", - "{/* cspell:ignore ncol circo Layerwise */}\n", - "{/* cspell:ignore markersize */}\n", - "\n", - "# Multi-product formulas to reduce Trotter error" + "# Multi-product formulas to reduce Trotter error\n", + "*Usage estimate: Four minutes on a Heron r2 processor (NOTE: This is an estimate only. Your runtime may vary.)*" ] }, { "cell_type": "markdown", - "id": "d08a5e51", + "id": "c4d0b2f2", + "metadata": {}, + "source": [ + "## Learning outcomes\n", + "After completing this tutorial, you can expect to understand the following information:\n", + "- How multi-product formulas (MPFs) reduce Trotter error in Hamiltonian simulation by combining expectation values from multiple shallow circuits\n", + "- When MPFs are beneficial over standard product formulas and when they are not the right tool\n", + "- How to compute static and dynamic MPF coefficients using the `qiskit_addon_mpf` package\n", + "- How to execute an MPF workflow end-to-end on IBM Quantum® hardware, including transpilation, error mitigation, and post-processing" + ] + }, + { + "cell_type": "markdown", + "id": "f5dfd316", + "metadata": {}, + "source": [ + "## Prerequisites\n", + "We suggest that users are familiar with the following topics before going through this tutorial:\n", + "- [Compilation methods for Hamiltonian simulation circuits](/docs/tutorials/compilation-methods-for-hamiltonian-simulation-circuits) — introduces Trotter (product formula) circuits in Qiskit.\n", + "- Product formulas in Qiskit, in particular the [`SuzukiTrotter`](/docs/api/qiskit/qiskit.synthesis.SuzukiTrotter) and [`LieTrotter`](/docs/api/qiskit/qiskit.synthesis.LieTrotter) synthesis classes.\n", + "- [Qiskit primitives and the Estimator interface](/docs/guides/primitives)." + ] + }, + { + "cell_type": "markdown", + "id": "07273b26", "metadata": {}, "source": [ - "*Estimated QPU usage:Four minutes on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n", - "\n", "## Background\n", - "This tutorial demonstrates how to use a Multi-Product Formula (MPF) to achieve a lower Trotter error on our observable compared to the one incurred by the deepest Trotter circuit that we will actually execute.\n", - "MPFs reduce the Trotter error of Hamiltonian dynamics through a weighted combination of several circuit executions. Consider the task of finding observable expectation values for the quantum state\n", - "$\\rho(t)=e^{-i H t} \\rho(0) e^{i H t}$ with the Hamiltonian $H$. One can use Product Formulas (PFs) to approximate the time-evolution $e^{-i H t}$ by doing the following:\n", "\n", - "- Write the Hamiltonian $H$ as $H=\\sum_{a=1}^d F_a,$ where $F_a$ are Hermitian operators such that each corresponding unitary can be efficiently implemented on a quantum device.\n", - "- Approximate the terms $F_a$ that do not commute with each other.\n", + "### What are multi-product formulas?\n", + "\n", + "When simulating quantum systems on a quantum computer, a central task is to approximate the time-evolution operator $e^{-iHt}$ for a Hamiltonian $H$. The standard approach uses *product formulas* (PFs), also known as Trotter-Suzuki decompositions. These decompose $H = \\sum_{a=1}^d F_a$ into terms whose individual unitaries $e^{-iF_a t}$ are efficient to implement, and then approximate the full evolution as an ordered product of these simpler unitaries.\n", + "\n", + "The first-order product formula (Lie-Trotter) is:\n", "\n", - "Then, the first-order PF (Lie-Trotter formula) is:\n", + "$$\n", + "S_1(t) := \\prod_{a=1}^d e^{-i F_a t},\n", + "$$\n", + "\n", + "which incurs a quadratic error: $S_1(t) = e^{-iHt} + \\mathcal{O}(t^2)$. Higher-order symmetric formulas $S_{2\\chi}(t)$, where $\\chi$ labels the order of the symmetric product formula (see Ref. [\\[1\\]](#references)), converge faster as $e^{-iHt} + \\mathcal{O}(t^{2\\chi+1})$, but at the cost of deeper circuits per step.\n", + "\n", + "To reduce the error at a *fixed* order $\\chi$, one usually splits the total evolution time $t$ into $k$ smaller Trotter steps. Each step approximates $e^{-iHt/k}$ by a product formula and the steps are concatenated:\n", "\n", - "$$S_1(t):=\\prod_{a=1}^d e^{-i F_a t},$$\n", + "$$\n", + "e^{-iHt} \\approx \\left[S_{2\\chi}(t/k)\\right]^k.\n", + "$$\n", "\n", - "which has a quadratic error term $S_1(t)=e^{-i H t}+\\mathcal{O}\\left(t^{2}\\right)$. One can also use higher-order PFs (Lie-Trotter-Suzuki formulas), which converge faster, and are defined recursively as:\n", + "For a $2\\chi$-order symmetric formula the residual Trotter error then scales as $\\mathcal{O}\\!\\left(t^{2\\chi+1} / k^{2\\chi}\\right)$. So increasing $k$ rapidly suppresses the Trotter error — but it also linearly deepens the circuit, and on noisy hardware that means more accumulated gate noise. This tension between **Trotter error (favors larger $k$)** and **hardware noise (favors smaller $k$)** is precisely what multi-product formulas are designed to resolve. Note that MPFs are about combining results from *different choices of $k$* at a fixed order $\\chi$ — they do not change the order of the underlying product formula.\n", "\n", - "$$S_2(t):=\\prod_{a=1}^d e^{-i F_a t/2}\\prod_{a=1}^d e^{-i F_a t/2}$$\n", + "**Multi-product formulas (MPFs)** [\\[1\\]](#references) construct a *weighted linear combination* of expectation values obtained from several shallower Trotter circuits, each using a different number of Trotter steps $k_1, k_2, \\ldots, k_r$ (a set of $r$ step counts):\n", "\n", - "$$S_{2 \\chi}(t):= S_{2 \\chi -2}(s_{\\chi}t)^2 S_{2 \\chi -2}((1-4s_{\\chi})t)S_{2 \\chi -2}(s_{\\chi}t)^2,$$\n", + "$$\n", + "\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_{j=1}^r x_j \\, \\langle A \\rangle_{k_j}(t),\n", + "$$\n", "\n", - "where $\\chi$ is the order of the symmetric PF and $s_p = \\left( 4 - 4^{1/(2p-1)} \\right)^{-1}$. For long time-evolutions, one can split the time interval $t$ into $k$ intervals, called Trotter steps, of duration $t/k$ and approximate the time-evolution in each interval with a $\\chi$ order product formula $S_{\\chi}$. Thus, the PF of order $\\chi$ for time-evolution operator over $k$ Trotter steps is:\n", + "where $\\langle A \\rangle_{k_j}(t)$ is the expectation value of an observable $A$ at time $t$ estimated from a Trotter circuit with $k_j$ steps, and the coefficients $\\{x_j\\}_{j=1}^r$ are chosen so that the leading Trotter-error terms in the combination cancel. We will revisit this expression in [Step 4](#small-scale-step-4), where we evaluate it explicitly to combine our Trotter results. The key practical point is that the deepest circuit in the MPF only needs $k_{\\max}$ steps, which is much smaller than the single $k$ that would be required to reach the same effective Trotter error directly. The shallower circuits make the MPF approach better suited to noisy hardware.\n", "\n", - "$$ S_{\\chi}^{k}(t) = \\left[ S_{\\chi} \\left( \\frac{t}{k} \\right)\\right]^k = e^{-i H t}+O\\left(t \\left( \\frac{t}{k} \\right)^{\\chi} \\right)$$\n", + "### How are the coefficients determined?\n", "\n", - "where the error term decreases with the number of Trotter steps $k$ and the order $\\chi$ of the PF.\n", + "There are two families of MPF coefficients:\n", "\n", - "Given an integer $k \\geq 1$ and a product formula $S_{\\chi}(t)$, the approximate time-evolved state $\\rho_k(t)$ can be obtained from $\\rho_0$ by applying $k$ iterations of the product formula $S_{\\chi}\\left(\\frac{t}{k}\\right)$.\n", + "**Static coefficients** are independent of the Hamiltonian, the initial state, and the evolution time. They are found by solving a linear system $Ax = b$ that enforces cancellation of the leading Trotter error terms. For a set of Trotter steps $\\{k_j\\}_{j=1}^r$ used with a $2\\chi$-order symmetric product formula, expanding the Trotter error in inverse powers of $k_j$ leads to constraint equations of the form:\n", "\n", "$$\n", - "\\rho_k(t)=S_{\\chi}\\left(\\frac{t}{k}\\right)^k \\rho_0 S_{\\chi}\\left(\\frac{t}{k}\\right)^{-k}\n", + "\\sum_{j=1}^r x_j = 1, \\quad \\sum_{j=1}^r \\frac{x_j}{k_j^{\\eta_n}} = 0 \\quad (n = 0, \\ldots, r-2),\n", "$$\n", "\n", - "$\\rho_k(t)$ is an approximation for $\\rho(t)$ with the Trotter approximation error ||$\\rho_k(t)-\\rho(t) ||$. If we consider a linear combination of Trotter approximations of $\\rho(t)$:\n", + "where the integer exponents $\\{\\eta_n\\}$ are the orders of the successive Trotter-error terms for the chosen product formula. For a *symmetric* $2\\chi$-order PF, the leading error in $\\left[S_{2\\chi}(t/k)\\right]^k$ scales as $1/k^{2\\chi}$, with subsequent corrections at $1/k^{2\\chi+2}, 1/k^{2\\chi+4}, \\ldots$ — so the exponents are $\\eta_n = 2\\chi + 2n$. For non-symmetric PFs, both odd and even powers contribute and $\\eta_n = 2\\chi + n$. See Ref. [\\[1\\]](#references) for the full derivation. The first equation in the system above ensures unbiasedness (the MPF reproduces the exact expectation value in the $k_j \\to \\infty$ limit), and the remaining $r-1$ equations successively cancel the first $r-1$ Trotter-error terms. When the resulting $L_1$-norm $\\|x\\|_1$ is too large (which amplifies sampling noise), you can instead solve an approximate optimization that caps $\\|x\\|_1$ while minimizing $\\|Ax - b\\|$.\n", "\n", + "**Dynamic coefficients** [\\[2\\]](#references), [\\[3\\]](#references) additionally depend on the Hamiltonian, initial state, and evolution time $t$. They minimize the Frobenius-norm distance between the true time-evolved state and the MPF approximation:\n", "\n", "$$\n", - "\\mu(t) = \\sum_{j}^{l} x_j \\rho^{k_j}_{j}\\left(\\frac{t}{k_j}\\right) + \\text{some remaining Trotter error} \\, ,\n", + "\\|\\rho(t) - \\mu^D(t)\\|_F^2 = 1 + \\sum_{i,j} M_{ij}(t)\\, x_i(t)\\, x_j(t) - 2\\sum_i L_i(t)\\, x_i(t),\n", "$$\n", "\n", - "where $x_j$ are our weighting coefficients, $\\rho^{k_j}_j$ is the density matrix corresponding to the pure state obtained by evolving the initial state with the product formula, $S^{k_j}_{\\chi}$, involving $k_j$ Trotter steps, and $j \\in {1, ..., l}$ indexes the number of PFs that make up the MPF. All the terms in $\\mu(t)$ use the same product formula $S_{\\chi}(t)$ as its base.\n", - "The goal is to improve upon ||$\\rho_k(t)-\\rho(t) \\|$ by finding $\\mu(t)$ with even lower $\\|\\mu(t)-\\rho(t)\\|$.\n", + "where $M_{ij}(t) = \\mathrm{Tr}[\\rho_{k_i}(t)\\,\\rho_{k_j}(t)]$ is the Gram matrix of overlaps between Trotter-evolved states for different step counts $k_i, k_j$, and $L_i(t) = \\mathrm{Tr}[\\rho(t)\\,\\rho_{k_i}(t)]$ measures overlap with the (approximate) exact state. In this tutorial these quantities are computed efficiently using tensor-network methods, specifically the TeNPy-based backends in `qiskit_addon_mpf`.\n", "\n", - "* $\\mu(t)$ needs not be a physical state as $x_i$ need not be positive. The goal here is to minimize the error in the expectation value of the observables and not to find a physical replacement for $\\rho(t)$.\n", - "* $k_j$ determines both the circuit depth and level of Trotter approximation. Smaller values of $k_j$ lead to shorter circuits, which incur fewer circuit errors but will be a less accurate approximation to the desired state.\n", + "### When to use MPFs\n", "\n", + "MPFs are most beneficial when:\n", "\n", - "The key here is that the remaining Trotter error given by $\\mu(t)$ is smaller than the Trotter error that one would obtain by simply using the largest $k_j$ value.\n", + "- **Circuit depth is the bottleneck.** If hardware noise limits the depth you can run, use MPFs to achieve higher effective Trotter accuracy from shallower circuits.\n", + "- **You need accurate expectation values, not full state preparation.** MPFs operate at the level of expectation values — they combine classical numbers, not quantum states. They are therefore ideal for observable estimation when using the Estimator primitive.\n", + "- **You combine a modest number of Trotter-step counts.** Typically combining $r = 3$–$5$ different step counts $k_j$ is sufficient to cancel several leading Trotter-error terms while keeping $\\|x\\|_1$ manageable.\n", "\n", - "You can view the usefulness of this from two perspectives:\n", + "### When MPFs might not help\n", "\n", - "1. For a fixed budget of Trotter steps that you can execute, you can obtain results with a Trotter error that is smaller in total.\n", - "2. Given some target number of Trotter steps that is too large to execute, you can use the MPF to find a collection of lower-depth circuits to run that results in a similar Trotter error." + "- **Very short evolution times.** When $t$ is small enough that a single low-order Trotter formula is already accurate, the overhead of running multiple circuits is unnecessary.\n", + "- **State-preparation tasks.** MPFs produce a corrected *expectation value*, not a corrected quantum state. If you need the actual time-evolved state (for example, as input to another quantum subroutine), MPFs do not apply.\n", + "- **Trotter-step counts that violate the convergence regime.** The static-coefficient derivation expands each individual $\\left[S_{2\\chi}(t/k_j)\\right]^{k_j}$ as a series in $t/k_j$; this expansion only converges well when $t/k_{\\min} \\lesssim 1$. If $k_{\\min}$ is chosen too small for the given $t$, the shallowest circuit is far outside the perturbative regime, the higher-order error terms that the MPF leaves uncanceled become large, and the cancellation can require large coefficients. The $L_1$-norm $\\|x\\|_1$ is the practical diagnostic: when $\\|x\\|_1 \\gg 1$, the sampling overhead $\\propto \\|x\\|_1^2$ might outweigh the Trotter-error reduction. See the [guide on choosing Trotter steps](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) for details.\n", + "\n", + "### What this tutorial covers\n", + "\n", + "This tutorial walks through an end-to-end MPF workflow in two stages. First, a **small-scale simulator example** (10-qubit Heisenberg chain) demonstrates how to set up the problem, compute static and dynamic MPF coefficients, and compare the resulting expectation values against exact diagonalization. Then, a **large-scale hardware example** (50-qubit XXZ chain) shows how to transpile, execute on IBM Quantum hardware with error mitigation, and post-process results using the MPF coefficients. Throughout, we use the `qiskit_addon_mpf` package alongside standard Qiskit tools." ] }, { "cell_type": "markdown", - "id": "152c479f", + "id": "b5d478ce", "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", - "* MPF Qiskit addons (`pip install qiskit_addon_mpf`)\n", - "* Qiskit addons utils (`pip install qiskit_addon_utils`)\n", - "* Quimb library (`pip install quimb`)\n", - "* Qiskit Quimb library (`pip install qiskit-quimb`)\n", - "* Numpy v0.21 for compatibility across packages (`pip install numpy==0.21`)" - ] - }, - { - "cell_type": "markdown", - "id": "32f44dc0", - "metadata": {}, - "source": [ - "## Part I. Small-scale example" + "- Qiskit SDK v2.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n", + "- Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n", + "- Qiskit Aer simulator (`pip install qiskit-aer`)\n", + "- MPF Qiskit addon with the TeNPy backend (`pip install \"qiskit-addon-mpf[tenpy]\"`)\n", + "- Qiskit addon utilities (`pip install qiskit-addon-utils`)\n", + "- SciPy (`pip install scipy`)" ] }, { "cell_type": "markdown", - "id": "717dc6b6", + "id": "2584c37e", "metadata": {}, "source": [ - "### Explore the stability of MPF\n", + "## Setup\n", "\n", - "There is no obvious restriction on the choice of number of Trotter steps $k_j$ that make up the MPF state $\\mu(t)$. However, these must be picked carefully to avoid instabilities in the resulting expectation values calculated from $\\mu(t)$. A good general rule is to set the smallest Trotter step $k_{\\text{min}}$ so that $t/k_{\\text{min}} \\lt 1$. If you want to learn more about this and how to choose your other $k_j$ values, refer to the [How to choose the Trotter steps for an MPF](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) guide.\n", - "\n", - "In the example below we explore the stability of the MPF solution by calculating the expectation value of the magnetization for a range of times using different time-evolved states. Specifically, we compare the expectation values calculated from each of the approximate time-evolutions implemented with the corresponding Trotter steps and the various MPF models (static and dynamic coefficients) with the exact values of the time-evolved observable. First, let's define the parameters for the Trotter formulas and the evolution times" + "Below we collect *all* the package imports used throughout this tutorial in a single cell. We also define a `CollectAndCollapse` transpiler pass that fuses adjacent `rxx` and `ryy` rotations into a single `XXPlusYYGate`. This pass is applied both during circuit construction in Step 1 (to keep the gate count low) and indirectly when we extract the layered structure for the dynamic MPF in Step 4 (TeNPy expects two-qubit gates, not pairs of unfused rotations)." ] }, { "cell_type": "code", "execution_count": 1, - "id": "bcf32912", + "id": "bf79f9e7", "metadata": {}, "outputs": [], "source": [ + "import warnings\n", + "\n", "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from functools import partial\n", + "from copy import deepcopy\n", "\n", - "mpf_trotter_steps = [1, 2, 4]\n", - "order = 2\n", - "symmetric = False\n", + "from qiskit import QuantumCircuit\n", + "from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector\n", + "from qiskit.synthesis import SuzukiTrotter\n", + "from qiskit.transpiler import CouplingMap, PassManager, Target\n", + "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", + "from qiskit.circuit.library import XXPlusYYGate\n", + "from qiskit.transpiler.passes.optimization.collect_and_collapse import (\n", + " CollectAndCollapse,\n", + " collect_using_filter_function,\n", + " collapse_to_operation,\n", + ")\n", "\n", - "trotter_times = np.arange(0.5, 1.55, 0.1)\n", - "exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)" + "from qiskit_aer import AerSimulator\n", + "from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService\n", + "\n", + "from qiskit_addon_utils.problem_generators import (\n", + " generate_xyz_hamiltonian,\n", + " generate_time_evolution_circuit,\n", + ")\n", + "from qiskit_addon_utils.slicing import slice_by_depth\n", + "from qiskit_addon_mpf.static import setup_static_lse\n", + "from qiskit_addon_mpf.dynamic import setup_dynamic_lse\n", + "from qiskit_addon_mpf.costs import (\n", + " setup_exact_problem,\n", + " setup_sum_of_squares_problem,\n", + " setup_frobenius_problem,\n", + ")\n", + "from qiskit_addon_mpf.backends.tenpy_layers import (\n", + " LayerModel,\n", + " LayerwiseEvolver,\n", + ")\n", + "from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state\n", + "\n", + "from scipy.linalg import expm\n", + "\n", + "# Suppress TeNPy's `unit_cell_width` future-API warning. The default\n", + "# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what\n", + "# `CouplingMap.from_line(...)` produces here, so the warning is informational.\n", + "warnings.filterwarnings(\n", + " \"ignore\",\n", + " message=r\".*unit_cell_width.*\",\n", + " category=UserWarning,\n", + ")\n", + "\n", + "\n", + "# --- Helper: collect XX + YY rotations into a single gate ---\n", + "def filter_function(node):\n", + " return node.op.name in {\"rxx\", \"ryy\"}\n", + "\n", + "\n", + "collect_function = partial(\n", + " collect_using_filter_function,\n", + " filter_function=filter_function,\n", + " split_blocks=True,\n", + " min_block_size=1,\n", + ")\n", + "\n", + "\n", + "def collapse_to_xx_plus_yy(block):\n", + " param = 0.0\n", + " for node in block.data:\n", + " param += node.operation.params[0]\n", + " return XXPlusYYGate(param)\n", + "\n", + "\n", + "collapse_function = partial(\n", + " collapse_to_operation,\n", + " collapse_function=collapse_to_xx_plus_yy,\n", + ")\n", + "\n", + "pm = PassManager()\n", + "pm.append(CollectAndCollapse(collect_function, collapse_function))" + ] + }, + { + "cell_type": "markdown", + "id": "24f08467", + "metadata": {}, + "source": [ + "## Small-scale simulator example" ] }, { "cell_type": "markdown", - "id": "095fa6cb", + "id": "378e82ba", "metadata": {}, "source": [ - "For this example we will use the Neel state as the initial state $\\vert \\text{Neel} \\rangle = \\vert 0101...01 \\rangle$ and the Heisenberg model on a line of 10 sites for the Hamiltonian governing the time-evolution\n", + "### Step 1: Map classical inputs to a quantum problem\n", + "\n", + "We begin with a 10-qubit Heisenberg model on a line, using the Néel state $\\vert 0101\\ldots01 \\rangle$ as the initial state. The Hamiltonian is:\n", "\n", "$$\n", - "\\hat{\\mathcal{H}}_{Heis} = J \\sum_{i=1}^{L-1} \\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ Z_i Z_{(i+1)} \\right) \\, ,\n", + "\\hat{\\mathcal{H}}_{\\text{Heis}} = J \\sum_{i=1}^{L-1} \\left(X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1}\\right),\n", "$$\n", "\n", - "where $J$ is the coupling strength for nearest-neighbor edges." + "where $J$ is the nearest-neighbor coupling strength. We measure the ZZ correlator $Z_{L/2-1} Z_{L/2}$ on a pair of qubits in the middle of the chain, and use Trotter steps $k_j = [1, 2, 4]$ with a second-order product formula." ] }, { "cell_type": "code", "execution_count": 2, - "id": "3a20e08d", + "id": "bdd0d4fc", "metadata": {}, "outputs": [ { @@ -163,41 +272,23 @@ } ], "source": [ - "from qiskit.transpiler import CouplingMap\n", - "from rustworkx.visualization import graphviz_draw\n", - "from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian\n", - "import numpy as np\n", - "\n", - "\n", "L = 10\n", "\n", - "# Generate some coupling map to use for this example\n", + "# Generate coupling map and Hamiltonian\n", "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n", - "graphviz_draw(coupling_map.graph, method=\"circo\")\n", "\n", - "# Get a qubit operator describing the Heisenberg field model\n", "hamiltonian = generate_xyz_hamiltonian(\n", " coupling_map,\n", " coupling_constants=(1.0, 1.0, 1.0),\n", " ext_magnetic_field=(0.0, 0.0, 0.0),\n", ")\n", - "\n", - "\n", "print(hamiltonian)" ] }, - { - "cell_type": "markdown", - "id": "69438c8c", - "metadata": {}, - "source": [ - "The observable that we will be measuring is magnetization on a pair of qubits in the middle of the chain." - ] - }, { "cell_type": "code", - "execution_count": null, - "id": "41fc0bd1", + "execution_count": 3, + "id": "fd3dc9c8", "metadata": {}, "outputs": [ { @@ -210,89 +301,46 @@ } ], "source": [ - "from qiskit.quantum_info import SparsePauliOp\n", - "\n", + "# Observable: ZZ on the middle pair of qubits\n", "observable = SparsePauliOp.from_sparse_list(\n", " [(\"ZZ\", (L // 2 - 1, L // 2), 1.0)], num_qubits=L\n", ")\n", "print(observable)" ] }, - { - "cell_type": "markdown", - "id": "13c67498", - "metadata": {}, - "source": [ - "We define a transpiler pass to collect the XX and YY rotations in the circuit as a single XX+YY gate. This will allow us to leverage TeNPy's spin conservation properties during the MPO computation, significantly speeding up the calculation." - ] - }, { "cell_type": "code", "execution_count": 4, - "id": "b0eec735", + "id": "398c33b2", "metadata": {}, "outputs": [], "source": [ - "from qiskit.circuit.library import XXPlusYYGate\n", - "from qiskit.transpiler import PassManager\n", - "from qiskit.transpiler.passes.optimization.collect_and_collapse import (\n", - " CollectAndCollapse,\n", - " collect_using_filter_function,\n", - " collapse_to_operation,\n", - ")\n", - "from functools import partial\n", - "\n", - "\n", - "def filter_function(node):\n", - " return node.op.name in {\"rxx\", \"ryy\"}\n", - "\n", - "\n", - "collect_function = partial(\n", - " collect_using_filter_function,\n", - " filter_function=filter_function,\n", - " split_blocks=True,\n", - " min_block_size=1,\n", - ")\n", - "\n", - "\n", - "def collapse_to_xx_plus_yy(block):\n", - " param = 0.0\n", - " for node in block.data:\n", - " param += node.operation.params[0]\n", - " return XXPlusYYGate(param)\n", - "\n", - "\n", - "collapse_function = partial(\n", - " collapse_to_operation,\n", - " collapse_function=collapse_to_xx_plus_yy,\n", - ")\n", + "# MPF parameters\n", + "mpf_trotter_steps = [1, 2, 4]\n", + "order = 2\n", + "symmetric = False\n", "\n", - "pm = PassManager()\n", - "pm.append(CollectAndCollapse(collect_function, collapse_function))" + "trotter_times = np.arange(0.5, 1.55, 0.1)\n", + "exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)" ] }, { "cell_type": "markdown", - "id": "46470a22", + "id": "1512ba0e", "metadata": {}, "source": [ - "Then we create the circuits implementing the approximate Trotter time-evolutions." + "#### Build Trotter circuits\n", + "\n", + "We create the circuits implementing the approximate Trotter time-evolutions for each time point and each Trotter step count. The `CollectAndCollapse` pass defined in the Setup section collects XX and YY rotations into single XX+YY gates, to prepare for more efficient tensor-network simulation later." ] }, { "cell_type": "code", "execution_count": 5, - "id": "8d7b9367", + "id": "1c194d2b", "metadata": {}, "outputs": [], "source": [ - "from qiskit.synthesis import SuzukiTrotter\n", - "from qiskit_addon_utils.problem_generators import (\n", - " generate_time_evolution_circuit,\n", - ")\n", - "from qiskit import QuantumCircuit\n", - "\n", - "\n", "# Initial Neel state preparation\n", "initial_state_circ = QuantumCircuit(L)\n", "initial_state_circ.x([i for i in range(L) if i % 2 != 0])\n", @@ -322,13 +370,13 @@ { "cell_type": "code", "execution_count": 6, - "id": "92dc20a7", + "id": "c7ee61e7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "\"Output" + "\"Output" ] }, "execution_count": 6, @@ -342,66 +390,84 @@ }, { "cell_type": "markdown", - "id": "aa95b6c8", + "id": "4cd6c782", "metadata": {}, "source": [ - "Next, we calculate the time-evolved expectation values from the Trotter circuits." + "### Step 2: Optimize problem for quantum hardware execution\n", + "\n", + "For the small-scale example we target the Aer simulator. Two transformations happen before the circuits are ready to execute:\n", + "\n", + "1. **Gate collection at the Hamiltonian-simulation level.** In the Setup cell we built a `CollectAndCollapse` pass that fuses adjacent `rxx` and `ryy` rotations into a single `XXPlusYYGate`. We already applied this pass when we built the Trotter circuits in Step 1 (the `pm.run(...)` call). This both reduces the two-qubit gate count and produces a structure that is more amenable to tensor-network simulation for the dynamic-coefficient computation later.\n", + "\n", + "2. **Lowering to the simulator's ISA.** Below we run Qiskit's preset pass manager at `optimization_level=3` to lower each Trotter circuit to the simulator's instruction-set architecture (ISA)." ] }, { "cell_type": "code", - "execution_count": 24, - "id": "69f60701", + "execution_count": 7, + "id": "03590a05", "metadata": {}, "outputs": [], "source": [ - "from copy import deepcopy\n", - "from qiskit_aer import AerSimulator\n", - "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n", - "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", + "aer_sim = AerSimulator()\n", + "pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)\n", "\n", + "isa_circs_all_times = [\n", + " pm_sim.run([deepcopy(c) for c in mpf_circuits])\n", + " for mpf_circuits in all_circs\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "ab6588d3", + "metadata": {}, + "source": [ + "### Step 3: Execute using Qiskit primitives\n", "\n", - "aer_sim = AerSimulator()\n", + "For the small-scale example we run the ISA-lowered Trotter circuits through the `EstimatorV2` primitive backed by Aer. Doing so gives us a *noiseless* reference value for each $(k_j, t)$ pair — these are the $\\langle A \\rangle_{k_j}(t)$ values that the MPF will combine in Step 4. We sweep over evolution times so that we can later plot the full time-series curve of each individual product formula and of the MPF." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "7225d782", + "metadata": {}, + "outputs": [], + "source": [ "estimator = Estimator(mode=aer_sim)\n", "\n", "mpf_expvals_all_times, mpf_stds_all_times = [], []\n", - "for t, mpf_circuits in zip(trotter_times, all_circs):\n", - " mpf_expvals = []\n", - " circuits = [deepcopy(circuit) for circuit in mpf_circuits]\n", - " pm_sim = generate_preset_pass_manager(\n", - " backend=aer_sim, optimization_level=3\n", - " )\n", - " isa_circuits = pm_sim.run(circuits)\n", + "for isa_circuits in isa_circs_all_times:\n", " result = estimator.run(\n", " [(circuit, observable) for circuit in isa_circuits], precision=0.005\n", " ).result()\n", - " mpf_expvals = [res.data.evs for res in result]\n", - " mpf_stds = [res.data.stds for res in result]\n", - " mpf_expvals_all_times.append(mpf_expvals)\n", - " mpf_stds_all_times.append(mpf_stds)" + " mpf_expvals_all_times.append([res.data.evs for res in result])\n", + " mpf_stds_all_times.append([res.data.stds for res in result])" ] }, { "cell_type": "markdown", - "id": "b84e4000", + "id": "a384f017", "metadata": {}, "source": [ - "We also calculate the exact expectation values for comparison." + "\n", + "### Step 4: Post-process and return result in desired classical format\n", + "\n", + "Step 4 is where the MPF is actually constructed. Even though the coefficients $x_j$ are *computed* here (and, for the dynamic variant, this computation can be intensive), conceptually they are a classical recipe for combining the quantum measurements from Step 3 into a single corrected expectation value — so we treat the whole coefficient and combination workflow as post-processing.\n", + "\n", + "To assess how well the MPF tracks the true dynamics, we first compute the exact time-evolved expectation values by directly exponentiating the Hamiltonian. This is only tractable because $L = 10$; in the large-scale hardware example below we will have to rely on tensor-network estimates instead." ] }, { "cell_type": "code", - "execution_count": null, - "id": "32ecef77", + "execution_count": 9, + "id": "a223a360", "metadata": {}, "outputs": [], "source": [ - "from scipy.linalg import expm\n", - "from qiskit.quantum_info import Statevector\n", - "\n", "exact_expvals = []\n", "for t in exact_evolution_times:\n", - " # Exact expectation values\n", " exp_H = expm(-1j * t * hamiltonian.to_matrix())\n", " initial_state = Statevector(initial_state_circ).data\n", " time_evolved_state = exp_H @ initial_state\n", @@ -416,132 +482,67 @@ }, { "cell_type": "markdown", - "id": "949130d1", + "id": "361e726e", "metadata": {}, "source": [ "#### Static MPF coefficients\n", "\n", - "Static MPFs are those where the $x_j$ values do not depend on the evolution time, $t$. Let's consider the order $\\chi = 1$ PF with $k_j$ Trotter steps, this can be written as:\n", - "\n", - "$$ S_1^{k_j}\\left( \\frac{t}{k_j} \\right)=e^{-i H t}+ \\sum_{n=1}^{\\infty} A_n \\frac{t^{n+1}}{k_j^n} $$\n", - "\n", - "where $A_n$ are matrices which depend on the commutators of $F_a$ terms in the decomposition of the Hamiltonian. It is important to note that $A_n$ themselves are independent of time and the number of Trotter steps $k_j$. Therefore, it is possible to cancel out lower-order error terms contributing to $\\mu(t)$ with a careful choice of the weights $x_j$ of the linear combination. To cancel the Trotter error for the first $l-1$ terms (these will give the largest contributions as they correspond to the lower number of Trotter steps) in the expression for $\\mu(t)$, the coefficients $x_j$ must satisfy the following equations:\n", - "\n", - "$$ \\sum_{j=1}^l x_j = 1 $$\n", - "$$ \\sum_{j=1}^{l-1} \\frac{x_j}{k_j^{n}} = 0 $$\n", - "\n", - "with $n=0, ... l-2$. The first equation guarantees that there is no bias in the constructed state $\\mu(t)$, while the second equation ensures the cancellation of the Trotter errors. For higher-order PF, the second equation becomes $ \\sum_{j=1}^{l-1} \\frac{x_j}{k_j^{\\eta}} = 0 $ where $\\eta = \\chi + 2n$ for symmetric PFs and $\\eta = \\chi + n$ otherwise, with $n=0, ..., l-2$. The resulting error (Refs. [\\[1\\]](#references),[\\[2\\]](#references)) is then\n", - "\n", - "$$ \\epsilon = \\mathcal{O} \\left( \\frac{t^{l+1}}{k_1^l} \\right).$$\n", - "\n", - "\n", - "Determining the static MPF coefficients for a given set of $k_j$ values amounts to solving the linear system of equations defined by the two equations above for the variables $x_j$: $Ax=b$. Where $x$ are our coefficients of interest, $A$ is a matrix which depends on $k_j$ and the type of PF we use ($S$), and $b$ is a vector of constraints. Specifically:\n", - "\n", - "$$A_{0,j} = 1 $$\n", - "$$A_{i>0,j} = k_{j}^{-(\\chi + s(i-1))}$$\n", - "$$b_0 = 1$$\n", - "$$b_{i>0} = 0 $$\n", - "\n", + "Static MPFs use coefficients $x_j$ that are independent of the evolution time, the Hamiltonian, and the initial state. We set up the linear system $Ax = b$ described in the Background and solve for the coefficients. The matrix $A$ is determined by the Trotter step counts $k_j$, the order $\\chi$ of the product formula, and whether the formula is symmetric (which controls the exponents $\\eta_n$).\n", "\n", - "where $\\chi$ is the ``order``, $s$ is $2$ if ``symmetric`` is ``True`` and $1$ otherwise, $k_{j}$ are the ``trotter_steps``, and $x$ are the variables to solve for. The indices $i$ and $j$ start at $0$. We can also visualize this in matrix form:\n", + "For our small-scale example we use $k_j = [1, 2, 4]$ with a non-symmetric order-$2\\chi=2$ Suzuki-Trotter formula (so $\\chi=1$ and $\\eta_n = 2 + n$, giving $\\eta_0 = 2,\\, \\eta_1 = 3$). The system becomes:\n", "\n", "$$\n", "A =\n", "\\begin{bmatrix}\n", - "A_{0,0} & A_{0,1} & A_{0,2} & ... \\\\\n", - "A_{1,0} & A_{1,1} & A_{1,2} & ... \\\\\n", - "A_{2,0} & A_{2,1} & A_{2,2} & ... \\\\\n", - "... & ... & ... & ...\n", - "\\end{bmatrix} =\n", - "\\begin{bmatrix}\n", - "1 & 1 & 1 & ... \\\\\n", - "k_{0}^{-(\\chi + s(1-1))} & k_{1}^{-(\\chi + s(1-1))} & k_{2}^{-(\\chi + s(1-1))} & ... \\\\\n", - "k_{0}^{-(\\chi + s(2-1))} & k_{1}^{-(\\chi + s(2-1))} & k_{2}^{-(\\chi + s(2-1))} & ... \\\\\n", - "... & ... & ... & ...\n", - "\\end{bmatrix}\n", - "$$\n", - "\n", - "and\n", - "\n", - "$$\n", + "1 & 1 & 1\\\\\n", + "1 & \\frac{1}{2^2} & \\frac{1}{4^2} \\\\\n", + "1 & \\frac{1}{2^3} & \\frac{1}{4^3} \\\\\n", + "\\end{bmatrix}, \\quad\n", "b =\n", "\\begin{bmatrix}\n", - "b_{0} \\\\\n", - "b_{1} \\\\\n", - "b_{2} \\\\\n", - "...\n", - "\\end{bmatrix} =\n", - "\\begin{bmatrix}\n", "1 \\\\\n", "0 \\\\\n", - "0 \\\\\n", - "...\n", - "\\end{bmatrix}\n", + "0\n", + "\\end{bmatrix}.\n", "$$\n", "\n", - "\n", - "\n", - "For more details, refer to the documentation of the Linear System of Equations ([LSE](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.LSE.html)).\n", - "\n", - "We can find a solution for $x$ analytically as $x = A^{-1}b$; see for example Refs. [\\[1\\]](#references) or [\\[2\\]](#references).\n", - "However, this exact solution can be \"ill-conditioned\", resulting in very large L1-norms of our coefficients, $x$, which can lead to bad performance of the MPF.\n", - "Instead, one can also obtain an approximate solution that minimizes the L1-norm of $x$ in order to attempt to optimize the MPF behavior." + "The first row enforces unbiasedness ($\\sum_j x_j = 1$); the second and third rows cancel the leading $1/k^2$ and next-order $1/k^3$ Trotter-error terms, respectively." ] }, { "cell_type": "markdown", - "id": "f16ac209", + "id": "4f2ca1e2", "metadata": {}, "source": [ "##### Set up the LSE\n", "\n", - "Now that we have chosen our $k_j$ values, we must first construct the LSE, $Ax=b$ as explained above.\n", - "The matrix $A$ depends not only on $k_j$ but also our choice of PF, in particular its _order_.\n", - "Additionally, you might take into account whether the PF is symmetric or not (see [\\[1\\]](#references)) by setting `symmetric=True/False`.\n", - "However, this is not required, as shown by Ref. [\\[2\\]](#references)." + "We use `setup_static_lse` from `qiskit_addon_mpf.static` to assemble the matrix $A$ and right-hand-side vector $b$ described above. The matrix $A$ depends not only on $k_j$ but also on our choice of product formula — in particular its *order* $\\chi$ and whether it is *symmetric*. The `symmetric` flag controls the exponent pattern $\\eta_n$ (symmetric formulas only produce even-power Trotter-error terms; see Ref. [\\[1\\]](#references)). Note that, as shown in Ref. [\\[2\\]](#references), setting `symmetric=True` is not strictly necessary even when the underlying PF is symmetric — the non-symmetric LSE remains valid (it enforces extra unneeded constraints).\n", + "\n", + "For our example we have already set `order = 2` and `symmetric = False` in Step 1." ] }, { "cell_type": "code", - "execution_count": 9, - "id": "509bee67", + "execution_count": 10, + "id": "827b0b42", "metadata": {}, "outputs": [], "source": [ - "from qiskit_addon_mpf.static import setup_static_lse\n", - "\n", "lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)" ] }, { "cell_type": "markdown", - "id": "d6f76196", + "id": "003e4bdd", "metadata": {}, "source": [ - "Let's work through the values chosen above to construct the $A$ matrix and the $b$ vector. With $j=0,1, 2$ Trotter steps $k_j = [1, 2, 4]$, order $\\chi = 2$ and choice of non-symmetric Trotter steps ($s=1$), we have that the matrix elements of $A$ below the first row are determined by the expression $A_{i>0,j} = k_{j}^{-(2 + 1(i-1))}$, specifically:\n", - "\n", - "$$ A_{0,0} = A_{0,1} = A_{0,2} = 1 $$\n", - "$$ A_{1,j} = k_{j}^{-1} \\rightarrow A_{1,0} = \\frac{1}{1^2}, \\;, A_{1,1} = \\frac{1}{2^2}, \\;, A_{1,2} = \\frac{1}{4^2}$$\n", - "$$ A_{2,j} = k_{j}^{-2} \\rightarrow A_{2,0} = \\frac{1}{1^3}, \\;, A_{2,1} = \\frac{1}{2^3}, \\;, A_{2,2} = \\frac{1}{4^3}$$\n", - "\n", - "or in matrix form:\n", - "\n", - "$$\n", - "A =\n", - "\\begin{bmatrix}\n", - "1 & 1 & 1\\\\\n", - "1 & \\frac{1}{2^2} & \\frac{1}{4^2} \\\\\n", - "1 & \\frac{1}{2^3} & \\frac{1}{4^3} \\\\\n", - "\\end{bmatrix}\n", - "$$\n", - "\n", - "This is possible to see by inspecting the `lse` object:" + "Inspect the constructed matrix $A$ and vector $b$ to confirm they match the system written above." ] }, { "cell_type": "code", - "execution_count": 10, - "id": "4982ce27", + "execution_count": 11, + "id": "6f879978", "metadata": {}, "outputs": [ { @@ -552,7 +553,7 @@ " [1. , 0.125 , 0.015625]])" ] }, - "execution_count": 10, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -562,62 +563,38 @@ ] }, { - "cell_type": "markdown", - "id": "c6f54a8f", + "cell_type": "code", + "execution_count": 12, + "id": "64fe7db9", "metadata": {}, - "source": [ - "While the vector of constraints $b$ has the following elements:\n", - "$$ b_{0} = 1 $$\n", - "$$ b_1 = b_2 = 0 $$\n", - "\n", - "Thus,\n", - "\n", - "$$\n", - "b =\n", - "\\begin{bmatrix}\n", - "1 \\\\\n", - "0 \\\\\n", - "0\n", - "\\end{bmatrix}\n", - "$$\n", - "\n", - "\n", - "And similarly in `lse`:" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "id": "c5b4f202", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([1., 0., 0.])" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [ + { + "data": { + "text/plain": [ + "array([1., 0., 0.])" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "lse.b" ] }, { "cell_type": "markdown", - "id": "fc47c29f", + "id": "16ab9f79", "metadata": {}, "source": [ - "The `lse` object has methods for finding the static coefficients $x_j$ satisfying the system of equations." + "With the LSE in hand, we solve for the static coefficients $x_j$ via `lse.solve()` (this is the direct $x = A^{-1}b$ solution)." ] }, { "cell_type": "code", - "execution_count": 12, - "id": "4459394c", + "execution_count": 13, + "id": "7b69192a", "metadata": {}, "outputs": [ { @@ -637,20 +614,18 @@ }, { "cell_type": "markdown", - "id": "6b228e6f", + "id": "1c238a59", "metadata": {}, "source": [ "##### Optimize for $x$ using an exact model\n", "\n", - "Alternatively to computing $x=A^{-1}b$, you can also use [setup_exact_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) to construct a [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance that uses the LSE as constraints and whose optimal solution will yield $x$.\n", - "\n", - "In the next section, it will be clear why this interface exists." + "Alternatively to computing $x = A^{-1}b$, you can use [setup\\_exact\\_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_exact_model.html) to construct a [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance that uses the LSE as constraints and whose optimal solution will yield $x$." ] }, { "cell_type": "code", - "execution_count": 13, - "id": "238ad3fe", + "execution_count": 14, + "id": "993465e9", "metadata": {}, "outputs": [ { @@ -662,25 +637,15 @@ } ], "source": [ - "from qiskit_addon_mpf.costs import setup_exact_problem\n", - "\n", "model_exact, coeffs_exact = setup_exact_problem(lse)\n", "model_exact.solve()\n", "print(coeffs_exact.value)" ] }, - { - "cell_type": "markdown", - "id": "31799d19", - "metadata": {}, - "source": [ - "As an indicator whether an MPF constructed with these coefficients will yield good results, we can use the L1-norm (see also Ref. [\\[1\\]](#references))." - ] - }, { "cell_type": "code", - "execution_count": 14, - "id": "ee706bcd", + "execution_count": 15, + "id": "eb61ea70", "metadata": {}, "outputs": [ { @@ -695,26 +660,23 @@ "print(\n", " \"L1 norm of the exact coefficients:\",\n", " np.linalg.norm(coeffs_exact.value, ord=1),\n", - ") # ord specifies the norm. ord=1 is for L1" + ")" ] }, { "cell_type": "markdown", - "id": "577dc339", + "id": "dac0472b", "metadata": {}, "source": [ "##### Optimize for $x$ using an approximate model\n", "\n", - "It might happen that the L1 norm for the chosen set of $k_j$ values is deemed too high.\n", - "If that is the case and you cannot choose a different set of $k_j$ values, you can use an approximate solution to the LSE instead of an exact one.\n", - "\n", - "To do so, simply use [setup_approximate_model](https://qiskit.github.io/qiskit-addon-mpf/stubs/qiskit_addon_mpf.static.setup_approximate_model.html) to construct a different [cvxpy.Problem](https://www.cvxpy.org/api_reference/cvxpy.problems.html#cvxpy.Problem) instance, which constrains the L1-norm to a chosen threshold while minimizing the difference of $Ax$ and $b$." + "It might happen that the $L_1$ norm for the chosen set of $k_j$ values is deemed too high. If that is the case and you cannot choose a different set of $k_j$ values, you can use an approximate solution that constrains the $L_1$-norm to a chosen threshold while minimizing $\\|Ax - b\\|$. Check out the guide on [How to use the approximate model](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html)." ] }, { "cell_type": "code", "execution_count": 16, - "id": "3e06be50", + "id": "0cd7dea4", "metadata": {}, "outputs": [ { @@ -727,8 +689,6 @@ } ], "source": [ - "from qiskit_addon_mpf.costs import setup_sum_of_squares_problem\n", - "\n", "model_approx, coeffs_approx = setup_sum_of_squares_problem(\n", " lse, max_l1_norm=1.5\n", ")\n", @@ -742,76 +702,29 @@ }, { "cell_type": "markdown", - "id": "214239bb", - "metadata": {}, - "source": [ - "Note that you have complete freedom over how to solve this optimization problem, which means that you can change the optimization solver, its convergence thresholds, and so on.\n", - "Check out the respective guide on [How to use the approximate model](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html)." - ] - }, - { - "cell_type": "markdown", - "id": "a3fca123", + "id": "10fd05eb", "metadata": {}, "source": [ "#### Dynamic MPF coefficients\n", "\n", - "In the previous section, we introduced a static MPF that improves on the standard Trotter approximation. However, this static version does not necessarily minimize the approximation error. Concretely, the static MPF, denoted $\\mu^S(t)$, is not the optimal projection of $\\rho(t)$ onto the subspace spanned by the product-formula states $\\{\\rho_{k_i}(t)\\}_{i=1}^r$.\n", - "\n", - "To address this, we consider a dynamic MPF (introduced in Ref. [\\[2\\]](#references) and experimentally demonstrated in Ref. [\\[3\\]](#references)) that does minimize the approximation error in the Frobenius norm. Formally, we focus on minimizing\n", - "\n", - "$$\n", - "\\|\\rho(t) - \\mu^D(t)\\|_F^2 \\;=\\; \\mathrm{Tr}\\bigl[ \\left( \\rho(t) - \\mu^D(t)\\right)^2 \\bigr],\n", - "$$\n", - "\n", - "with respect to some coefficients $x_i(t)$ at each time $t$. The *optimal* projector in the Frobenius norm is then $\\mu^D(t) \\;=\\; \\sum_{i=1}^r x_i(t)\\,\\rho_{k_i}(t)$, and we call $\\mu^D(t)$ the *dynamic* MPF. Plugging in the definitions above:\n", - "\n", + "The static MPF cancels Trotter-error terms in a Hamiltonian- and state-agnostic way, so it does not necessarily produce the smallest possible approximation error for a given Hamiltonian and initial state. The dynamic MPF (Refs. [\\[2\\]](#references), [\\[3\\]](#references)) instead finds time-dependent coefficients $x_i(t)$ that minimize the Frobenius-norm distance $\\|\\rho(t) - \\mu^D(t)\\|_F^2$ at each time $t$. As shown in the Background, this requires the overlap matrix $M_{ij}(t)$ between Trotter-evolved states and the overlap $L_i(t)$ with the exact state — both of which we estimate using tensor-network (TeNPy) backends in `qiskit_addon_mpf`.\n", "\n", + "To set up the dynamic LSE we need three ingredients:\n", "\n", - "$$\n", - "\\|\\rho(t) - \\mu^D(t)\\|_F^2\n", - "\\;=\\; \\\\\n", - "= \\mathrm{Tr}\\bigl[ \\left( \\rho(t) - \\mu^D(t)\\right)^2 \\bigr]\n", - "\\;=\\; \\\\\n", - "= \\mathrm{Tr}\\bigl[ \\left( \\rho(t) - \\sum_{i=1}^r x_i(t)\\,\\rho_{k_i}(t) \\right) \\left( \\rho(t) - \\sum_{j=1}^r x_j(t)\\,\\rho_{k_j}(t) \\right) \\bigr]\n", - "\\;=\\; \\\\\n", - "= 1 \\;+\\; \\sum_{i,j=1}^r M_{i,j}(t)\\,x_i(t)\\,x_j(t)\n", - "\\;-\\;\n", - "2 \\sum_{i=1}^r L_i^{\\mathrm{exact}}(t)\\,x_i(t),\n", - "$$\n", - "\n", - "where $M_{i,j}(t)$ is the *Gram matrix*, defined by\n", - "\n", - "$$\n", - "M_{i,j}(t) \\;=\\; \\mathrm{Tr}\\bigl[\\rho_{k_i}(t)\\,\\rho_{k_j}(t)\\bigr]\n", - "\\;=\\;\n", - "\\bigl|\\langle \\psi_{\\mathrm{in}} \\!\\mid S\\bigl(t/k_i\\bigr)^{-k_i}\\,S\\bigl(t/k_j\\bigr)^{k_j} \\!\\mid \\psi_{\\mathrm{in}} \\rangle \\bigr|^2.\n", - "$$\n", - "\n", + "1. An **approximate evolver factory** that the addon will run for each $k_j$ to produce $\\rho_{k_j}(t)$ as an MPS/MPO. We build it from the layered structure of the order-$2$ Trotter circuit (one layer per `slice_by_depth`), wrapped as a `LayerwiseEvolver` with TeNPy truncation parameters.\n", + "2. An **exact evolver factory** that produces a high-accuracy reference $\\rho(t)$. We use a small-time-step fourth-order Suzuki-Trotter circuit (`dt=0.1`, `order=4`) as a proxy for exact evolution.\n", + "3. An **identity factory** and an **initial-state MPS** that seed the TeNPy simulation.\n", "\n", - "and\n", - "\n", - "$$\n", - "L_i^{\\mathrm{exact}}(t) = \\mathrm{Tr}[\\rho(t)\\,\\rho_{k_i}(t)]\n", - "$$\n", - "\n", - "represents the overlap between the exact state $\\rho(t)$ and each product-formula approximation $\\rho_{k_i}(t)$. In practical scenarios, these overlaps may only be measured approximately due to noise or partial access to $\\rho(t)$.\n", - "\n", - "Here, $\\lvert\\psi_{\\mathrm{in}}\\rangle$ is the initial state, and $S(\\cdot)$ is the operation applied in the product formula. By choosing the coefficients $x_i(t)$ that minimize this expression (and handling approximate overlap data when $\\rho(t)$ is not fully known), we obtain the “best” (in a Frobenius-norm sense) dynamic approximation of $\\rho(t)$ within the MPF subspace. The quantities $L_i(t)$ and $M_{i,j}(t)$ can be calculated efficiently using tensor network methods [\\[3\\]](#references). The MPF Qiskit addon provides several \"backends\" to carry out the calculation. The example below shows the most flexible way to do so, and the [TeNPy layer-based backend](https://qiskit.github.io/qiskit-addon-mpf/apidocs/qiskit_addon_mpf.backends.tenpy_layers.html#module-qiskit_addon_mpf.backends.tenpy_layers) docs also explain in great detail. To use this method, start from the circuit implementing the desired time-evolution and create models that represent these operations from the layers of the corresponding circuit. Finally, an `Evolver` object is created that can be used to generate the time-evolved quantities $M_{i,j}(t)$ and $L_i(t)$. We start by creating the `Evolver` object corresponding to the approximate time-evolution ([`ApproxEvolverFactory`](https://qiskit.github.io/qiskit-addon-mpf/apidocs/qiskit_addon_mpf.dynamic.html#qiskit_addon_mpf.dynamic.ApproxEvolverFactory)) implemented by the circuits. In particular, pay extra attention to the `order` variable so that they match. Note that in generating the circuits corresponding to the approximate time-evolution, we use placeholder values for the `time = 1.0` and the number of Trotter steps (`reps=1`). The correct approximating circuits are then produced by the dynamic problem solver in `setup_dynamic_lse`." + "The cell below constructs the approximate evolver factory." ] }, { "cell_type": "code", "execution_count": 17, - "id": "87323ed2", + "id": "52d78403", "metadata": {}, "outputs": [], "source": [ - "from qiskit_addon_utils.slicing import slice_by_depth\n", - "from qiskit_addon_mpf.backends.tenpy_layers import LayerModel\n", - "from qiskit_addon_mpf.backends.tenpy_layers import LayerwiseEvolver\n", - "from functools import partial\n", - "\n", "# Create approximate time-evolution circuits\n", "single_2nd_order_circ = generate_time_evolution_circuit(\n", " hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)\n", @@ -844,26 +757,26 @@ }, { "cell_type": "markdown", - "id": "fd39f15a", + "id": "91d2783e", "metadata": {}, "source": [ "\n", - "The options of `LayerwiseEvolver` that determine the details of the tensor network simulation must be chosen carefully to avoid setting up an ill-defined optimization problem.\n", + " The options of `LayerwiseEvolver` that determine the details of the tensor network simulation must be chosen carefully to avoid setting up an ill-defined optimization problem.\n", "" ] }, { "cell_type": "markdown", - "id": "d6b611e3", + "id": "1c702f67", "metadata": {}, "source": [ - "We then set up the exact evolver (for example, [`ExactEvolverFactory`](https://qiskit.github.io/qiskit-addon-mpf/apidocs/qiskit_addon_mpf.dynamic.html#qiskit_addon_mpf.dynamic.ExactEvolverFactory)), which returns an [`Evolver`](https://qiskit.github.io/qiskit-addon-mpf/apidocs/qiskit_addon_mpf.backends.html#qiskit_addon_mpf.backends.Evolver) object computing the true or “reference” time-evolution. Realistically, we would approximate the exact evolution by using a higher-order Suzuki–Trotter formula or another reliable method with a small time step. Below, we approximate the exact time-evolved state with a fourth-order Suzuki-Trotter formula using a small time step `dt=0.1`, which means the number of Trotter steps used at time $t$ is $k=t/dt$. We also specify some TeNPy-specific truncation options to bound the maximum bond dimension of the underlying tensor network, as well as the minimum singular values of the split tensor network bonds. These parameters can affect the accuracy of the expectation value calculated with the dynamic MPF coefficients, so it is important to explore a range of values to find the optimal balance between computational time and accuracy. Note that the calculation of the MPF coefficients does not rely of the expectation value of the PF obtained from hardware execution, and therefore it can be tuned in post-processing." + "We approximate the exact time-evolved state with a fourth-order Suzuki-Trotter formula using a small time step `dt=0.1`. The TeNPy truncation parameters can affect accuracy, so it is important to explore a range of values." ] }, { "cell_type": "code", - "execution_count": 20, - "id": "10ca6c85", + "execution_count": 18, + "id": "abab8bfc", "metadata": {}, "outputs": [], "source": [ @@ -894,23 +807,19 @@ }, { "cell_type": "markdown", - "id": "1eb4d98c", + "id": "2486594c", "metadata": {}, "source": [ - "Next, create your system’s initial state in a format compatible with TeNPy (for example, an `MPS_neel_state`=$\\vert 0101...01 \\rangle$). This sets up the many-body wavefunction you’ll evolve in time $\\lvert\\psi_{\\mathrm{in}}\\rangle$ as a tensor." + "Finally, we define an `identity_factory` that yields the initial MPO state and prepare the Néel initial state as an MPS that matches the lattice used by the layered Trotter model." ] }, { "cell_type": "code", - "execution_count": null, - "id": "69758099", + "execution_count": 19, + "id": "1e216575", "metadata": {}, "outputs": [], "source": [ - "from qiskit_addon_mpf.backends.tenpy_tebd import MPOState\n", - "from qiskit_addon_mpf.backends.tenpy_tebd import MPS_neel_state\n", - "\n", - "\n", "def identity_factory():\n", " return MPOState.initialize_from_lattice(models[0].lat, conserve=True)\n", "\n", @@ -920,16 +829,16 @@ }, { "cell_type": "markdown", - "id": "057d479e", + "id": "5658315c", "metadata": {}, "source": [ - "For each time step $t$ we set up the dynamic linear system of equations with the [`setup_dynamic_lse`](https://qiskit.github.io/qiskit-addon-mpf/apidocs/qiskit_addon_mpf.dynamic.html) method. The corresponding object contains the information about the dynamic MPF problem: `lse.A` gives the Gram matrix $M$ while `lse.b` gives the overlap $L$. We can then solve the LSE (when not ill-defined) to find the dynamic coefficients using the `setup_frobenius_problem`. It is important to note the difference with the static coefficients, which only depend on the details of the product formula used and are independent of the details of the time-evolution (Hamiltonian and initial state)." + "With the factories in place we now compute the dynamic coefficients at each evolution time. For each $t$, `setup_dynamic_lse` builds the relevant overlap matrices via TeNPy, and `setup_frobenius_problem` returns a `cvxpy.Problem` that minimizes the Frobenius-norm cost. The solver returns coefficients $x_j(t)$ tailored to that time; we collect them in `mpf_dynamic_coeffs_list`. If the solver fails for a given $t$, we fall back to zero coefficients so the loop continues." ] }, { "cell_type": "code", - "execution_count": 22, - "id": "bed976ac", + "execution_count": 20, + "id": "b05dc012", "metadata": {}, "outputs": [ { @@ -962,9 +871,6 @@ } ], "source": [ - "from qiskit_addon_mpf.dynamic import setup_dynamic_lse\n", - "from qiskit_addon_mpf.costs import setup_frobenius_problem\n", - "\n", "mpf_dynamic_coeffs_list = []\n", "for t in trotter_times:\n", " print(f\"Computing dynamic coefficients for time={t}\")\n", @@ -988,22 +894,24 @@ }, { "cell_type": "markdown", - "id": "bf4a82e0", + "id": "4f8814e3", "metadata": {}, "source": [ - "Finally, plot these expectation values throughout the evolution time." + "#### Combine Trotter expectation values with the MPF coefficients\n", + "\n", + "Now we evaluate $\\langle A \\rangle_{\\text{MPF}}(t) = \\sum_j x_j \\, \\langle A \\rangle_{k_j}(t)$ for each set of coefficients (static-exact, static-approximate, and dynamic), propagate the per-circuit standard errors, and plot the resulting time series against the exact-diagonalization curve." ] }, { "cell_type": "code", - "execution_count": null, - "id": "2da9c948", + "execution_count": 21, + "id": "35042576", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "\"Output" + "\"Output" ] }, "metadata": {}, @@ -1011,8 +919,6 @@ } ], "source": [ - "import matplotlib.pyplot as plt\n", - "\n", "sym = {1: \"^\", 2: \"s\", 4: \"p\"}\n", "# Get expectation values at all times for each Trotter step\n", "for k, step in enumerate(mpf_trotter_steps):\n", @@ -1032,7 +938,7 @@ " marker=sym[step],\n", " color=\"grey\",\n", " label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n", - " ) # , , )\n", + " )\n", "\n", "# Get expectation values at all times for the static MPF with exact coeffs\n", "exact_mpf_curve, exact_mpf_curve_error = [], []\n", @@ -1061,7 +967,7 @@ ")\n", "\n", "\n", - "# Get expectation values at all times for the static MPF with approximate\n", + "# Get expectation values at all times for the static MPF with approximate coeffs\n", "approx_mpf_curve, approx_mpf_curve_error = [], []\n", "for trotter_expvals, trotter_stds in zip(\n", " mpf_expvals_all_times, mpf_stds_all_times\n", @@ -1083,11 +989,12 @@ " yerr=approx_mpf_curve_error,\n", " markersize=4,\n", " marker=\"o\",\n", + " label=\"Static MPF - Approx\",\n", " color=\"orange\",\n", - " label=\"Static MPF - Approximate\",\n", ")\n", "\n", - "# # Get expectation values at all times for the dynamic MPF\n", + "\n", + "# Get expectation values at all times for the dynamic MPF\n", "dynamic_mpf_curve, dynamic_mpf_curve_error = [], []\n", "for trotter_expvals, trotter_stds, dynamic_coeffs in zip(\n", " mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list\n", @@ -1109,594 +1016,22 @@ " yerr=dynamic_mpf_curve_error,\n", " markersize=4,\n", " marker=\"o\",\n", - " color=\"pink\",\n", " label=\"Dynamic MPF\",\n", + " color=\"pink\",\n", ")\n", "\n", "\n", + "# Exact expectation values\n", "plt.plot(\n", " exact_evolution_times,\n", " exact_expvals,\n", - " lw=3,\n", " color=\"red\",\n", - " label=\"Exact time-evolution\",\n", - ")\n", - "\n", - "\n", - "plt.title(\n", - " f\"Expectation values for (ZZ,{(L//2-1, L//2)}) as a function of time\"\n", - ")\n", - "plt.xlabel(\"Time\")\n", - "plt.ylabel(\"Expectation Value\")\n", - "plt.legend()\n", - "plt.grid()" - ] - }, - { - "cell_type": "markdown", - "id": "3a4c6912", - "metadata": {}, - "source": [ - "In the cases like the example above, where the $k=1$ PF behaves poorly at all times, the quality of the dynamic MPF results is also heavily affected. In such situations, it is useful to investigate the possibility of using individual PFs with higher number of Trotter steps to improve the overall quality of the results. In these simulations, we see the interplay of different types of errors: error from finite sampling, and Trotter error from the product formulas. MPF helps to reduce the Trotter error due to the product formulas but incur in higher sampling error compared to the product formulas. This can be advantageous, as product formulas can reduce the sampling error with increased sampling, but the systematic error due to the Trotter approximation remains untouched.\n", - "\n", - "Another interesting behavior that we can observe from the plot is that the expectation value for the PF for $k=1$ starts to behave erratically (on top of not being a good approximation for the exact one) at times for which $t/k > 1 $, as explained in the [guide](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) on how to choose the number of Trotter steps." - ] - }, - { - "cell_type": "markdown", - "id": "2a16f346", - "metadata": {}, - "source": [ - "### Step 1: Map classical inputs to a quantum problem\n", - "Let's now consider a single time $t=1.0$ and calculate the expectation value of the magnetization with the various methods using one QPU. The particular choice of $t$ was done so to maximize the difference between the various methods and observe their relative efficacy. To determine the window of time for which dynamic MPF is guaranteed to produce observables with lower error than any of the individual Trotter formulas within the multi-product, we can implement the “MPF test” - see equation (17) and surrounding text in [\\[3\\]](#references)." - ] - }, - { - "cell_type": "markdown", - "id": "623cc755", - "metadata": {}, - "source": [ - "#### Set up the Trotter circuits\n", - "\n", - "At this point, we have found our expansion coefficients, $x$, and all that is left to do is to generate the Trotterized quantum circuits.\n", - "Once again, the [qiskit_addon_utils.problem_generators](/docs/api/qiskit-addon-utils/problem-generators) module comes to the rescue with a useful function to do this:" - ] - }, - { - "cell_type": "code", - "execution_count": 26, - "id": "3dc0dcf2", - "metadata": {}, - "outputs": [], - "source": [ - "from qiskit.synthesis import SuzukiTrotter\n", - "from qiskit_addon_utils.problem_generators import (\n", - " generate_time_evolution_circuit,\n", - ")\n", - "from qiskit import QuantumCircuit\n", - "\n", - "\n", - "total_time = 1.0\n", - "mpf_circuits = []\n", - "for k in mpf_trotter_steps:\n", - " # Initial Neel state preparation\n", - " circuit = QuantumCircuit(L)\n", - " circuit.x([i for i in range(L) if i % 2 != 0])\n", - "\n", - " trotter_circ = generate_time_evolution_circuit(\n", - " hamiltonian,\n", - " synthesis=SuzukiTrotter(order=order, reps=k),\n", - " time=total_time,\n", - " )\n", - "\n", - " circuit.compose(trotter_circ, inplace=True)\n", - "\n", - " mpf_circuits.append(pm.run(circuit))" - ] - }, - { - "cell_type": "code", - "execution_count": 27, - "id": "87d2ac0c", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "\"Output" - ] - }, - "execution_count": 27, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "mpf_circuits[-1].draw(\"mpl\", fold=-1, scale=0.4)" - ] - }, - { - "cell_type": "markdown", - "id": "97eb3019", - "metadata": {}, - "source": [ - "### Step 2: Optimize problem for quantum hardware execution\n", - "Let's return to the calculation of the expectation value for a single time point. We'll pick a backend for executing the experiment on hardware." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c9a3020e", - "metadata": {}, - "outputs": [], - "source": [ - "from qiskit_ibm_runtime import QiskitRuntimeService\n", - "\n", - "\n", - "service = QiskitRuntimeService()\n", - "backend = service.least_busy(min_num_qubits=127)\n", - "print(backend)\n", - "\n", - "qubits = list(range(backend.num_qubits))" - ] - }, - { - "cell_type": "markdown", - "id": "6eeaac1b", - "metadata": {}, - "source": [ - "Then we remove outlier qubits from the coupling map to ensure that the layout stage of the transpiler doesn't include them. Below we use the reported backend properties stored in the `target` object and remove qubits that either have measurement error or a two-qubit gate above a certain threshold (`max_meas_err`, `max_twoq_err`) or $T_2$ time (which determines the loss of coherence) below a certain threshold (`min_t2`)." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8bb2c619", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "size of largest component 10\n" - ] - } - ], - "source": [ - "import copy\n", - "from qiskit.transpiler import Target, CouplingMap\n", - "\n", - "target = backend.target\n", - "instruction_2q = \"cz\"\n", - "\n", - "cmap = target.build_coupling_map(filter_idle_qubits=True)\n", - "cmap_list = list(cmap.get_edges())\n", - "\n", - "max_meas_err = 0.012\n", - "min_t2 = 40\n", - "max_twoq_err = 0.005\n", - "\n", - "# Remove qubits with bad measurement or t2\n", - "cust_cmap_list = copy.deepcopy(cmap_list)\n", - "for q in range(target.num_qubits):\n", - " meas_err = target[\"measure\"][(q,)].error\n", - " if target.qubit_properties[q].t2 is not None:\n", - " t2 = target.qubit_properties[q].t2 * 1e6\n", - " else:\n", - " t2 = 0\n", - " if meas_err > max_meas_err or t2 < min_t2:\n", - " # print(q)\n", - " for q_pair in cmap_list:\n", - " if q in q_pair:\n", - " try:\n", - " cust_cmap_list.remove(q_pair)\n", - " except ValueError:\n", - " continue\n", - "\n", - "# Remove qubits with bad 2q gate or t2\n", - "for q in cmap_list:\n", - " twoq_gate_err = target[instruction_2q][q].error\n", - " if twoq_gate_err > max_twoq_err:\n", - " # print(q)\n", - " for q_pair in cmap_list:\n", - " if q == q_pair:\n", - " try:\n", - " cust_cmap_list.remove(q_pair)\n", - " except ValueError:\n", - " continue\n", - "\n", - "\n", - "cust_cmap = CouplingMap(cust_cmap_list)\n", - "\n", - "cust_target = Target.from_configuration(\n", - " basis_gates=backend.configuration().basis_gates\n", - " + [\"measure\"], # or whatever new set of gates\n", - " coupling_map=cust_cmap,\n", - ")\n", - "\n", - "sorted_components = sorted(\n", - " [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],\n", - " reverse=True,\n", - ")\n", - "print(\"size of largest component\", len(sorted_components[0]))" - ] - }, - { - "cell_type": "markdown", - "id": "2aba48dd", - "metadata": {}, - "source": [ - "We want to set `max_meas_err`, `min_t2`, and `max_twoq_err` such that we find a large enough subset of qubits that supports the circuit to run. In our case it's enough to find a 10-qubit 1D chain." - ] - }, - { - "cell_type": "code", - "execution_count": 34, - "id": "c5d8e90b", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "\"Output" - ] - }, - "execution_count": 34, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "cust_cmap.draw()" - ] - }, - { - "cell_type": "markdown", - "id": "4d11a75a", - "metadata": {}, - "source": [ - "We can then map the circuit and observable onto physical qubits of the device." - ] - }, - { - "cell_type": "code", - "execution_count": 72, - "id": "bbbec0ec", - "metadata": {}, - "outputs": [], - "source": [ - "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", - "\n", - "transpiler = generate_preset_pass_manager(\n", - " optimization_level=3, target=cust_target\n", - ")\n", - "\n", - "transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]\n", - "\n", - "qubits_layouts = [\n", - " [\n", - " idx\n", - " for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()\n", - " if qb._register.name != \"ancilla\"\n", - " ]\n", - " for circuit in transpiled_circuits\n", - "]\n", - "\n", - "transpiled_circuits = []\n", - "for circuit, layout in zip(mpf_circuits, qubits_layouts):\n", - " transpiler = generate_preset_pass_manager(\n", - " optimization_level=3, backend=backend, initial_layout=layout\n", - " )\n", - " transpiled_circuit = transpiler.run(circuit)\n", - " transpiled_circuits.append(transpiled_circuit)\n", - "\n", - "\n", - "# transform the observable defined on virtual qubits to\n", - "# an observable defined on all physical qubits\n", - "isa_observables = [\n", - " observable.apply_layout(circ.layout) for circ in transpiled_circuits\n", - "]" - ] - }, - { - "cell_type": "code", - "execution_count": 73, - "id": "25ce07a6", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "51\n", - "OrderedDict([('sx', 310), ('rz', 232), ('cz', 132), ('x', 19)])\n" - ] - }, - { - "data": { - "text/plain": [ - "\"Output" - ] - }, - "execution_count": 73, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "print(transpiled_circuits[-1].depth(lambda x: x.operation.num_qubits == 2))\n", - "print(transpiled_circuits[-1].count_ops())\n", - "transpiled_circuits[-1].draw(\"mpl\", idle_wires=False, fold=False)" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "id": "c5edec73", - "metadata": {}, - "source": [ - "### Step 3: Execute using Qiskit primitives\n", - "With the Estimator primitive we can obtain the estimation of expectation value from the QPU. We execute the optimized AQC circuits with additional error mitigation and suppression techniques." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f91a669f", - "metadata": {}, - "outputs": [], - "source": [ - "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n", - "\n", - "\n", - "estimator = Estimator(mode=backend)\n", - "estimator.options.default_shots = 30000\n", - "\n", - "# Set simple error suppression/mitigation options\n", - "estimator.options.dynamical_decoupling.enable = True\n", - "estimator.options.twirling.enable_gates = True\n", - "estimator.options.twirling.enable_measure = True\n", - "estimator.options.twirling.num_randomizations = \"auto\"\n", - "estimator.options.twirling.strategy = \"active-accum\"\n", - "estimator.options.resilience.measure_mitigation = True\n", - "estimator.options.experimental.execution_path = \"gen3-turbo\"\n", - "\n", - "estimator.options.resilience.zne_mitigation = True\n", - "estimator.options.resilience.zne.noise_factors = (1, 3, 5)\n", - "estimator.options.resilience.zne.extrapolator = (\"exponential\", \"linear\")\n", - "\n", - "estimator.options.environment.job_tags = [\"mpf small\"]\n", - "\n", - "\n", - "job = estimator.run(\n", - " [\n", - " (circ, observable)\n", - " for circ, observable in zip(transpiled_circuits, isa_observables)\n", - " ]\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "87484abf", - "metadata": {}, - "source": [ - "### Step 4: Post-process and return result in desired classical format\n", - "\n", - "\n", - "The only post-processing step is to combine the expectation value obtained from the Qiskit Runtime primitives at different Trotter steps using the respective MPF coefficients. For an observable $A$ we have:\n", - "\n", - "$$ \\langle A \\rangle_{\\text{mpf}} = \\text{Tr} [A \\mu(t)] = \\sum_{j} x_j \\text{Tr} [A \\rho_{k_j}] = \\sum_{j} x_j \\langle A \\rangle_j$$\n", - "\n", - "First, we extract the individual expectation values obtained for each of the Trotter circuits:" - ] - }, - { - "cell_type": "code", - "execution_count": 39, - "id": "a5ef9913", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[array(-0.06361607), array(-0.23820448), array(-0.50271805)]\n" - ] - } - ], - "source": [ - "result_exp = job.result()\n", - "evs_exp = [res.data.evs for res in result_exp]\n", - "evs_std = [res.data.stds for res in result_exp]\n", - "\n", - "print(evs_exp)" - ] - }, - { - "cell_type": "markdown", - "id": "446ab041", - "metadata": {}, - "source": [ - "Next, we simply recombine them with our MPF coefficients to yield the total expectation values of the MPF. Below, we do so for each of the different ways by which we have computed $x$." - ] - }, - { - "cell_type": "code", - "execution_count": 29, - "id": "3ee35ecb", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Exact static MPF expectation value: -0.6329590442738475 +- 0.012798249760406036\n", - "Approximate static MPF expectation value: -0.5690390035339492 +- 0.010459559917168473\n", - "Dynamic MPF expectation value: -0.4655579758795695 +- 0.007639139186720507\n" - ] - } - ], - "source": [ - "exact_mpf_std = np.sqrt(\n", - " sum(\n", - " [\n", - " (coeff**2) * (std**2)\n", - " for coeff, std in zip(coeffs_exact.value, evs_std)\n", - " ]\n", - " )\n", - ")\n", - "print(\n", - " \"Exact static MPF expectation value: \",\n", - " evs_exp @ coeffs_exact.value,\n", - " \"+-\",\n", - " exact_mpf_std,\n", - ")\n", - "approx_mpf_std = np.sqrt(\n", - " sum(\n", - " [\n", - " (coeff**2) * (std**2)\n", - " for coeff, std in zip(coeffs_approx.value, evs_std)\n", - " ]\n", - " )\n", - ")\n", - "print(\n", - " \"Approximate static MPF expectation value: \",\n", - " evs_exp @ coeffs_approx.value,\n", - " \"+-\",\n", - " approx_mpf_std,\n", - ")\n", - "dynamic_mpf_std = np.sqrt(\n", - " sum(\n", - " [\n", - " (coeff**2) * (std**2)\n", - " for coeff, std in zip(mpf_dynamic_coeffs_list[7], evs_std)\n", - " ]\n", - " )\n", - ")\n", - "print(\n", - " \"Dynamic MPF expectation value: \",\n", - " evs_exp @ mpf_dynamic_coeffs_list[7],\n", - " \"+-\",\n", - " dynamic_mpf_std,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "8a99f6b2", - "metadata": {}, - "source": [ - "Finally, for this small problem we can compute the exact reference value using [scipy.linalg.expm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.expm.html) as follows:" - ] - }, - { - "cell_type": "code", - "execution_count": 30, - "id": "4a7dd8fc", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Exact expectation value -0.39909900734489434\n" - ] - } - ], - "source": [ - "from scipy.linalg import expm\n", - "from qiskit.quantum_info import Statevector\n", - "\n", - "exp_H = expm(-1j * total_time * hamiltonian.to_matrix())\n", - "\n", - "initial_state_circuit = QuantumCircuit(L)\n", - "initial_state_circuit.x([i for i in range(L) if i % 2 != 0])\n", - "initial_state = Statevector(initial_state_circuit).data\n", - "\n", - "time_evolved_state = exp_H @ initial_state\n", - "\n", - "exact_obs = (\n", - " time_evolved_state.conj() @ observable.to_matrix() @ time_evolved_state\n", - ")\n", - "print(\"Exact expectation value \", exact_obs.real)" - ] - }, - { - "cell_type": "code", - "execution_count": 31, - "id": "a3eefe73", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "\"Output" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "sym = {1: \"^\", 2: \"s\", 4: \"p\"}\n", - "# Get expectation values at all times for each Trotter step\n", - "for k, step in enumerate(mpf_trotter_steps):\n", - " plt.errorbar(\n", - " k,\n", - " evs_exp[k],\n", - " yerr=evs_std[k],\n", - " alpha=0.5,\n", - " markersize=4,\n", - " marker=sym[step],\n", - " color=\"grey\",\n", - " label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n", - " ) # , , )\n", - "\n", - "\n", - "plt.errorbar(\n", - " 3,\n", - " evs_exp @ coeffs_exact.value,\n", - " yerr=exact_mpf_std,\n", - " markersize=4,\n", - " marker=\"o\",\n", - " color=\"purple\",\n", - " label=\"Static MPF\",\n", - ")\n", - "\n", - "plt.errorbar(\n", - " 4,\n", - " evs_exp @ coeffs_approx.value,\n", - " yerr=approx_mpf_std,\n", - " markersize=4,\n", - " marker=\"o\",\n", - " color=\"orange\",\n", - " label=\"Approximate static MPF\",\n", - ")\n", - "\n", - "plt.errorbar(\n", - " 5,\n", - " evs_exp @ mpf_dynamic_coeffs_list[7],\n", - " yerr=dynamic_mpf_std,\n", - " markersize=4,\n", - " marker=\"o\",\n", - " color=\"pink\",\n", - " label=\"Dynamic MPF\",\n", - ")\n", - "\n", - "plt.axhline(\n", - " y=exact_obs.real,\n", " linestyle=\"--\",\n", - " color=\"red\",\n", " label=\"Exact time-evolution\",\n", ")\n", "\n", - "\n", - "plt.title(\n", - " f\"Expectation values for (ZZ,{(L//2-1, L//2)}) at time {total_time} for the different methods \"\n", - ")\n", - "plt.xlabel(\"Method\")\n", + "plt.title(f\"$\\\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\\\rangle$ vs time\")\n", + "plt.xlabel(\"Time\")\n", "plt.ylabel(\"Expectation Value\")\n", "plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2), ncol=2)\n", "plt.grid(alpha=0.1)\n", @@ -1706,164 +1041,69 @@ }, { "cell_type": "markdown", - "id": "921fb0c1", + "id": "34923748", "metadata": {}, "source": [ - "In the example above, the dynamic MPF method performs best in terms of expectation value, improving over what we would have obtained by using the highest number of Trotter steps alone. Even though the various MPF techniques are not always achieving an improved expectation value compared to the highest number of Trotter steps (like the exact and the approximate model in the plot above), the standard deviation of these values captures well the increased variance incurred when using the MPF technique. This highlights uncertainty around the obtained expectation value, which always includes the expectation value we would expect from an exact time-evolution of the system. On the other hand, the expectation values calculated with the lower number of Trotter steps fails to capture the exact expectation value within their uncertainty, thus confidently returning the wrong result." - ] - }, - { - "cell_type": "code", - "execution_count": 32, - "id": "7d639f4a", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "relative error for each trotter steps [0.33548293650112293, 0.16089452939226306, 0.10361904247828346]\n", - "relative error with MPF exact coeffs 0.2338600369291003\n", - "relative error with MPF approx coeffs 0.16993999618905486\n", - "relative error with MPF dynamic coeffs 0.06645896853467514\n" - ] - } - ], - "source": [ - "def relative_error(ev, exact_ev):\n", - " return abs(ev - exact_ev)\n", + "The plot above illustrates the interplay between Trotter error and sampling error.\n", "\n", + "- **Trotter error.** The individual product formulas (grey markers) deviate from the exact curve more and more as time grows. The $k=1$ circuit has the largest deviation and is the shallowest, but it is also already in the regime where $t/k \\gtrsim 1$, so the leading $1/k^{2}$ error term is large. The MPF combinations (colored markers) cancel several of these leading Trotter-error terms, so they track the exact curve much more closely than any single $k_j$ circuit. The remaining gap reflects the higher-order Trotter terms that the MPF does *not* cancel: an order-$2$, $r=3$ static MPF only kills the first two error orders, and at large $t/k_{\\min}$ the uncanceled tail eventually dominates — so MPF does not guarantee that very shallow circuits remain accurate at arbitrary times.\n", "\n", - "relative_error_k = [relative_error(ev, exact_obs.real) for ev in evs_exp]\n", - "relative_error_mpf = relative_error(evs_exp @ mpf_coeffs, exact_obs.real)\n", - "relative_error_approx_mpf = relative_error(\n", - " evs_exp @ coeffs_approx.value, exact_obs.real\n", - ")\n", - "relative_error_dynamic_mpf = relative_error(\n", - " evs_exp @ mpf_dynamic_coeffs_list[7], exact_obs.real\n", - ")\n", + "- **Sampling error.** The wider error bars on the MPF curves are a direct consequence of the linear combination: propagating independent per-circuit standard errors $\\sigma_{k_j}$ gives a total variance $\\sigma_{\\text{MPF}}^2 = \\sum_j x_j^2 \\, \\sigma_{k_j}^2$. Therefore, the larger the $\\|x\\|_2$ (and in practice the $\\|x\\|_1$, which is what we control), the more shots are required to reach a given target uncertainty. This is the trade-off behind the approximate-solver option in the Background: we cap $\\|x\\|_1$ to keep this overhead manageable. Crucially, unlike Trotter error, sampling error shrinks with $1/\\sqrt{N_{\\text{shots}}}$, so it can always be reduced by spending more shots.\n", "\n", - "print(\"relative error for each trotter steps\", relative_error_k)\n", - "print(\"relative error with MPF exact coeffs\", relative_error_mpf)\n", - "print(\"relative error with MPF approx coeffs\", relative_error_approx_mpf)\n", - "print(\"relative error with MPF dynamic coeffs\", relative_error_dynamic_mpf)" + "In the large-scale hardware example below, hardware noise enters as an additional error source on each $\\langle A \\rangle_{k_j}$, which similarly gets amplified by the MPF coefficients. We will see how error mitigation interacts with MPFs in that section." ] }, { "cell_type": "markdown", - "id": "0130299b", + "id": "6fa763ff", "metadata": {}, "source": [ - "## Part II: scale it up\n", + "## Large-scale hardware example\n", "\n", - "Let's scale the problem up beyond what is possible to simulate exactly. In this section we will focus on reproducing some of the results shown in Ref. [\\[3\\]](#references)." + "In this section we scale the problem up beyond what is possible to simulate exactly. We reproduce some of the results shown in Ref. [\\[3\\]](#references), using a 50-qubit XXZ chain at time $t = 3$. We use the same four-step workflow as above, but pull out the hardware-specific details (qubit filtering, transpilation, and error mitigation) into their own subsections after Step 1." ] }, { "cell_type": "markdown", - "id": "85a5e83a", + "id": "481fea70", "metadata": {}, "source": [ - "### Step 1: Map classical inputs to a quantum problem" - ] - }, - { - "cell_type": "markdown", - "id": "951205d9", - "metadata": {}, - "source": [ - "#### Hamiltonian\n", + "### Step 1: Map classical inputs to a quantum problem\n", "\n", - "For the large-scale example, we use the XXZ model on a line of 50 sites:\n", + "The mapping mirrors the small-scale example: define a Hamiltonian, choose Trotter parameters, compute MPF coefficients (static and dynamic), and build circuits. The key differences are:\n", "\n", - "$$\n", - "\\hat{\\mathcal{H}}_{XXZ} = \\sum_{i=1}^{L-1} J_{i,(i+1)}\\left(X_i X_{(i+1)}+Y_i Y_{(i+1)}+ 2\\cdot Z_i Z_{(i+1)} \\right) \\, ,\n", - "$$\n", + "- An **XXZ Hamiltonian** on 50 sites with random couplings drawn from $\\mathcal{U}(0.5, 1.5)$ (Ref. [\\[3\\]](#references)).\n", + "- A **symmetric** second-order Trotter formula with $k_j = [2, 3, 4]$ (so $\\chi=1$, `symmetric=True`).\n", + "- A single fixed evolution time $t = 3$.\n", + "- An additional **single-circuit comparison run with $k = 6$ Trotter steps**, used as a baseline. We chose $k = 6$ because its two-qubit depth on hardware is comparable to the deepest MPF constituent ($k_{\\max}=4$) plus the overhead of running multiple MPF circuits — so $k=6$ is a fair \"single deep circuit\" comparison against the MPF combination, not a circuit that targets the MPF's effective Trotter error (which would require many more steps).\n", "\n", - "where $J_{i,(i+1)}$ is a random coefficient corresponding to edge $(i, i+1)$. This is the Hamiltonian considered in the demonstration presented in Ref. [\\[3\\]](#references)." - ] - }, - { - "cell_type": "code", - "execution_count": 33, - "id": "34bf68ac", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "\"Output" - ] - }, - "execution_count": 33, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "L = 50\n", - "# Generate some coupling map to use for this example\n", - "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n", - "graphviz_draw(coupling_map.graph, method=\"circo\")" + "Note that even though we are still in Step 1 here (mapping and circuit construction), we also pre-compute the dynamic coefficients alongside the static ones in this cell. Dynamic coefficients depend on $H$ and $t$ but not on the quantum measurements, so they can be computed any time before Step 4. We do it now to keep all the MPF-specific setup in one place." ] }, { "cell_type": "code", - "execution_count": null, - "id": "06e49110", + "execution_count": 22, + "id": "a019ac32", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXX', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYY', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZ', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'XXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'YYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZI', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IXXIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IYYIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],\n", - " coeffs=[1. +0.j, 2.09762701+0.j, 2.09762701+0.j, 4.19525402+0.j,\n", - " 2.43037873+0.j, 2.43037873+0.j, 4.86075747+0.j, 2.20552675+0.j,\n", - " 2.20552675+0.j, 4.4110535 +0.j, 2.08976637+0.j, 2.08976637+0.j,\n", - " 4.17953273+0.j, 1.8473096 +0.j, 1.8473096 +0.j, 3.6946192 +0.j,\n", - " 2.29178823+0.j, 2.29178823+0.j, 4.58357645+0.j, 1.87517442+0.j,\n", - " 1.87517442+0.j, 3.75034885+0.j, 2.783546 +0.j, 2.783546 +0.j,\n", - " 5.567092 +0.j, 2.92732552+0.j, 2.92732552+0.j, 5.85465104+0.j,\n", - " 1.76688304+0.j, 1.76688304+0.j, 3.53376608+0.j, 2.58345008+0.j,\n", - " 2.58345008+0.j, 5.16690015+0.j, 2.05778984+0.j, 2.05778984+0.j,\n", - " 4.11557968+0.j, 2.13608912+0.j, 2.13608912+0.j, 4.27217824+0.j,\n", - " 2.85119328+0.j, 2.85119328+0.j, 5.70238655+0.j, 1.14207212+0.j,\n", - " 1.14207212+0.j, 2.28414423+0.j, 1.1742586 +0.j, 1.1742586 +0.j,\n", - " 2.3485172 +0.j, 1.04043679+0.j, 1.04043679+0.j, 2.08087359+0.j,\n", - " 2.66523969+0.j, 2.66523969+0.j, 5.33047938+0.j, 2.5563135 +0.j,\n", - " 2.5563135 +0.j, 5.112627 +0.j, 2.7400243 +0.j, 2.7400243 +0.j,\n", - " 5.48004859+0.j, 2.95723668+0.j, 2.95723668+0.j, 5.91447337+0.j,\n", - " 2.59831713+0.j, 2.59831713+0.j, 5.19663426+0.j, 1.92295872+0.j,\n", - " 1.92295872+0.j, 3.84591745+0.j, 2.56105835+0.j, 2.56105835+0.j,\n", - " 5.12211671+0.j, 1.23654885+0.j, 1.23654885+0.j, 2.4730977 +0.j,\n", - " 2.27984204+0.j, 2.27984204+0.j, 4.55968409+0.j, 1.28670657+0.j,\n", - " 1.28670657+0.j, 2.57341315+0.j, 2.88933783+0.j, 2.88933783+0.j,\n", - " 5.77867567+0.j, 2.04369664+0.j, 2.04369664+0.j, 4.08739329+0.j,\n", - " 1.82932388+0.j, 1.82932388+0.j, 3.65864776+0.j, 1.52911122+0.j,\n", - " 1.52911122+0.j, 3.05822245+0.j, 2.54846738+0.j, 2.54846738+0.j,\n", - " 5.09693476+0.j, 1.91230066+0.j, 1.91230066+0.j, 3.82460133+0.j,\n", - " 2.1368679 +0.j, 2.1368679 +0.j, 4.2737358 +0.j, 1.0375796 +0.j,\n", - " 1.0375796 +0.j, 2.0751592 +0.j, 2.23527099+0.j, 2.23527099+0.j,\n", - " 4.47054199+0.j, 2.22419145+0.j, 2.22419145+0.j, 4.44838289+0.j,\n", - " 2.23386799+0.j, 2.23386799+0.j, 4.46773599+0.j, 2.88749616+0.j,\n", - " 2.88749616+0.j, 5.77499231+0.j, 2.3636406 +0.j, 2.3636406 +0.j,\n", - " 4.7272812 +0.j, 1.7190158 +0.j, 1.7190158 +0.j, 3.4380316 +0.j,\n", - " 1.87406391+0.j, 1.87406391+0.j, 3.74812782+0.j, 2.39526239+0.j,\n", - " 2.39526239+0.j, 4.79052478+0.j, 1.12045094+0.j, 1.12045094+0.j,\n", - " 2.24090189+0.j, 2.33353343+0.j, 2.33353343+0.j, 4.66706686+0.j,\n", - " 2.34127574+0.j, 2.34127574+0.j, 4.68255148+0.j, 1.42076512+0.j,\n", - " 1.42076512+0.j, 2.84153024+0.j, 1.2578526 +0.j, 1.2578526 +0.j,\n", - " 2.51570519+0.j, 1.6308567 +0.j, 1.6308567 +0.j, 3.2617134 +0.j])\n" + "Static coefficients: [ 0.26666667 -2.31428571 3.04761905]\n", + "L1 norm: 5.628571428571431\n", + "Approximate coefficients: [-0.24255546 -0.25744454 1.5 ]\n", + "L1 norm (approx): 2.0\n", + "Computing dynamic coefficients for time=3\n" ] } ], "source": [ - "import numpy as np\n", - "from qiskit.quantum_info import SparsePauliOp, Pauli\n", - "\n", + "# -------------------------Step 1: Map classical inputs-------------------------\n", + "L = 50\n", + "coupling_map = CouplingMap.from_line(L, bidirectional=False)\n", "\n", - "# Generate random coefficients for our XXZ Hamiltonian\n", + "# XXZ Hamiltonian with random couplings (Ref. [3])\n", "np.random.seed(0)\n", "even_edges = list(coupling_map.get_edges())[::2]\n", "odd_edges = list(coupling_map.get_edges())[1::2]\n", @@ -1880,174 +1120,49 @@ " num_qubits=L,\n", " )\n", "\n", - "print(hamiltonian)" - ] - }, - { - "cell_type": "markdown", - "id": "5f87dc5d", - "metadata": {}, - "source": [ - "For an observable we choose $Z_{24}Z_{25}$, as seen in the lower panel Fig. 5 of Ref. [\\[3\\]](#references)." - ] - }, - { - "cell_type": "code", - "execution_count": 35, - "id": "3bb3c684", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIZZIIIIIIIIIIIIIIIIIIIIIIII'],\n", - " coeffs=[1.+0.j])\n" - ] - } - ], - "source": [ "observable = SparsePauliOp.from_sparse_list(\n", " [(\"ZZ\", (L // 2 - 1, L // 2), 1.0)], num_qubits=L\n", ")\n", - "print(observable)" - ] - }, - { - "cell_type": "markdown", - "id": "8a0fb376", - "metadata": {}, - "source": [ - "#### Choose Trotter steps\n", - "The experiment showcased in Fig. 4 of Ref. [\\[3\\]](#references) uses $k_j = [2, 3, 4]$ symmetric Trotter steps of order $2$. We focus on the results for time $t=3$, where the MPF and a PF with a higher number of Trotter steps (6 in this case) have the same Trotter error. However, the MPF expectation value is calculated from circuits corresponding to the lower number of Trotter steps and thus shallower. In practice, even if the MPF and the deeper Trotter steps circuit have the same Trotter error, we expect the experimental expectation value calculated from the MPF circuits to be closer to the theory one, as it entails running shallower circuits less exposed to hardware noise compared to the circuit corresponding to the higher Trotter step PF." - ] - }, - { - "cell_type": "code", - "execution_count": 36, - "id": "06e08219", - "metadata": {}, - "outputs": [], - "source": [ + "\n", "total_time = 3\n", "mpf_trotter_steps = [2, 3, 4]\n", "order = 2\n", - "symmetric = True" - ] - }, - { - "cell_type": "markdown", - "id": "85ea15da", - "metadata": {}, - "source": [ - "#### Set up the LSE\n", - "Here we look at the static MPF coefficients for this problem." - ] - }, - { - "cell_type": "code", - "execution_count": 37, - "id": "7fee3d13", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "The static coefficients associated with the ansatze are: [ 0.26666667 -2.31428571 3.04761905]\n", - "L1 norm: 5.628571428571431\n" - ] - } - ], - "source": [ + "symmetric = True\n", + "\n", + "# Static coefficients\n", "lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)\n", "mpf_coeffs = lse.solve()\n", - "print(\n", - " f\"The static coefficients associated with the ansatze are: {mpf_coeffs}\"\n", - ")\n", - "print(\"L1 norm:\", np.linalg.norm(mpf_coeffs, ord=1))" - ] - }, - { - "cell_type": "code", - "execution_count": 38, - "id": "7ee6bdf3", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[-0.24255546 -0.25744454 1.5 ]\n", - "L1 norm of the approximate coefficients: 2.0\n" - ] - } - ], - "source": [ + "print(f\"Static coefficients: {mpf_coeffs}\")\n", + "print(f\"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}\")\n", + "\n", "model_approx, coeffs_approx = setup_sum_of_squares_problem(\n", " lse, max_l1_norm=2.0\n", ")\n", "model_approx.solve()\n", - "print(coeffs_approx.value)\n", - "print(\n", - " \"L1 norm of the approximate coefficients:\",\n", - " np.linalg.norm(coeffs_approx.value, ord=1),\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "87262bc4", - "metadata": {}, - "source": [ - "#### Dynamic coefficients" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b2fec4a5", - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - } - ], - "source": [ - "# Create approximate time-evolution circuits\n", + "print(f\"Approximate coefficients: {coeffs_approx.value}\")\n", + "print(f\"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}\")\n", + "\n", + "# -------------------------Dynamic coefficients-------------------------\n", "single_2nd_order_circ = generate_time_evolution_circuit(\n", " hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)\n", ")\n", - "single_2nd_order_circ = pm.run(single_2nd_order_circ) # collect XX and YY\n", + "single_2nd_order_circ = pm.run(single_2nd_order_circ)\n", "\n", - "# Find layers in the circuit\n", "layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)\n", - "\n", - "# Create tensor network models\n", "models = [\n", " LayerModel.from_quantum_circuit(layer, conserve=\"Sz\") for layer in layers\n", "]\n", "\n", - "# Create the time-evolution object\n", "approx_factory = partial(\n", " LayerwiseEvolver,\n", " layers=models,\n", " options={\n", " \"preserve_norm\": False,\n", - " \"trunc_params\": {\n", - " \"chi_max\": 64,\n", - " \"svd_min\": 1e-8,\n", - " \"trunc_cut\": None,\n", - " },\n", + " \"trunc_params\": {\"chi_max\": 64, \"svd_min\": 1e-8, \"trunc_cut\": None},\n", " \"max_delta_t\": 4,\n", " },\n", ")\n", "\n", - "# Create exact time-evolution circuits\n", "single_4th_order_circ = generate_time_evolution_circuit(\n", " hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)\n", ")\n", @@ -2057,18 +1172,13 @@ " for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)\n", "]\n", "\n", - "# Create the time-evolution object\n", "exact_factory = partial(\n", " LayerwiseEvolver,\n", " layers=exact_model_layers,\n", " dt=0.1,\n", " options={\n", " \"preserve_norm\": False,\n", - " \"trunc_params\": {\n", - " \"chi_max\": 64,\n", - " \"svd_min\": 1e-8,\n", - " \"trunc_cut\": None,\n", - " },\n", + " \"trunc_params\": {\"chi_max\": 64, \"svd_min\": 1e-8, \"trunc_cut\": None},\n", " \"max_delta_t\": 3,\n", " },\n", ")\n", @@ -2080,8 +1190,8 @@ "\n", "mps_initial_state = MPS_neel_state(models[0].lat)\n", "\n", - "\n", - "lse = setup_dynamic_lse(\n", + "print(f\"Computing dynamic coefficients for time={total_time}\")\n", + "lse_dyn = setup_dynamic_lse(\n", " mpf_trotter_steps,\n", " total_time,\n", " identity_factory,\n", @@ -2089,113 +1199,75 @@ " approx_factory,\n", " mps_initial_state,\n", ")\n", - "problem, coeffs = setup_frobenius_problem(lse)\n", + "problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)\n", "try:\n", " problem.solve()\n", - " mpf_dynamic_coeffs = coeffs.value\n", + " mpf_dynamic_coeffs = coeffs_dyn.value\n", "except Exception as error:\n", - " print(error, \"Calculation Failed for time\", total_time)\n", - "print(\"\")" - ] - }, - { - "cell_type": "markdown", - "id": "43829572", - "metadata": {}, - "source": [ - "#### Construct each of the Trotter circuits in our MPF decomposition" - ] - }, - { - "cell_type": "code", - "execution_count": 41, - "id": "b79c5547", - "metadata": {}, - "outputs": [], - "source": [ - "from qiskit.synthesis import SuzukiTrotter\n", - "from qiskit_addon_utils.problem_generators import (\n", - " generate_time_evolution_circuit,\n", - ")\n", - "from qiskit import QuantumCircuit\n", - "\n", + " mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))\n", + " print(error, \"Calculation Failed\")\n", "\n", + "# -------------------------Step 1 (cont): Build circuits-------------------------\n", "mpf_circuits = []\n", "for k in mpf_trotter_steps:\n", - " # Initial state preparation |1010..>\n", " circuit = QuantumCircuit(L)\n", " circuit.x([i for i in range(L) if i % 2])\n", - "\n", " trotter_circ = generate_time_evolution_circuit(\n", " hamiltonian,\n", " synthesis=SuzukiTrotter(reps=k, order=order),\n", " time=total_time,\n", " )\n", - "\n", " circuit.compose(trotter_circ, qubits=range(L), inplace=True)\n", + " mpf_circuits.append(circuit)\n", "\n", - " mpf_circuits.append(circuit)" - ] - }, - { - "cell_type": "markdown", - "id": "ffd915d7", - "metadata": {}, - "source": [ - "#### Construct Trotter circuit with comparable Trotter error to MPF" - ] - }, - { - "cell_type": "code", - "execution_count": 42, - "id": "b22143ff", - "metadata": {}, - "outputs": [], - "source": [ - "k = 6\n", - "\n", - "# Initial state preparation |1010..>\n", + "# Baseline \"single deep circuit\" comparison run with k=6 Trotter steps.\n", + "# (Its two-qubit depth is comparable to the deepest MPF constituent plus the\n", + "# overhead of running multiple circuits; it does NOT target the MPF's\n", + "# effective Trotter error.)\n", "comp_circuit = QuantumCircuit(L)\n", "comp_circuit.x([i for i in range(L) if i % 2])\n", - "\n", - "\n", "trotter_circ = generate_time_evolution_circuit(\n", " hamiltonian,\n", - " synthesis=SuzukiTrotter(reps=k, order=order),\n", + " synthesis=SuzukiTrotter(reps=6, order=order),\n", " time=total_time,\n", ")\n", - "\n", "comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)\n", - "\n", - "\n", "mpf_circuits.append(comp_circuit)" ] }, { "cell_type": "markdown", - "id": "25f03e39", + "id": "b5d388b1", "metadata": {}, "source": [ - "### Step 2: Optimize problem for quantum hardware execution" + "### Step 2: Optimize problem for quantum hardware execution\n", + "\n", + "At 50 qubits we need to be deliberate about which physical qubits we use. We filter the backend's coupling map to exclude qubits with high measurement error, low $T_2$ times, or high two-qubit gate error, then transpile onto the filtered topology. Because calibration data varies across backends and over time, we use an adaptive approach: start with strict thresholds and progressively relax them until we find a connected component large enough for our circuit." ] }, { "cell_type": "code", - "execution_count": null, - "id": "00b29fa6", + "execution_count": 23, + "id": "05bad997", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "size of largest component 73\n" + "\n", + "Thresholds (meas=0.055, t2=30, 2q=0.01): largest component = 15 qubits\n", + "Thresholds (meas=0.08, t2=20, 2q=0.015): largest component = 89 qubits\n", + "\n", + "Using connected component with 89 qubits\n" ] } ], "source": [ - "import copy\n", - "from qiskit.transpiler import Target, CouplingMap\n", + "service = QiskitRuntimeService()\n", + "# backend = service.least_busy(min_num_qubits=127)\n", + "backend = service.backend(\"ibm_fez\")\n", + "print(backend)\n", "\n", "target = backend.target\n", "instruction_2q = \"cz\"\n", @@ -2203,71 +1275,96 @@ "cmap = target.build_coupling_map(filter_idle_qubits=True)\n", "cmap_list = list(cmap.get_edges())\n", "\n", - "max_meas_err = 0.055\n", - "min_t2 = 30\n", - "max_twoq_err = 0.01\n", - "\n", - "# Remove qubits with bad measurement or t2\n", - "cust_cmap_list = copy.deepcopy(cmap_list)\n", - "for q in range(target.num_qubits):\n", - " meas_err = target[\"measure\"][(q,)].error\n", - " if target.qubit_properties[q].t2 is not None:\n", - " t2 = target.qubit_properties[q].t2 * 1e6\n", - " else:\n", - " t2 = 0\n", - " if meas_err > max_meas_err or t2 < min_t2:\n", - " # print(q)\n", - " for q_pair in cmap_list:\n", - " if q in q_pair:\n", - " try:\n", - " cust_cmap_list.remove(q_pair)\n", - " except ValueError:\n", - " continue\n", - "\n", - "# Remove qubits with bad 2q gate or t2\n", - "for q in cmap_list:\n", - " twoq_gate_err = target[instruction_2q][q].error\n", - " if twoq_gate_err > max_twoq_err:\n", - " # print(q)\n", - " for q_pair in cmap_list:\n", - " if q == q_pair:\n", - " try:\n", - " cust_cmap_list.remove(q_pair)\n", - " except ValueError:\n", - " continue\n", - "\n", - "\n", - "cust_cmap = CouplingMap(cust_cmap_list)\n", + "min_qubits_needed = (\n", + " L # We need at least this many qubits in a connected component\n", + ")\n", + "\n", + "\n", + "def filter_coupling_map(\n", + " cmap_list, target, max_meas_err, min_t2, max_twoq_err\n", + "):\n", + " \"\"\"Filter coupling map edges based on qubit quality thresholds.\"\"\"\n", + " filtered = deepcopy(cmap_list)\n", + " for q in range(target.num_qubits):\n", + " meas_err = target[\"measure\"][(q,)].error\n", + " t2 = target.qubit_properties[q].t2\n", + " t2 = t2 * 1e6 if t2 is not None else 0\n", + " if meas_err > max_meas_err or t2 < min_t2:\n", + " for q_pair in cmap_list:\n", + " if q in q_pair:\n", + " try:\n", + " filtered.remove(q_pair)\n", + " except ValueError:\n", + " continue\n", + " for q in cmap_list:\n", + " twoq_gate_err = target[instruction_2q][q].error\n", + " if twoq_gate_err > max_twoq_err:\n", + " for q_pair in cmap_list:\n", + " if q == q_pair:\n", + " try:\n", + " filtered.remove(q_pair)\n", + " except ValueError:\n", + " continue\n", + " return filtered\n", + "\n", + "\n", + "# Start with strict thresholds and relax until we find a large enough component\n", + "threshold_configs = [\n", + " (0.055, 30, 0.010), # strict\n", + " (0.080, 20, 0.015), # moderate\n", + " (0.120, 10, 0.020), # relaxed\n", + "]\n", + "\n", + "for max_meas_err, min_t2, max_twoq_err in threshold_configs:\n", + " cust_cmap_list = filter_coupling_map(\n", + " cmap_list, target, max_meas_err, min_t2, max_twoq_err\n", + " )\n", + " cust_cmap = CouplingMap(cust_cmap_list)\n", + " sorted_components = sorted(\n", + " [\n", + " list(comp.physical_qubits)\n", + " for comp in cust_cmap.connected_components()\n", + " ],\n", + " key=len,\n", + " reverse=True,\n", + " )\n", + " largest = len(sorted_components[0]) if sorted_components else 0\n", + " print(\n", + " f\"Thresholds (meas={max_meas_err}, t2={min_t2}, 2q={max_twoq_err}): \"\n", + " f\"largest component = {largest} qubits\"\n", + " )\n", + " if largest >= min_qubits_needed:\n", + " break\n", + "\n", + "if largest < min_qubits_needed:\n", + " # Fall back to the full coupling map if no threshold works\n", + " print(\n", + " f\"Warning: filtering leaves no component >= {min_qubits_needed} qubits. Using full coupling map.\"\n", + " )\n", + " cust_cmap = cmap\n", "\n", "cust_target = Target.from_configuration(\n", - " basis_gates=backend.configuration().basis_gates\n", - " + [\"measure\"], # or whatever new set of gates\n", + " basis_gates=backend.configuration().basis_gates + [\"measure\"],\n", " coupling_map=cust_cmap,\n", ")\n", "\n", - "sorted_components = sorted(\n", - " [list(comp.physical_qubits) for comp in cust_cmap.connected_components()],\n", - " reverse=True,\n", - ")\n", - "print(\"size of largest component\", len(sorted_components[0]))" + "print(f\"\\nUsing connected component with {largest} qubits\")" ] }, { "cell_type": "code", - "execution_count": 104, - "id": "3caa3939", + "execution_count": 24, + "id": "0f4303ee", "metadata": {}, "outputs": [], "source": [ - "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", - "\n", "transpiler = generate_preset_pass_manager(\n", " optimization_level=3, target=cust_target\n", ")\n", "\n", "transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]\n", "\n", - "qubits_layouts = [\n", + "qubit_layouts = [\n", " [\n", " idx\n", " for idx, qb in circuit.layout.initial_layout.get_physical_bits().items()\n", @@ -2277,16 +1374,13 @@ "]\n", "\n", "transpiled_circuits = []\n", - "for circuit, layout in zip(mpf_circuits, qubits_layouts):\n", + "for circuit, layout in zip(mpf_circuits, qubit_layouts):\n", " transpiler = generate_preset_pass_manager(\n", " optimization_level=3, backend=backend, initial_layout=layout\n", " )\n", " transpiled_circuit = transpiler.run(circuit)\n", " transpiled_circuits.append(transpiled_circuit)\n", "\n", - "\n", - "# transform the observable defined on virtual qubits to\n", - "# an observable defined on all physical qubits\n", "isa_observables = [\n", " observable.apply_layout(circ.layout) for circ in transpiled_circuits\n", "]" @@ -2294,25 +1388,27 @@ }, { "cell_type": "markdown", - "id": "1a983973", + "id": "b578c285", "metadata": {}, "source": [ - "### Step 3: Execute using Qiskit primitives" + "### Step 3: Execute using Qiskit primitives\n", + "\n", + "Running deeper circuits on real hardware requires aggressive error mitigation. We enable dynamical decoupling, gate and measurement twirling, measurement error mitigation, and zero-noise extrapolation (ZNE). Note that the ZNE noise factors we use here (`1, 1.2, 1.4`) are smaller than in a shallow-circuit scenario, since the deeper MPF constituents are already close to the noise threshold and large noise amplifications would push them past the point where ZNE extrapolation is reliable.\n", + "\n", + "We submit all four circuits (three MPF constituents at $k_j = [2, 3, 4]$ plus the $k = 6$ baseline) in a single Estimator job." ] }, { "cell_type": "code", - "execution_count": null, - "id": "83883725", + "execution_count": 25, + "id": "e2722b61", "metadata": {}, "outputs": [], "source": [ - "from qiskit_ibm_runtime import EstimatorV2 as Estimator\n", - "\n", "estimator = Estimator(mode=backend)\n", "estimator.options.default_shots = 30000\n", "\n", - "# Set simple error suppression/mitigation options\n", + "# Error suppression/mitigation\n", "estimator.options.dynamical_decoupling.enable = True\n", "estimator.options.twirling.enable_gates = True\n", "estimator.options.twirling.enable_measure = True\n", @@ -2325,7 +1421,7 @@ "estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)\n", "estimator.options.resilience.zne.extrapolator = \"linear\"\n", "\n", - "estimator.options.environment.job_tags = [\"mpf large\"]\n", + "estimator.options.environment.job_tags = [\"TUT_MPF\"]\n", "\n", "job_50 = estimator.run(\n", " [\n", @@ -2337,24 +1433,26 @@ }, { "cell_type": "markdown", - "id": "0fee3da4", + "id": "afc0c029", "metadata": {}, "source": [ - "### Step 4: Post-process and return result in desired classical format" + "### Step 4: Post-process and return result in desired classical format\n", + "\n", + "We pull the per-circuit expectation values and standard deviations from the job result, then combine them with each set of MPF coefficients exactly as in the small-scale example: $\\langle A \\rangle_{\\text{MPF}} = \\sum_j x_j \\, \\langle A \\rangle_{k_j}$, with propagated variance $\\sigma^2 = \\sum_j x_j^2 \\sigma_{k_j}^2$." ] }, { "cell_type": "code", - "execution_count": 46, - "id": "5d2fc41c", + "execution_count": 26, + "id": "a924d79c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "[array(-0.08034071), array(-0.00605026), array(-0.15345759), array(-0.18127293)]\n", - "[array(0.04482517), array(0.03438413), array(0.21540776), array(0.21520829)]\n" + "[array(-0.17503803), array(-0.12294363), array(-0.33742479), array(-0.15413359)]\n", + "[array(0.01520694), array(0.05557148), array(0.01454295), array(0.10680865)]\n" ] } ], @@ -2369,17 +1467,17 @@ }, { "cell_type": "code", - "execution_count": 52, - "id": "bff48c28", + "execution_count": 27, + "id": "1071de0d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Exact static MPF expectation value: -0.47510243192011536 +- 0.6613940032465087\n", - "Approximate static MPF expectation value: -0.20914170384216998 +- 0.32341567460419135\n", - "Dynamic MPF expectation value: -0.07994951978722761 +- 0.07423091963310202\n" + "Exact static MPF expectation value: -0.7904923379680913 +- 0.13609157812004166\n", + "Approximate static MPF expectation value: -0.4320295863244896 +- 0.026346772829133787\n", + "Dynamic MPF expectation value: -0.2118021477340412 +- 0.01980723126846985\n" ] } ], @@ -2425,14 +1523,14 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "d751af7c", + "execution_count": 28, + "id": "64360d85", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "\"Output" + "\"Output" ] }, "metadata": {}, @@ -2441,7 +1539,6 @@ ], "source": [ "sym = {2: \"^\", 3: \"s\", 4: \"p\"}\n", - "# Get expectation values at all times for each Trotter step\n", "for k, step in enumerate(mpf_trotter_steps):\n", " plt.errorbar(\n", " k,\n", @@ -2454,7 +1551,6 @@ " label=f\"{mpf_trotter_steps[k]} Trotter steps\",\n", " )\n", "\n", - "\n", "plt.errorbar(\n", " 3,\n", " evs[-1],\n", @@ -2466,7 +1562,6 @@ " label=\"6 Trotter steps\",\n", ")\n", "\n", - "\n", "plt.errorbar(\n", " 4,\n", " evs[:3] @ mpf_coeffs,\n", @@ -2502,9 +1597,8 @@ " y=exact_obs, linestyle=\"--\", color=\"red\", label=\"Exact time-evolution\"\n", ")\n", "\n", - "\n", "plt.title(\n", - " f\"Expectation values for (ZZ,{(L//2-1, L//2)}) at time {total_time} for the different methods \"\n", + " f\"$\\\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\\\rangle$ at time {total_time} for the different methods\"\n", ")\n", "plt.xlabel(\"Method\")\n", "plt.ylabel(\"Expectation Value\")\n", @@ -2516,28 +1610,49 @@ }, { "cell_type": "markdown", - "id": "04200168", + "id": "31effe5f", "metadata": {}, "source": [ - "When executing circuits on hardware, we might encounter additional challenges in obtaining accurate expectation values due to the presence of hardware noise. This is not accounted for in the MPF formalism and could play against the MPF solution. For example, this could be the reason for the failure of the dynamic coefficients to provide a better estimate of the expectation value compared to the approximate static coefficient in the plot. That is, the approximate evolver, which simulates the approximate circuit, does not accurately reflect the results obtained by executing the approximate circuits in the presence of hardware noise. For these reasons, it is recommended to combine different error mitigation techniques in order to get results as close as possible to the ideal values for each of the product formulas. This will show consistent benefits from the MPF approach.\n", + "A few observations about the hardware results above:\n", + "\n", + "- **Coefficient norm matters more than mathematical optimality.** The exact-static MPF has $\\|x\\|_1 = 5.63$, and its final estimate ends up far from the true value — the large coefficient norm amplifies the residual gate noise, decoherence, and ZNE error on each $\\langle A \\rangle_{k_j}$ by roughly the same factor, overwhelming the Trotter-error cancellation we gained. The approximate-static MPF (capped at $\\|x\\|_1 = 2$) and the dynamic MPF (whose coefficients here are similarly modest in magnitude) avoid this overwhelm and land much closer to the reference.\n", + "\n", + "- **The dynamic MPF is the best estimator in this run.** Tailoring the coefficients to this specific Hamiltonian, initial state, and time gives an answer closer to the exact reference than any single Trotter circuit — including the deeper $k = 6$ baseline — and is substantially closer than either static MPF. This is the regime MPFs are designed for: shallow constituent circuits combined with well-chosen, small-norm coefficients.\n", + "\n", + "- **Individual product formulas are still competitive on hardware.** Note that the single-circuit baselines ($k = 2, 3, 6$) all sit within roughly $0.1$ of the exact reference, while the static-exact MPF is off by more than $0.5$. This is the failure mode the small-scale Background discussion warned about: an MPF that mathematically cancels Trotter error can be *worse* than a single product formula once hardware noise is amplified by $\\|x\\|_1$.\n", + "\n", + "- **The chosen Trotter steps are partly outside the safe regime.** At $t = 3$ with $k_{\\min} = 2$, the ratio $t/k_{\\min} = 1.5 > 1$, which is past the convergence threshold discussed in the Background. The leading-error cancellation the static MPF relies on is therefore less effective, and the dynamic MPF's tensor-network reference for the \"exact\" state is also strained — even so, the dynamic MPF performs well here because its coefficients adapt to whatever errors the approximate evolver does see.\n", "\n", - "Overall, the approximate static coefficients still give a more accurate solution than the product formula with higher number of Trotter steps with the same amount of Trotter error in the noiseless setting.\n", + "The practical takeaway is that on hardware, MPFs should be paired with strong error mitigation on each individual $\\langle A \\rangle_{k_j}$, the coefficient $L_1$-norm should be kept modest (use the approximate solver, or prefer the dynamic MPF when its tensor-network setup is tractable), and the Trotter steps $k_j$ should be chosen so that $t/k_{\\min} \\lesssim 1$. With those choices, as we see here for the dynamic MPF, the MPF can recover the depth-versus-accuracy advantage shown in Ref. [\\[3\\]](#references). Note also that individual runs are noisy — on a different submission of the same job (or a different backend), the relative ordering of these estimators can shift; the qualitative trend that small-$\\|x\\|_1$ MPFs do well and large-$\\|x\\|_1$ MPFs are amplified by hardware noise is the robust one." + ] + }, + { + "cell_type": "markdown", + "id": "5ac2f8a8", + "metadata": {}, + "source": [ + "## Next steps\n", + "\n", "\n", - "It is also important to note that in the example that reproduces the experiment in Ref. [\\[3\\]](#references), the time point $t=3$ is beyond the limit at which it is expected that the PF with $k=2$ behaves well, which is $t/k>1$ as discussed in this [guide](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html)." + "If you found this work interesting, you might be interested in the following material:\n", + "- [How to choose the Trotter steps for an MPF](https://qiskit.github.io/qiskit-addon-mpf/how_tos/choose_trotter_steps.html) — practical guidance on selecting $k_j$ values to avoid instabilities\n", + "- [How to use the approximate model](https://qiskit.github.io/qiskit-addon-mpf/how_tos/using_approximate_model.html) — tuning the $L_1$-norm constraint and solver options for the approximate static MPF\n", + "- [`qiskit-addon-mpf` API reference](https://qiskit.github.io/qiskit-addon-mpf/) — full documentation for static, dynamic, and backend modules\n", + "" ] }, { "cell_type": "markdown", - "id": "78e91ddc", + "id": "70be41e1", "metadata": {}, "source": [ "## References\n", "\n", - "[1] [Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. (2023). Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 7, 1067](https://quantum-journal.org/papers/q-2023-07-25-1067/).\n", + "\\[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. [Quantum, 7, 1067 (2023)](https://quantum-journal.org/papers/q-2023-07-25-1067/)\n", "\n", - "[2] [Zhuk, S., Robertson, N. F., & Bravyi, S. (2024). Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 6(3), 033309](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309).\n", + "\\[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. [Physical Review Research, 6(3), 033309 (2024)](https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.6.033309)\n", "\n", - "[3] [Robertson, N. F., et al. (2024). Tensor network enhanced dynamic multiproduct formulas. arXiv preprint arXiv:2407.17405](https://arxiv.org/abs/2407.17405)." + "\\[3] Robertson, N. F., et al. Tensor network enhanced dynamic multiproduct formulas. [arXiv:2407.17405 (2024)](https://arxiv.org/abs/2407.17405)" ] } ], @@ -2558,11 +1673,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" - }, - "vscode": { - "interpreter": { - "hash": "72460b19b35c51a2148375679c8e282bc3cd3c5550c209f0a38b7d09cded82d9" - } } }, "nbformat": 4, diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/25ce07a6-1.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/25ce07a6-1.avif deleted file mode 100644 index 692112f6ed2..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/25ce07a6-1.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/2da9c948-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/2da9c948-0.avif deleted file mode 100644 index 68d7536c73e..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/2da9c948-0.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/34bf68ac-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/34bf68ac-0.avif deleted file mode 100644 index 213fcfe1a31..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/34bf68ac-0.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/35042576-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/35042576-0.avif new file mode 100644 index 00000000000..c26b6ed1c13 Binary files /dev/null and b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/35042576-0.avif differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/64360d85-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/64360d85-0.avif new file mode 100644 index 00000000000..80440a4c320 Binary files /dev/null and b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/64360d85-0.avif differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/87d2ac0c-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/87d2ac0c-0.avif deleted file mode 100644 index 7a180dca513..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/87d2ac0c-0.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/92dc20a7-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/92dc20a7-0.avif deleted file mode 100644 index 3331e1f5616..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/92dc20a7-0.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/a3eefe73-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/a3eefe73-0.avif deleted file mode 100644 index 3c6fc0e54a6..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/a3eefe73-0.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/c5d8e90b-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/c5d8e90b-0.avif deleted file mode 100644 index 7ee9f11f065..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/c5d8e90b-0.avif and /dev/null differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/c7ee61e7-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/c7ee61e7-0.avif new file mode 100644 index 00000000000..cc832a46d24 Binary files /dev/null and b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/c7ee61e7-0.avif differ diff --git a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/d751af7c-0.avif b/public/docs/images/tutorials/multi-product-formula/extracted-outputs/d751af7c-0.avif deleted file mode 100644 index 8364ca7036a..00000000000 Binary files a/public/docs/images/tutorials/multi-product-formula/extracted-outputs/d751af7c-0.avif and /dev/null differ