{ "cells": [ { "cell_type": "markdown", "id": "d2c31ae8", "metadata": {}, "source": [ "---\n", "title: Shadow tomography of the Greenberger-Horne-Zeilinger state\n", "description: This tutorial demonstrates the shadow tomography protocol by estimating observables of the Greenberger-Horne-Zeilinger state with up to 100 qubits\n", "---\n", "\n", "{/* cspell:ignore matchgate binom bitstr elementwise correlator vmin vmax ests postselects Kueng Zhao Miyake Bertoni Haferkamp Hinsche */}\n", "\n", "# Shadow tomography of the Greenberger-Horne-Zeilinger state\n", "\n", "*Usage estimate: Under 10 seconds on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "8bf80006", "metadata": {}, "source": [ "## Learning outcomes\n", "\n", "* How shadow tomography constructs a classical shadow of a quantum state and why it requires exponentially fewer measurements than full state tomography\n", "* How to implement local Pauli shadow tomography using Qiskit to estimate single-qubit magnetizations, two-point correlators, and two-qubit reduced density matrices\n", "* How to use [qiskit-paulice](/docs/addons/qiskit-addon-paulice) for low-overhead error detection using spacetime codes on a large-scale hardware experiment\n", "\n", "## Prerequisites\n", "\n", "* [Measurements](/learning/courses/general-formulation-of-quantum-information/general-measurements/introduction)\n", "* [Quantum state tomography](/learning/courses/general-formulation-of-quantum-information/general-measurements/discrimination-and-tomography)\n", "* [Error detection using spacetime codes](/docs/tutorials/spacetime-codes)\n", "\n", "## Background\n", "\n", "Full quantum state tomography of an $N$-qubit system requires exponentially many measurements to estimate all $4^N -1$ independent parameters of the density matrix. For a system of 50 qubits, for example, this process requires on the order of $10^{30}$ shots, which is more than a billion times the number of stars in the observable universe. Shadow tomography [\\[1\\]](#references) is a quantum information protocol that constructs a “classical shadow” of the quantum state, from which target observables can be estimated by using a number of shots that scales logarithmically in the number of observables and is independent of system size $N$, making it relevant across many areas in quantum computing, including energy estimation, entanglement entropy and correlation function calculations, kernel estimation in quantum machine learning, and others.\n", "\n", "This tutorial shows how to use shadow tomography to estimate local observables, such as the single-qubit magnetization $\\langle Z_i \\rangle$ for qubit $i$, the two-point correlation function $\\langle Z_i Z_j\\rangle$, and the two-qubit reduced density matrix (2-RDM), of the highly entangled Greenberger-Horne-Zeilinger (GHZ) state. The tutorial first illustrates a small-scale shadow tomography workflow on a 10-qubit GHZ state by using statevector simulation. Then it shows a 100-qubit demonstration on real quantum hardware.\n", "\n", "### Shadow tomography protocol\n", "\n", "The shadow tomography protocol has the following steps:\n", "\n", "1. Apply a random unitary $U$, drawn from a unitary ensemble, to the initial quantum state.\n", "2. Measure the outcome in the computational basis, obtaining a bitstring $\\ket{b}$.\n", "3. Store the classical snapshot, $\\hat{\\rho} = \\mathcal{M}^{-1}(U^{\\dagger}\\ket{b}\\bra{b}U)$, where $\\mathcal{M}^{-1}$ is the inverse of the measurement channel $\\mathcal{M}(\\rho)=\\mathbb{E}[U^{\\dagger}\\ket{b}\\bra{b}U]$.\n", "4. Repeat the process $T$ times, where the collection of $T$ snapshots $\\{\\hat{\\rho}_t\\}$ forms a classical shadow $S(\\rho, T)$.\n", "5. Estimate observable $o$ of operator $O$ as $\\frac{1}{T}\\sum_{t=1}^{T}\\text{Tr}(O\\hat{\\rho}_t)$.\n", "\n", "### Shadow tomography variants\n", "\n", "The choice of the unitary ensemble, that is, basis, in step 1 of the protocol, affects the sample complexity, variance of estimates, and information completeness of a given shadow tomography. Many ensemble groups have been considered in the literature, each with trade-offs between complexity and observable estimation efficiency. The properties of the most commonly used variants are summarized in the following table, although this list is not exhaustive.\n", "\n", "| Variant | Ensemble | Sample complexity | Circuit depth | Suitable observables |\n", "| ------------------------------------------------------------- | --------------------------------- | ------------------------------------------------------ | --------------------------- | ----------------------------- |\n", "| local Pauli [\\[1\\]](#references) | single-qubit Paulis $\\{X, Y, Z\\}$ | $\\mathcal{O}(3^k\\text{log}M/\\varepsilon^2)$ | 1 | local/few-body |\n", "| global Clifford [\\[1\\]](#references) | $N$-qubit Clifford | $\\mathcal{O}(\\text{Tr}(O^2)\\text{log}M/\\varepsilon^2)$ | $\\mathcal{O}(N\\text{log}N)$ | low-rank, fidelity |\n", "| fermionic Gaussian [\\[2\\]](#references), [\\[3\\]](#references) | matchgate circuits | $\\mathcal{O}(\\binom{\\mu}{k}\\frac{1}{\\varepsilon^2})$ | $\\mathcal{O}(N^2)$ | fermionic |\n", "| shallow shadows [\\[4\\]](#references), [\\[5\\]](#references) | low-depth circuits | $\\mathcal{O}(k2^k\\text{log}M/\\varepsilon^2)$ | $\\mathcal{O}(\\text{log}N)$ | non-local (hardware-friendly) |\n", "\n", "The terms used in the table are defined as follows:\n", "\n", "* \"sample complexity\" refers to the number of measurements needed for a given shadow tomography variant.\n", "* \"circuit depth\" describes how expensive the measurements are to implement on quantum hardware.\n", "* $N$ is the number of qubits.\n", "* $k$ is the locality of the observable (for example, $\\langle Z_i \\rangle$ has $k=1$ and $\\langle Z_i Z_j\\rangle$ has $k=2$).\n", "* $M$ is the number of observables to be estimated simultaneously.\n", "* $\\varepsilon$ is the additive error tolerance.\n", "* $\\mu$ is the number of fermions.\n", "\n", "Though the local Pauli ensemble is the most hardware-friendly, it pays a $3^k$ sample complexity penalty, which makes it ill-suited for non-local property estimation. More generally, the local-shadow sample-complexity bound scales as $4^k$, improving to $3^k$ in the special case of tensor-product observables.\n", "\n", "Conversely, while the global Clifford ensemble can be used to estimate non-local observables, its circuit depth is not favorable. Fermionic Gaussian shadows offer utility in simulations of many-body physics, chemistry, and materials, but require matchgate circuits, which might carry significant overhead on general purpose gate-based devices.\n", "\n", "Shallow shadows has a circuit depth between the local Pauli and global Clifford, making it a highly promising variant for practical implementation, though there are still some open questions regarding optimal circuit depth for a given observable.\n", "\n", "As the frontier of shadow tomography research continues expanding, some of the most relevant problems to address include optimally choosing the measurement ensemble for a target set of observables, improving sample complexity bounds, and extending the protocol to time-evolved observables. Nonetheless, shadow tomography still presents one of the leading workflows for extracting maximal classical information from a quantum state using a minimum number of measurements, which makes it ideal for characterizing quantum states at scale where both full tomography and classical simulation are impossible.\n", "\n", "### The local Pauli ensemble\n", "\n", "This tutorial uses the random Pauli basis, that is, single-qubit Clifford gates, to collect shadows and estimate properties of the GHZ state. Note that the Pauli measurements are a great choice for estimating local observables (such as spin density, local correlation functions, and position probability density), but not global ones (such as total energy, momentum, and charge).\n", "\n", "Here, a random Pauli basis $\\{X, Y, Z\\}$ per qubit is independently chosen and the corresponding single-qubit rotation is applied, followed by measurement:\n", "\n", "$$\n", "\n", "\\begin{equation}\n", "U = U_0 \\otimes U_1 \\otimes \\cdots \\otimes U_{n-1}, U_q \\in \\{I, H, S^{\\dagger}H\\}\n", "\\end{equation}\n", "\\tag{1}\n", "\n", "$$\n", "\n", "For this ensemble, the measurement channel's inverse is:\n", "\n", "$$\n", "\n", "\\begin{equation}\n", "\\hat{\\rho}_q=3U_q^{\\dagger}\\ket{b_q}\\bra{b_q}U_q - I\n", "\\end{equation}\n", "\\tag{2}\n", "$$\n", "\n", "The factor of $3$ comes from averaging over the three Pauli bases, which produces a depolarizing channel with parameter $\\frac{1}{3}$. The full estimator then factorizes as a product of single-qubit estimators. As shown in the table, for $M$ observables, each of locality $k$, the number of Pauli shadows needed to estimate all observables to additive error $\\varepsilon$ with high probability is $\\mathcal{O}(\\frac{3^k\\text{log}(M)}{\\varepsilon^2})$. The system size does not affect the measurement cost as long as the observables remain local, so exponentially many observables can be estimated simultaneously with only logarithmic cost.\n", "\n", "### Median of means\n", "\n", "A plain mean over $T$ snapshots reduces the variance of a single-shot estimate of an observable as $\\frac{1}{T}$, but is sensitive to outliers. Median of means can be used to remedy this. The process involves the following:\n", "\n", "* Splitting the $T$ snapshots into $K$ equally-sized groups\n", "* Computing the mean within each group\n", "* Returning the median across group means\n", "\n", "This tutorial performs local Pauli shadow tomography to estimate properties of the GHZ state, first using a 10-qubit simulated example and then demonstrating a 100-qubit hardware experiment.\n", "\n" ] }, { "cell_type": "markdown", "id": "55b94021", "metadata": {}, "source": [ "## Requirements\n", "\n", "Before starting this tutorial, ensure that you have the following installed:\n", "\n", "* Qiskit SDK v2.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n", "* IBM Quantum Compute v0.28 or later (`pip install qiskit-ibm-runtime`)\n", "* Qiskit Paulice (`pip install qiskit-paulice`)\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "7db2e559", "metadata": {}, "source": [ "## Setup\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "bc380c46", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from collections import deque, defaultdict\n", "from qiskit import QuantumCircuit, ClassicalRegister\n", "from qiskit.quantum_info import (\n", " Statevector,\n", " partial_trace,\n", " Operator,\n", " Pauli,\n", " SparsePauliOp,\n", ")\n", "from qiskit.transpiler import generate_preset_pass_manager\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "from qiskit_ibm_runtime import SamplerV2 as Sampler\n", "from qiskit_paulice.layout import get_low_overhead_ancillas" ] }, { "cell_type": "code", "execution_count": 2, "id": "20f94695", "metadata": {}, "outputs": [], "source": [ "seed = 42 # for reproducibility\n", "rng = np.random.default_rng(seed)" ] }, { "cell_type": "markdown", "id": "f1926db9", "metadata": {}, "source": [ "The following helper functions are used by the large-scale hardware example later in the tutorial:\n", "\n", "* `bad_cz`, `bad_readout`, `bad_coherence` — Query the backend calibration data to identify gates, qubits, and edges that exceed error or coherence thresholds, so they can be excluded from the circuit layout.\n", "* `parallel_ghz` — Builds a GHZ circuit of a given size using a breadth-first search over the coupling map, spreading CNOT layers in parallel to minimize circuit depth.\n", "* `weighted_coverage`, `active_wires`, `z_trace_backward` — Compute how much of the GHZ circuit's error surface is protected by a given set of spacetime parity checks, using the backward-propagation technique from [\\[6\\]](#references).\n", "* `add_measurements` — Appends GHZ-qubit and check-qubit measurements into separate classical registers so they can be postselected independently.\n", "* `build_shadow_circuits` — Produces one circuit per random Pauli basis by appending the appropriate single-qubit rotations before measurement.\n", "* `parse_hw_results_w_checks` — Unpacks the Sampler output into `(ghz_bits, check_bits, bases)` triples, one per shot.\n", "* `adapt_to_flat_bits`, `estimate_zz_nn_hw` — Map the GHZ-register bits back into the full-device index space and compute the median-of-means $\\langle Z_i Z_{i+1}\\rangle$ estimate.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "4ed737e3", "metadata": {}, "outputs": [], "source": [ "def bad_cz(target, threshold: float = 0.01) -> list:\n", " \"Return edges whose CZ error exceeds threshold.\"\n", " edges = []\n", " for edge in backend.target.build_coupling_map().get_edges():\n", " if (edge[1], edge[0]) not in edges:\n", " edges.append(edge)\n", " cz_errors = {edge: target[\"cz\"][edge].error for edge in edges}\n", " worst = sorted(cz_errors.items(), key=lambda x: x[1], reverse=True)\n", " return [list(edge) for edge, error in worst if error > threshold]\n", "\n", "\n", "def bad_readout(target, threshold: float = 0.01) -> list:\n", " \"Return nodes whose measurement error exceeds threshold.\"\n", " meas_errors = {\n", " node: target[\"measure\"][(node,)].error\n", " for node in range(backend.num_qubits)\n", " }\n", " worst = sorted(meas_errors.items(), key=lambda x: x[1], reverse=True)\n", " return [node for node, error in worst if error > threshold]\n", "\n", "\n", "def bad_coherence(target, threshold: float = 60) -> list:\n", " \"Return nodes whose T2 value is below threshold.\"\n", " t2s = {\n", " node: (target.qubit_properties[node].t2 or 0) * 1e6\n", " for node in range(backend.num_qubits)\n", " }\n", " return [\n", " node\n", " for node, val in sorted(t2s.items(), key=lambda x: x[1])\n", " if val < threshold\n", " ]\n", "\n", "\n", "def parallel_ghz(\n", " root: int, num_qubits: int, backend, bad_edges: list, skip: list\n", ") -> tuple:\n", " \"Build a GHZ circuit of size num_qubits using BFS from root, one neighbor per layer.\"\n", " edges = [\n", " e\n", " for e in [list(e) for e in backend.coupling_map.get_edges()]\n", " if e not in bad_edges\n", " and [e[1], e[0]] not in bad_edges\n", " and e[0] not in skip\n", " and e[1] not in skip\n", " ]\n", " adj = defaultdict(list)\n", " for u, v in edges:\n", " adj[u].append(v)\n", " adj[v].append(u)\n", "\n", " qc = QuantumCircuit(backend.num_qubits)\n", " visited = [root]\n", " queue = deque([root])\n", " explored = defaultdict(set)\n", " layers = []\n", " qc.h(root)\n", "\n", " while queue and len(visited) < num_qubits:\n", " layer = []\n", " for node in list(queue):\n", " queue.popleft()\n", " unvisited = [\n", " nb\n", " for nb in adj[node]\n", " if nb not in visited and nb not in explored[node]\n", " ]\n", " if unvisited:\n", " nb = unvisited[0]\n", " visited.append(nb)\n", " queue.append(node)\n", " queue.append(nb)\n", " explored[node].add(nb)\n", " layer.append((node, nb))\n", " if len(visited) == num_qubits:\n", " break\n", " if layer:\n", " layers.append(layer)\n", " else:\n", " break\n", "\n", " for layer in layers:\n", " for q1, q2 in layer:\n", " qc.cx(q1, q2)\n", " qc.barrier()\n", "\n", " return qc, visited, layers\n", "\n", "\n", "def weighted_coverage(\n", " layers: list, parities: list, w_idle: float = 0.2, w_gate: float = 0.8\n", ") -> float:\n", " \"Compute weighted fraction of wires covered by at least one parity check.\"\n", " wires = active_wires(layers)\n", " covered_by_any = {n_layer: set() for n_layer in range(len(layers))}\n", " for parity in parities:\n", " for n_layer, qs in z_trace_backward(layers, parity).items():\n", " covered_by_any[n_layer] |= qs\n", " covered_weight = total_weight = 0\n", " for n_layer in range(len(layers)):\n", " idle = wires[n_layer][\"idle\"]\n", " gate = wires[n_layer][\"gate\"]\n", " total_weight += w_idle * len(idle) + w_gate * len(gate)\n", " covered_weight += w_idle * len(\n", " covered_by_any[n_layer] & idle\n", " ) + w_gate * len(covered_by_any[n_layer] & gate)\n", " return covered_weight / total_weight if total_weight > 0 else 0\n", "\n", "\n", "def active_wires(layers: list) -> dict:\n", " \"\"\"\n", " Returns per-layer dict with two sets:\n", " - 'idle': activated wires that are idle in this layer\n", " - 'gate': activated wires that are control/target of a CNOT at this layer\n", " \"\"\"\n", " first_activation = {}\n", " for n_layer, layer in enumerate(layers):\n", " for c, t in layer:\n", " first_activation.setdefault(c, n_layer)\n", " first_activation.setdefault(t, n_layer)\n", " result = {}\n", " for n_layer in range(len(layers)):\n", " active = {\n", " q\n", " for q, n_layer0 in first_activation.items()\n", " if n_layer >= n_layer0\n", " }\n", " gate = {q for c, t in layers[n_layer] for q in (c, t)}\n", " result[n_layer] = {\"idle\": active - gate, \"gate\": gate}\n", " return result\n", "\n", "\n", "def z_trace_backward(layers: list, initial_Zs: list) -> dict:\n", " \"\"\"\n", " Backward propagate Zs with parity cancellation.\n", " Returns {layer: set of qubits with odd parity Z at that layer}.\n", " \"\"\"\n", " wires = active_wires(layers)\n", " support = set(initial_Zs)\n", " trace = {}\n", " for n_layer in range(len(layers) - 1, -1, -1):\n", " active = wires[n_layer][\"idle\"] | wires[n_layer][\"gate\"]\n", " trace[n_layer] = support & active\n", " new_support = set()\n", " for q in support:\n", " hit = False\n", " for c, t in layers[n_layer]:\n", " if q == t:\n", " new_support ^= {t, c}\n", " hit = True\n", " break\n", " elif q == c:\n", " new_support ^= {c}\n", " hit = True\n", " break\n", " if not hit:\n", " new_support ^= {q}\n", " support = new_support\n", " return trace\n", "\n", "\n", "def add_measurements(\n", " qc: QuantumCircuit, ghz_qubits: list, checks: list\n", ") -> QuantumCircuit:\n", " \"Add GHZ and check qubit measurements into separate classical registers.\"\n", " c1 = ClassicalRegister(len(ghz_qubits), \"c1\")\n", " qc.add_register(c1)\n", " for q, c in zip(ghz_qubits, c1):\n", " qc.measure(q, c)\n", " c2 = ClassicalRegister(len(checks), \"c2\")\n", " qc.add_register(c2)\n", " for q, c in zip(checks, c2):\n", " qc.measure(q, c)\n", " return qc\n", "\n", "\n", "def build_shadow_circuits(\n", " base_circuit: QuantumCircuit,\n", " n_measure: int,\n", " ghz_qubits: list,\n", " checks: list,\n", ") -> tuple:\n", " \"Build n_measure circuits by appending random Pauli basis rotations to each GHZ qubit.\"\n", " n = base_circuit.num_qubits\n", " all_bases = rng.integers(0, 3, size=(n_measure, n))\n", " circuits = []\n", " for m in range(n_measure):\n", " qc = base_circuit.copy()\n", " qc.barrier()\n", " for q in ghz_qubits:\n", " b = all_bases[m, q]\n", " if b == 1:\n", " qc.h(q)\n", " elif b == 2:\n", " qc.sdg(q)\n", " qc.h(q)\n", " add_measurements(qc, ghz_qubits, checks)\n", " circuits.append(qc)\n", " return circuits, all_bases\n", "\n", "\n", "def parse_hw_results_w_checks(\n", " result, bases: np.ndarray, ghz_qubits: list, checks: list\n", ") -> list:\n", " \"Parse hardware bitstrings into (ghz_bits, check_bits, bases) triples.\"\n", " snapshots = []\n", " n_ghz = len(ghz_qubits)\n", " n_checks = len(checks)\n", " for t, pub_result in enumerate(result):\n", " ghz_strings = (\n", " pub_result.data.c1.get_bitstrings()\n", " ) # one entry per shot\n", " chk_strings = pub_result.data.c2.get_bitstrings()\n", " for ghz_bitstring, check_bitstring in zip(ghz_strings, chk_strings):\n", " ghz_bits = [\n", " int(ghz_bitstring[n_ghz - 1 - i]) for i in range(n_ghz)\n", " ]\n", " check_bits = [\n", " int(check_bitstring[n_checks - 1 - i])\n", " for i in range(n_checks)\n", " ]\n", " snapshots.append((ghz_bits, check_bits, bases[t]))\n", " return snapshots\n", "\n", "\n", "def adapt_to_flat_bits(\n", " ghz_bits: list, ghz_qubits: list, n_total: int\n", ") -> list:\n", " \"Map GHZ-register bits back to a full n_total-length bit vector.\"\n", " bits = [0] * n_total\n", " for i, q in enumerate(ghz_qubits):\n", " bits[q] = ghz_bits[i]\n", " return bits\n", "\n", "\n", "def estimate_zz_nn_hw(\n", " snapshots: list, ghz_qubits: list, n_total: int, n_groups: int\n", ") -> np.ndarray:\n", " \"Estimate for each neighboring GHZ qubit pair from flat-bit snapshots.\"\n", " pairs = list(zip(ghz_qubits[:-1], ghz_qubits[1:]))\n", " all_est = np.column_stack(\n", " [\n", " snapshot_estimate(\n", " [\"Z\" if q in pair else \"I\" for q in range(n_total)], snapshots\n", " )\n", " for pair in pairs\n", " ]\n", " )\n", " return mom(all_est, n_groups).real" ] }, { "cell_type": "markdown", "id": "431a5bd2-e6ed-471b-ad9e-c4edd27784a8", "metadata": {}, "source": [ "## Small-scale simulator example\n", "\n", "First, prepare a 10-qubit GHZ state using statevector simulation and estimate its properties using local Pauli shadows. The GHZ state is a maximally entangled quantum state with three or more qubits serving as a benchmark for multipartite entanglement. It is a common building block for quantum communication and cryptography protocols. Additionally, the GHZ circuit has fully known analytical properties, making it ideal for verification in this tutorial.\n", "\n" ] }, { "attachments": {}, "cell_type": "markdown", "id": "988ee237", "metadata": {}, "source": [ "### Step 1: Map classical inputs to a quantum problem\n", "\n", "Prepare an ideal 10-qubit GHZ state.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "a3debf65-06df-4277-933e-14b6f6170756", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10-qubit GHZ state ready (norm = 1.000000)\n" ] } ], "source": [ "def ghz_statevector(num_qubits):\n", " qc = QuantumCircuit(num_qubits)\n", " qc.h(0)\n", " for q in range(num_qubits - 1):\n", " qc.cx(q, q + 1)\n", " return Statevector(qc)\n", "\n", "\n", "num_qubits = 10 # number of qubits\n", "psi = ghz_statevector(num_qubits)\n", "psi.seed(seed)\n", "print(\n", " f\"{num_qubits}-qubit GHZ state ready (norm = {np.linalg.norm(psi):.6f})\"\n", ")" ] }, { "cell_type": "markdown", "id": "ac6f36e3", "metadata": {}, "source": [ "### Step 2: Optimize for quantum hardware execution\n", "\n", "Next, define a function that collects shadow snapshots by rotating the GHZ state into random Pauli bases. This mirrors what the Sampler primitive does on hardware in the large-scale example.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "1834cb22", "metadata": {}, "outputs": [], "source": [ "def basis_circuit(b):\n", " \"\"\"0 -> Z basis (I), 1 -> X basis (H), 2 -> Y basis (Sdg;H).\"\"\"\n", " qc = QuantumCircuit(1)\n", " if b == 1:\n", " qc.h(0)\n", " elif b == 2:\n", " qc.sdg(0)\n", " qc.h(0)\n", " # b == 0: Z basis, no rotation\n", " return qc\n", "\n", "\n", "_basis = {\n", " b: basis_circuit(b) for b in (0, 1, 2)\n", "} # built once, shared with the estimator's _U\n", "\n", "\n", "def collect_shadow_snapshots(psi, n_snapshots, num_qubits):\n", " snapshots = []\n", " for _ in range(n_snapshots):\n", " bases = rng.integers(0, 3, size=num_qubits)\n", " rot = QuantumCircuit(num_qubits)\n", " for q, b in enumerate(bases):\n", " rot.compose(_basis[b], qubits=[q], inplace=True)\n", " evolved = psi.evolve(rot)\n", " evolved.seed(int(rng.integers(0, 2**31)))\n", " bitstr = next(\n", " iter(evolved.sample_counts(shots=1))\n", " ) # Born rule internal\n", " bits = [int(x) for x in bitstr[::-1]] # little-endian -> qubit\n", " snapshots.append((bits, bases))\n", " return snapshots" ] }, { "cell_type": "markdown", "id": "ad027ffc", "metadata": {}, "source": [ "### Step 3: Execute using Qiskit primitives\n", "\n", "For this simulation, use the statevector directly.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "4e2f32f4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Collected 5000 snapshots.\n" ] } ], "source": [ "n_snapshots = 5000 # number of shadow snapshots\n", "n_groups = 20 # groups/batches for median-of-means\n", "\n", "snapshots = collect_shadow_snapshots(psi, n_snapshots, num_qubits)\n", "print(f\"Collected {n_snapshots} snapshots.\")" ] }, { "attachments": {}, "cell_type": "markdown", "id": "b4d480b3", "metadata": {}, "source": [ "### Step 4: Post-process and return result in desired classical format\n", "\n", "Define the shadow estimator and the median-of-means (MOM) estimator for the random measurement bases used to collect the shadow snapshots. Use these functions to estimate properties of the GHZ state, including single-qubit magnetization, two-qubit correlators, and the two-qubit RDM, from the collected measurements. Then compare the estimates with the exact values.\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "b22a1b00", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Observable Estimated Ideal Max error\n", "---------------------------------------------------------\n", " (all qubits) - 0 0.0540\n", " mean 1.0040 1 0.1520\n", "2-RDM Frobenius err - 0 0.1080\n" ] } ], "source": [ "_U = {\n", " b: Operator(basis_circuit(b)) for b in (0, 1, 2)\n", "} # rotation Uᵢ per basis\n", "_P = {p: Operator(Pauli(p)) for p in \"IXYZ\"} # single-qubit Paulis\n", "\n", "\n", "def shadow_estimator(P, b, bit): # P is a label char, e.g. \"X\"\n", " UPU = _U[b] @ _P[P] @ _U[b].adjoint() # Uᵢ P Uᵢ†\n", " return (3 * UPU.data[bit, bit] - np.trace(_P[P].data)).real\n", "\n", "\n", "def mom(estimates: np.ndarray, n_groups: int) -> np.ndarray:\n", " \"Median-of-means outlier suppression.\"\n", " T = len(estimates)\n", " spg = T // n_groups # shots per group\n", " trimmed = estimates[\n", " : spg * n_groups\n", " ] # drop any leftover shots that don't cleanly fit into a group\n", " groups = trimmed.reshape(n_groups, spg, *estimates.shape[1:])\n", " means = groups.mean(axis=1) # take the mean of each group\n", " return np.median(means.real, axis=0) + 1j * np.median(\n", " means.imag, axis=0\n", " ) # take the median of means for real and imaginary parts separately\n", "\n", "\n", "def snapshot_estimate(\n", " pauli_string, snapshots\n", "): # full-width label, len == num_qubits\n", " out = np.empty(len(snapshots), dtype=complex)\n", " for t, (bits, bases) in enumerate(snapshots):\n", " e = 1.0 + 0j\n", " for q, P in enumerate(pauli_string):\n", " e *= shadow_estimator(P, bases[q], bits[q])\n", " out[t] = e\n", " return out\n", "\n", "\n", "Z_est = np.column_stack(\n", " [\n", " snapshot_estimate(\n", " [\"Z\" if q == i else \"I\" for q in range(num_qubits)], snapshots\n", " )\n", " for i in range(num_qubits)\n", " ]\n", ")\n", "Z_vals = mom(Z_est, n_groups).real\n", "\n", "ZZ_est = np.column_stack(\n", " [\n", " snapshot_estimate(\n", " [\"Z\" if q in (i, i + 1) else \"I\" for q in range(num_qubits)],\n", " snapshots,\n", " )\n", " for i in range(num_qubits - 1)\n", " ]\n", ")\n", "ZZ_nn = mom(ZZ_est, n_groups).real\n", "\n", "ZZ_full = np.zeros((n_snapshots, num_qubits, num_qubits), dtype=complex)\n", "for i in range(num_qubits):\n", " ZZ_full[:, i, i] = 1.0 # diagonal is always 1\n", " for j in range(i + 1, num_qubits):\n", " col = snapshot_estimate(\n", " [\"Z\" if q in (i, j) else \"I\" for q in range(num_qubits)],\n", " snapshots,\n", " )\n", " ZZ_full[:, i, j] = col\n", " ZZ_full[:, j, i] = col # symmetric\n", "ZZ_mat = mom(ZZ_full, n_groups).real\n", "\n", "# 2-qubit reduced density matrix (RDM) for pair (qi, qj)\n", "# the elementwise complex median (real/imag split) for the RDM coefficients generalizes the scalar median in Ref.[1]\n", "qi, qj = 0, num_qubits - 1\n", "terms = []\n", "for a in \"IXYZ\":\n", " for b in \"IXYZ\":\n", " label = \"\".join({qi: a, qj: b}.get(q, \"I\") for q in range(num_qubits))\n", " c = mom(snapshot_estimate(label, snapshots).reshape(-1, 1), n_groups)[\n", " 0\n", " ]\n", " terms.append((b + a, c / 4))\n", "rho_shadow = SparsePauliOp.from_list(terms).to_matrix()\n", "\n", "traced_out = [q for q in range(num_qubits) if q not in (qi, qj)]\n", "rho_exact = partial_trace(\n", " psi, traced_out\n", ").data # diag(0.5, 0, 0, 0.5) for GHZ\n", "\n", "print(f\"{'Observable':<22} {'Estimated':>12} {'Ideal':>8} {'Max error':>12}\")\n", "print(\"-\" * 57)\n", "print(\n", " f\"{' (all qubits)':<22} {'-':>12} {'0':>8} {np.max(np.abs(Z_vals)):>12.4f}\"\n", ")\n", "print(\n", " f\"{' mean':<22} {np.mean(ZZ_nn):>12.4f} {'1':>8} {np.max(np.abs(ZZ_nn - 1)):>12.4f}\"\n", ")\n", "print(\n", " f\"{'2-RDM Frobenius err':<22} {'-':>12} {'0':>8} {np.linalg.norm(rho_shadow - rho_exact):>12.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "1cff0147", "metadata": {}, "source": [ "#### Visualize the results\n", "\n", "The first plot shows the estimated single-qubit magnetization values for the different qubits. The second plot shows the distribution of estimated two-point correlator values for all qubit pairs. The final plot shows the magnitude of the two-qubit reduced density matrix entries. A fairly large collection of snapshots is used for the given system size, so the results are close to the exact ones.\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "31dc35ea-6554-4ca7-9c3b-0b5394c46e4e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[Text(0, 0, '00'), Text(0, 1, '01'), Text(0, 2, '10'), Text(0, 3, '11')]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(\n", " 1, 3, figsize=(14, 4), gridspec_kw={\"width_ratios\": [1, 1, 1]}\n", ")\n", "\n", "z_exact = psi.expectation_value(\n", " SparsePauliOp.from_sparse_list([(\"Z\", [0], 1.0)], num_qubits)\n", ").real # 0.0\n", "zz_exact = psi.expectation_value(\n", " SparsePauliOp.from_sparse_list([(\"ZZ\", [0, 1], 1.0)], num_qubits)\n", ").real # 1.0\n", "\n", "# \n", "axes[0].bar(range(num_qubits), Z_vals, color=\"#4a90d9\", alpha=0.85)\n", "axes[0].axhline(0, color=\"black\", lw=0.8)\n", "axes[0].set_ylim(-0.5, 0.5)\n", "axes[0].set_title(r\"$\\langle Z_i \\rangle$ (ideal = %0.2f)\" % z_exact)\n", "axes[0].set_xlabel(\"Qubit index\")\n", "axes[0].grid(axis=\"y\", alpha=0.3)\n", "\n", "# \n", "im = axes[1].imshow(ZZ_mat, vmin=-1, vmax=1, cmap=\"magma\")\n", "plt.colorbar(im, ax=axes[1])\n", "axes[1].set_title(r\"$\\langle Z_i Z_j \\rangle$ (ideal = %0.2f)\" % zz_exact)\n", "axes[1].set_xlabel(\"Qubit j\")\n", "axes[1].set_ylabel(\"Qubit i\")\n", "\n", "# 2-RDM magnitude\n", "im2 = axes[2].imshow(np.abs(rho_shadow), vmin=0, vmax=0.5, cmap=\"Blues\")\n", "plt.colorbar(im2, ax=axes[2])\n", "axes[2].set_title(r\"$|\\rho_2|$ for qubits (0, N-1)\")\n", "axes[2].set_xticks(range(4))\n", "axes[2].set_yticks(range(4))\n", "axes[2].set_xticklabels([\"00\", \"01\", \"10\", \"11\"])\n", "axes[2].set_yticklabels([\"00\", \"01\", \"10\", \"11\"])" ] }, { "cell_type": "markdown", "id": "7b12da12", "metadata": {}, "source": [ "#### Assess shot scaling\n", "\n", "In this tutorial, 5,000 snapshots are used to estimate the observables. The shot scaling of shadow tomography is assessed by plotting the mean error of one of the estimated functions against different shot budgets. The error in the two-qubit correlator decays as $1/\\sqrt{T}$ as expected for the standard statistical convergence of a mean estimator. The importance of shadow tomography is that the constant prefactor in front of $1/\\sqrt{T}$ does not grow with the number of qubits for local observables, such as those examined in this tutorial.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "id": "5da0364f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "T= 100 error = 0.1856\n", "T= 250 error = 0.1613\n", "T= 500 error = 0.1027\n", "T= 1000 error = 0.0827\n", "T= 2000 error = 0.0513\n", "T= 5000 error = 0.0311\n", "T=10000 error = 0.0162\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "shot_counts = [100, 250, 500, 1000, 2000, 5000, 10000]\n", "errors = []\n", "\n", "for N_test in shot_counts:\n", " snaps = collect_shadow_snapshots(psi, N_test, num_qubits)\n", " zz_vals = []\n", " for i in range(num_qubits - 1):\n", " ests = []\n", " for bits, bases in snaps:\n", " e = 1.0 + 0j\n", " for q in range(num_qubits):\n", " op = \"Z\" if q in (i, i + 1) else \"I\"\n", " e *= shadow_estimator(op, bases[q], bits[q])\n", " ests.append(e.real)\n", " zz_vals.append(np.mean(ests))\n", " errors.append(float(np.mean(np.abs(np.array(zz_vals) - 1.0))))\n", " print(f\"T={N_test:5d} error = {errors[-1]:.4f}\")\n", "\n", "fig, ax = plt.subplots(figsize=(6, 4))\n", "T_arr = np.array(shot_counts, dtype=float)\n", "ax.loglog(T_arr, errors, \"o-\", color=\"#4a90d9\", lw=2, label=\"Shadow estimate\")\n", "ax.loglog(\n", " T_arr,\n", " errors[0] * np.sqrt(shot_counts[0] / T_arr),\n", " \"k--\",\n", " alpha=0.5,\n", " label=r\"$1/\\sqrt{T}$ reference\",\n", ")\n", "ax.set_xlabel(\"Shots T\")\n", "ax.set_ylabel(r\"Mean $|\\langle Z_i Z_{i+1} \\rangle - 1|$\")\n", "ax.set_title(f\"Convergence for {num_qubits}-qubit GHZ\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0d6db390-e7a8-4efe-902c-8d9a312170c6", "metadata": {}, "source": [ "## Large-scale hardware example\n", "\n", "The preceding simulation confirms that local Pauli shadow tomography recovers the correct observables on a noiseless statevector. Running the same protocol using 100 qubits on real hardware introduces two challenges that do not appear in simulation: gate and readout errors accumulate across the deep GHZ circuit, and finding the optimal qubit layout requires searching the device's coupling graph rather than choosing a layout by hand.\n", "\n", "To address the first challenge, use low-overhead error detection with spacetime codes [\\[6\\]](#references). A small number of ancilla qubits adjacent to the GHZ tree are used as parity checks. Each ancilla is entangled with two of its GHZ neighbors before measurement, so a single bit-flip error on either neighbor flips the ancilla's outcome from `0` to `1`. Snapshots where any check ancilla reads `1` are discarded, leaving a cleaner subset from which to estimate $\\langle Z_i Z_{i+1}\\rangle$.\n", "\n", "The second challenge is handled by searching over candidate GHZ tree layouts and selecting the one that maximizes parity-check coverage, that is, the fraction of circuit wires monitored by at least one check, within an acceptable increase in circuit depth.\n", "\n" ] }, { "cell_type": "markdown", "id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25", "metadata": {}, "source": [ "### Steps 1-4 (performed in a single code block)\n", "\n", "The same four-step workflow is now applied to a 100-qubit GHZ state on real quantum hardware, adding two hardware-specific features that do not arise in simulation:\n", "\n", "* **Step 1** characterizes the backend noise (gate errors, readout errors, and T2 times) so that unreliable qubits and edges are excluded from the circuit layout.\n", "* **Step 2** searches over possible GHZ tree layouts to find one that maximizes parity-check coverage within an acceptable circuit depth — replacing the trivial chain used in simulation.\n", "* **Step 3** builds the shadow circuits (one per random Pauli basis), transpiles them to hardware-native gates, and submits them to the Sampler primitive with `shots=8` (each circuit produces eight independent bitstrings).\n", "* **Step 4** postselects snapshots whose check qubits all read 0 (no detected error), then estimates $\\langle Z_i Z_{i+1}\\rangle$ with and without that postselection to demonstrate the improvement from error detection.\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "502e92dd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GHZ qubits: [57, 47, 67, 46, 66, 48, 45, 65, 49, 44, 64, 38, 43, 77, 29, 42, 85, 28, 56, 41, 84, 30, 27, 36, 86, 31, 21, 87, 32, 20, 88, 33, 22, 97, 89, 34, 23, 107, 78, 39, 35, 16, 106, 90, 69, 53, 19, 24, 3, 108, 105, 91, 70, 52, 15, 25, 2, 109, 104, 92, 71, 54, 51, 14, 4, 1, 110, 103, 98, 93, 58, 55, 13, 5, 0, 118, 96, 72, 59, 12, 6, 129, 102, 73, 11, 7, 128, 101, 74, 10, 8, 127, 100, 9, 126, 116, 125, 121, 124, 120] (100)\n", "Check qubits: [68, 37, 50, 63, 83, 17, 26, 18, 117, 111, 75] (11)\n", "Covered fraction (no idle): 0.5454545454545454\n", "Shadow snapshots collected: 8000\n", "Shadow snapshots kept post error-detection: 566\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Before detection: = 0.573749999999951\n", "After detection: = 0.8817377312951776\n" ] } ], "source": [ "# ── Step 1: Map classical inputs to a quantum problem ─────────────────────────\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=133\n", ")\n", "coupling_map = backend.target.build_coupling_map()\n", "n_total = backend.num_qubits\n", "\n", "thresh_cz, thresh_meas, thresh_t2 = 0.025, 0.15, 10\n", "bad_edges = bad_cz(backend.target, threshold=thresh_cz)\n", "bad_nodes_readout = bad_readout(backend.target, threshold=thresh_meas)\n", "dead_qubits = bad_readout(backend.target, threshold=0.4)\n", "bad_nodes = list(\n", " set(bad_nodes_readout)\n", " | set(bad_coherence(backend.target, threshold=thresh_t2))\n", ")\n", "\n", "# ── Step 2: Optimize for quantum hardware execution ───────────────────────────\n", "# Search over starting qubits and random node-skip sets to find the GHZ layout\n", "# with the best parity-check coverage within an acceptable circuit depth.\n", "ghz_size = 100\n", "max_skips, shuffles, max_depth_increase = 10, 200, 10\n", "w_idle, w_gate = 0.2, 0.8\n", "\n", "# Find the root that yields the shallowest GHZ — exclude bad_nodes from the start\n", "best_root = None\n", "base_depth = float(\"inf\")\n", "for root in range(backend.num_qubits):\n", " if root in bad_nodes:\n", " continue\n", " qc, ghz_qubits, _ = parallel_ghz(\n", " root, ghz_size, backend, bad_edges, bad_nodes\n", " )\n", " if len(ghz_qubits) != ghz_size:\n", " continue\n", " depth = qc.depth(lambda x: x.operation.num_qubits == 2)\n", " if depth < base_depth:\n", " best_root, base_depth = root, depth\n", "if best_root is None:\n", " raise Exception(\n", " \"No qualifying GHZ root layout found. Relax error thresholds.\"\n", " )\n", "\n", "qc, ghz_qubits, layers = parallel_ghz(\n", " best_root, ghz_size, backend, bad_edges, bad_nodes\n", ")\n", "\n", "degree_two_nodes = [\n", " i\n", " for i in ghz_qubits\n", " if all(n in ghz_qubits for n in coupling_map.neighbors(i))\n", " and len(coupling_map.neighbors(i)) >= 2\n", "]\n", "\n", "best_covered_fraction = -1\n", "best_qc = best_checks = best_parities = best_layers = best_ghz_qubits = None\n", "\n", "for num_skips in range(max_skips):\n", " for _ in range(shuffles):\n", " skip_trial = list(\n", " rng.choice(degree_two_nodes, num_skips, replace=False)\n", " )\n", " qc, ghz_qubits, layers = parallel_ghz(\n", " best_root, ghz_size, backend, bad_edges, skip_trial + bad_nodes\n", " )\n", " if len(ghz_qubits) != ghz_size:\n", " continue\n", " depth = qc.depth(lambda x: x.operation.num_qubits == 2)\n", " checks, parities = [], []\n", " for ancilla, neighbors in get_low_overhead_ancillas(\n", " coupling_map, ghz_qubits\n", " ).items():\n", " if ancilla in dead_qubits:\n", " continue\n", " good_neighbors = [\n", " n\n", " for n in neighbors\n", " if [ancilla, n] not in bad_edges\n", " and [n, ancilla] not in bad_edges\n", " ]\n", " if len(good_neighbors) < 2:\n", " continue\n", " checks.append(ancilla)\n", " parities.append((good_neighbors[0], good_neighbors[1]))\n", " qc.cx(good_neighbors[0], ancilla)\n", " qc.cx(good_neighbors[1], ancilla)\n", " covered_fraction = weighted_coverage(\n", " layers=layers, parities=parities, w_idle=w_idle, w_gate=w_gate\n", " )\n", " if (\n", " covered_fraction > best_covered_fraction\n", " and depth <= base_depth + max_depth_increase\n", " ):\n", " best_covered_fraction = covered_fraction\n", " (\n", " best_qc,\n", " best_ghz_qubits,\n", " best_checks,\n", " best_parities,\n", " best_layers,\n", " ) = qc, ghz_qubits, checks, parities, layers\n", "\n", "qc, ghz_qubits, checks, parities, layers = (\n", " best_qc,\n", " best_ghz_qubits,\n", " best_checks,\n", " best_parities,\n", " best_layers,\n", ")\n", "if len(ghz_qubits) != ghz_size:\n", " raise Exception(\"No GHZ found. Relax error thresholds.\")\n", "\n", "print(f\"GHZ qubits: {ghz_qubits} ({len(ghz_qubits)})\")\n", "print(f\"Check qubits: {checks} ({len(checks)})\")\n", "print(\n", " \"Covered fraction (no idle): \",\n", " weighted_coverage(\n", " layers=layers, parities=parities, w_idle=0.0, w_gate=1.0\n", " ),\n", ")\n", "\n", "# ── Step 3: Execute using Qiskit primitives ───────────────────────────────────\n", "# Build one shadow circuit per random Pauli basis, each repeated for shots_per_basis\n", "# independent bitstrings. Transpile to hardware-native gates, then submit using the\n", "# Sampler primitive. Each circuit yields shots_per_basis snapshots in that basis.\n", "n_hw = 1000\n", "shots_per_basis = 8 # multiple shots per basis circuit\n", "shadow_circs, hw_bases = build_shadow_circuits(qc, n_hw, ghz_qubits, checks)\n", "sampler = Sampler(backend)\n", "sampler.options.environment.job_tags = [\"TUT_STGHZS\"]\n", "pm = generate_preset_pass_manager(optimization_level=1, backend=backend)\n", "isa_circuit = pm.run(shadow_circs)\n", "job = sampler.run(isa_circuit, shots=shots_per_basis)\n", "result = job.result()\n", "\n", "# ── Step 4: Post-process and return result in desired classical format ─────────\n", "# Filter snapshots whose check qubits all read 0 (no detected error), then\n", "# estimate with and without the error-detection postselection.\n", "snapshots_hw = parse_hw_results_w_checks(result, hw_bases, ghz_qubits, checks)\n", "snapshots_kept = [\n", " (ghz_bits, check_bits, bases)\n", " for (ghz_bits, check_bits, bases) in snapshots_hw\n", " if all(b == 0 for b in check_bits)\n", "]\n", "snapshots_kept_adapted = [\n", " (adapt_to_flat_bits(ghz_bits, ghz_qubits, n_total), bases)\n", " for (ghz_bits, check_bits, bases) in snapshots_kept\n", "]\n", "snapshots_hw_adapted = [\n", " (adapt_to_flat_bits(ghz_bits, ghz_qubits, n_total), bases)\n", " for (ghz_bits, check_bits, bases) in snapshots_hw\n", "]\n", "\n", "print(\"Shadow snapshots collected:\", len(snapshots_hw))\n", "print(\"Shadow snapshots kept post error-detection:\", len(snapshots_kept))\n", "\n", "n_groups = min(5, len(snapshots_kept)) # groups for median-of-means\n", "ZZ_hw = estimate_zz_nn_hw(snapshots_hw_adapted, ghz_qubits, n_total, n_groups)\n", "ZZ_kept = estimate_zz_nn_hw(\n", " snapshots_kept_adapted, ghz_qubits, n_total, n_groups\n", ")\n", "\n", "# Plot before and after error detection\n", "fig, axes = plt.subplots(1, 2, figsize=(13, 4))\n", "\n", "axes[0].plot(ZZ_hw, \"o-\", color=\"#e05c5c\", markersize=4)\n", "axes[0].axhline(1.0, color=\"gray\", ls=\"--\", alpha=0.7, label=\"Ideal\")\n", "axes[0].set_ylim(-0.6, 3.0)\n", "axes[0].set_title(r\"$\\langle Z_i Z_{i+1} \\rangle$ (without error detection)\")\n", "axes[0].set_xlabel(\"GHZ qubit pair\")\n", "axes[0].legend()\n", "axes[0].grid(alpha=0.3)\n", "\n", "axes[1].plot(ZZ_kept, \"o-\", color=\"#e05c5c\", markersize=4)\n", "axes[1].axhline(1.0, color=\"gray\", ls=\"--\", alpha=0.7, label=\"Ideal\")\n", "axes[1].set_ylim(-0.6, 3.0)\n", "axes[1].set_title(r\"$\\langle Z_i Z_{i+1} \\rangle$ (with error detection)\")\n", "axes[1].set_xlabel(\"GHZ qubit pair\")\n", "axes[1].legend()\n", "axes[1].grid(alpha=0.3)\n", "\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(\"Before detection: = \", float(np.average(ZZ_hw)))\n", "print(\"After detection: = \", float(np.average(ZZ_kept)))" ] }, { "cell_type": "markdown", "id": "dde3cae9", "metadata": {}, "source": [ "#### Interpret the results\n", "\n", "The two plots show $\\langle Z_i Z_{i+1}\\rangle$ estimated from shadow tomography for each neighboring qubit pair in the 100-qubit GHZ chain, before and after postselecting on the parity checks.\n", "\n", "For an ideal GHZ state, all nearest-neighbor correlators equal $1$. On hardware, noise depresses the raw estimates and introduces qubit-to-qubit scatter. After error detection, snapshots in which at least one detectable error occurred are removed, and the remaining estimates are pulled closer to $1$ with reduced variance, which is consistent with a cleaner effective state.\n", "\n", "The improvement after error detection is consistent with the checks detecting real errors.\n", "\n" ] }, { "cell_type": "markdown", "id": "75f48e6a-c7e4-46f3-9d39-a7a877427a04", "metadata": {}, "source": [ "## Next steps\n", "\n", "Explore the following resources to learn more:\n", "\n", "\n", " * [General formulation of quantum information](/learning/courses/general-formulation-of-quantum-information)\n", " * [Low-overhead error detection using spacetime codes tutorial](/docs/tutorials/spacetime-codes)\n", " * [Qiskit Paulice add-on API documentation](/docs/api/qiskit-addon-paulice)\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "9d876315", "metadata": {}, "source": [ "## References\n", "\n", "* \\[1] Huang, H. Y., Kueng, R., & Preskill, J. (2020). [Predicting many properties of a quantum system from very few measurements](https://www.nature.com/articles/s41567-020-0932-7). *Nature Physics*, 16, 1050–1057.\n", "* \\[2] Zhao, A., Rubin, N. C., & Miyake, A. (2021). [Fermionic partial tomography via classical shadows](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.127.110504). *Phys. Rev. Lett.*, 127, 110504.\n", "* \\[3] Low, G. H. (2022). [Classical shadows of fermions with particle number symmetry](https://arxiv.org/abs/2208.08964). *arXiv:2208.08964*.\n", "* \\[4] Bertoni, C., Haferkamp, J., Hinsche, M. et al. (2024). [Shallow Shadows: Expectation Estimation Using Low-Depth Random Clifford Circuits](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.133.020602). *Phys. Rev. Lett.*, 133, 020602.\n", "* \\[5] Hu, H.-Y., Gu, A., Majumder, S. et al. (2025). [Demonstration of robust and efficient quantum property learning with shallow shadows](https://doi.org/10.1038/s41467-025-57349-w). *Nature Communications*, 16, 2943.\n", "* \\[6] Martiel, S. & Javadi-Abhari, A. (2025). [Low-overhead error detection with spacetime codes](https://arxiv.org/abs/2504.15725). *arXiv:2504.15725*.\n", "\n" ] }, { "cell_type": "markdown", "id": "d74f7bbb", "metadata": {}, "source": [ "© IBM Corp. 2026\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "id": "a1b8767d", "source": "© IBM Corp., 2017-2026" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3" } }, "nbformat": 4, "nbformat_minor": 5 }