{
"cells": [
{
"cell_type": "markdown",
"id": "d2c31ae8",
"metadata": {},
"source": [
"---\n",
"title: \"Tomografia d'ombra dello stato di Greenberger-Horne-Zeilinger\"\n",
"description: \"Questo tutorial illustra il protocollo di tomografia d’ombra attraverso la stima delle grandezze osservabili dello stato di Greenberger-Horne-Zeilinger con un massimo di 100 qubit\"\n",
"---\n",
"\n",
"{/* cspell:ignore matchgate binom bitstr elementwise correlator vmin vmax ests postselects Kueng Zhao Miyake Bertoni Haferkamp Hinsche */}\n",
"\n",
"\n",
"\n",
"# Tomografia d'ombra dello stato di Greenberger-Horne-Zeilinger\n",
"\n",
"Stima dei tempi di *esecuzione: meno di 10 secondi su un processore Heron r2 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8bf80006",
"metadata": {},
"source": [
"\n",
"\n",
"## Risultati di apprendimento\n",
"\n",
"* In che modo la tomografia d’ombra ricostruisce un’ombra classica di uno stato quantistico e perché richiede un numero esponenzialmente inferiore di misurazioni rispetto alla tomografia completa dello stato\n",
"* Come implementare la tomografia locale dell’ombra di Pauli utilizzando Qiskit per stimare le magnetizzazioni di singoli qubit, i correlatori a due punti e le matrici di densità ridotta a due qubit\n",
"* Come utilizzare [qiskit-paulice](/docs/addons/qiskit-addon-paulice) per il rilevamento degli errori a basso overhead tramite codici spazio-temporali in un esperimento hardware su larga scala\n",
"\n",
"\n",
"\n",
"## Prerequisiti\n",
"\n",
"* [Misurazioni](/learning/courses/general-formulation-of-quantum-information/general-measurements/introduction)\n",
"* [Tomografia dello stato quantistico](/learning/courses/general-formulation-of-quantum-information/general-measurements/discrimination-and-tomography)\n",
"* [Rilevamento degli errori mediante codici spazio-temporali](/docs/tutorials/spacetime-codes)\n",
"\n",
"\n",
"\n",
"## Sfondo\n",
"\n",
"La tomografia completa dello stato quantistico di un sistema a $N$ qubit richiede un numero esponenziale di misurazioni per stimare tutti i parametri indipendenti $4^N -1$ della matrice di densità. Per un sistema di 50 qubit, ad esempio, questo processo richiede circa $10^{30}$, ovvero più di un miliardo di volte il numero di stelle presenti nell’universo osservabile. La tomografia d’ombra [\\[1\\]](#references) è un protocollo di informazione quantistica che costruisce un’“ombra classica” dello stato quantistico, dalla quale è possibile stimare le osservabili di interesse utilizzando un numero di misurazioni che varia in modo logaritmico rispetto al numero di osservabili ed è indipendente dalle dimensioni del sistema $N$; ciò la rende rilevante in molti ambiti dell’informatica quantistica, tra cui la stima dell’energia, il calcolo dell’entropia di entanglement e delle funzioni di correlazione, la stima dei kernel nell’apprendimento automatico quantistico e altri ancora.\n",
"\n",
"Questo tutorial illustra come utilizzare la tomografia d'ombra per stimare grandezze osservabili locali, quali l' $\\langle Z_i \\rangle$ e di magnetizzazione di un singolo qubit per il qubit $i$, la funzione di correlazione a due punti $\\langle Z_i Z_j\\rangle$ e la matrice di densità ridotta a due qubit (2-RDM), dello stato di Greenberger-Horne-Zeilinger (GHZ) altamente intrecciato. Il tutorial illustra innanzitutto un flusso di lavoro di tomografia d'ombra su piccola scala su uno stato GHZ a 10 qubit, utilizzando la simulazione con vettori di stato. Viene poi mostrata una dimostrazione con 100 qubit su hardware quantistico reale.\n",
"\n",
"\n",
"\n",
"### Protocollo di tomografia delle ombre\n",
"\n",
"Il protocollo di tomografia dell'ombra prevede le seguenti fasi:\n",
"\n",
"1. Applicare allo stato quantistico iniziale un’ $U$ e unitaria casuale, estratta da un insieme unitario.\n",
"2. Misurare il risultato nella base computazionale, ottenendo una stringa di bit $\\ket{b}$.\n",
"3. Salvare lo snapshot classico, $\\hat{\\rho} = \\mathcal{M}^{-1}(U^{\\dagger}\\ket{b}\\bra{b}U)$, dove $\\mathcal{M}^{-1}$ è l'inverso del canale di misura $\\mathcal{M}(\\rho)=\\mathbb{E}[U^{\\dagger}\\ket{b}\\bra{b}U]$.\n",
"4. Ripetere il processo $T$ volte, dove l’insieme delle istantanee $T$ $\\{\\hat{\\rho}_t\\}$ forma una classica “ombra” $S(\\rho, T)$.\n",
"5. Stima l' $o$ e osservabile dell'operatore $O$ come $\\frac{1}{T}\\sum_{t=1}^{T}\\text{Tr}(O\\hat{\\rho}_t)$.\n",
"\n",
"\n",
"\n",
"### Varianti della tomografia delle ombre\n",
"\n",
"La scelta dell'insieme unitario, ovvero della base, nella fase 1 del protocollo, influisce sulla complessità del campione, sulla varianza delle stime e sulla completezza delle informazioni di una data tomografia ombra. Nella letteratura sono stati presi in esame numerosi gruppi di modelli, ciascuno caratterizzato da un compromesso tra complessità ed efficienza nella stima delle variabili osservabili. Le proprietà delle varianti più comunemente utilizzate sono riassunte nella tabella seguente, sebbene tale elenco non sia esaustivo.\n",
"\n",
"| Variante | Insieme | Complessità di campionamento | Profondità del circuito | Grandezze osservabili idonee |\n",
"| --------------------------------------------------------------- | ------------------------------------ | ------------------------------------------------------ | --------------------------- | --------------------------------------- |\n",
"| Pauli locale [\\[1\\]](#references) | Paulis a singolo qubit $\\{X, Y, Z\\}$ | $\\mathcal{O}(3^k\\text{log}M/\\varepsilon^2)$ | 1 | locale/a pochi corpi |\n",
"| Clifford globale [\\[1\\]](#references) | $N$ -qubit di Clifford | $\\mathcal{O}(\\text{Tr}(O^2)\\text{log}M/\\varepsilon^2)$ | $\\mathcal{O}(N\\text{log}N)$ | basso rango, fedeltà |\n",
"| Gaussiano fermionico [\\[2\\]](#references), [\\[3\\]](#references) | circuiti matchgate | $\\mathcal{O}(\\binom{\\mu}{k}\\frac{1}{\\varepsilon^2})$ | $\\mathcal{O}(N^2)$ | fermionico |\n",
"| ombre superficiali [\\[4\\]](#references), [\\[5\\]](#references) | circuiti a bassa profondità | $\\mathcal{O}(k2^k\\text{log}M/\\varepsilon^2)$ | $\\mathcal{O}(\\text{log}N)$ | non locale (compatibile con l'hardware) |\n",
"\n",
"I termini utilizzati nella tabella sono definiti come segue:\n",
"\n",
"* Per \"complessità del campione\" si intende il numero di misurazioni necessarie per una determinata variante della tomografia d'ombra.\n",
"* Il termine \"profondità del circuito\" descrive quanto sia oneroso implementare le misurazioni sull'hardware quantistico.\n",
"* $N$ è il numero di qubit.\n",
"* $k$ è l'indirizzo della risorsa osservabile (ad esempio, $\\langle Z_i \\rangle$ ha $k=1$ e $\\langle Z_i Z_j\\rangle$ ha $k=2$ ).\n",
"* $M$ è il numero di grandezze osservabili da stimare contemporaneamente.\n",
"* $\\varepsilon$ è la tolleranza sull'errore additivo.\n",
"* $\\mu$ è il numero di fermioni.\n",
"\n",
"Sebbene l’insieme locale di Pauli sia quello che richiede meno risorse hardware, comporta una penalizzazione in termini di complessità di campion $3^k$, il che lo rende poco adatto alla stima delle proprietà non locali. Più in generale, il limite di complessità campionaria dell’ombra locale varia proporzionalmente a $4^k$, migliorando fino a $3^k$ nel caso particolare degli osservabili a prodotto tensoriale.\n",
"\n",
"Al contrario, sebbene l’insieme di Clifford globale possa essere utilizzato per stimare osservabili non locali, la sua profondità di circuito non è favorevole. Le ombre gaussiane fermioniche sono utili nelle simulazioni della fisica dei sistemi a molti corpi, della chimica e dei materiali, ma richiedono circuiti matchgate, che potrebbero comportare un sovraccarico significativo sui dispositivi generici basati su gate.\n",
"\n",
"\"Shallow shadows\" presenta una profondità di circuito compresa tra quella locale di Pauli e quella globale di Clifford, il che lo rende una variante molto promettente per l'implementazione pratica, sebbene permangano ancora alcune questioni aperte riguardo alla profondità ottimale del circuito per un dato osservabile.\n",
"\n",
"Man mano che i confini della ricerca sulla tomografia d’ombra continuano ad ampliarsi, alcuni dei problemi più rilevanti da affrontare includono la scelta ottimale dell’insieme di misurazioni per un insieme target di osservabili, il miglioramento dei limiti di complessità del campione e l’estensione del protocollo agli osservabili che evolvono nel tempo. Ciononostante, la tomografia d’ombra rappresenta ancora uno dei principali flussi di lavoro per estrarre il massimo delle informazioni classiche da uno stato quantistico utilizzando un numero minimo di misurazioni, il che la rende ideale per caratterizzare stati quantistici su larga scala, laddove sia la tomografia completa che la simulazione classica risultano impossibili.\n",
"\n",
"\n",
"\n",
"### L'ensemble locale Pauli\n",
"\n",
"Questo tutorial utilizza la base casuale di Pauli, ovvero porte di Clifford su un singolo qubit, per raccogliere ombre e stimare le proprietà dello stato GHZ. Si noti che le misure di Pauli rappresentano un'ottima scelta per stimare gli osservabili locali (come la densità di spin, le funzioni di correlazione locali e la densità di probabilità di posizione), ma non quelli globali (come l'energia totale, la quantità di moto e la carica).\n",
"\n",
"In questo caso, per ogni qubit viene scelta in modo indipendente una base di Pauli casuale $\\{X, Y, Z\\}$ e viene applicata la corrispondente rotazione del singolo qubit, seguita dalla misurazione:\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",
"Per questo insieme, l'inverso del canale di misura è:\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",
"Il fattore $3$ deriva dal calcolo della media sulle tre basi di Pauli, che produce un canale depolarizzante con parametro $\\frac{1}{3}$. Lo stimatore completo si scompone quindi in un prodotto di stimatori a singolo qubit. Come illustrato nella tabella, per $M$ osservabili, ciascuno di località $k$, il numero di ombre di Pauli necessarie per stimare tutti gli osservabili con un errore additivo $\\varepsilon$ con elevata probabilità è pari a $\\mathcal{O}(\\frac{3^k\\text{log}(M)}{\\varepsilon^2})$. La dimensione del sistema non influisce sul costo della misurazione fintanto che gli osservabili rimangono locali; pertanto, è possibile stimare simultaneamente un numero esponenziale di osservabili con un costo solo logaritmico.\n",
"\n",
"\n",
"\n",
"### Mediana delle medie\n",
"\n",
"Una media semplice calcolata su una serie di istantanee dell’ $T$ e riduce la varianza di una stima a campionamento singolo di una grandezza osservabile secondo la formula $\\frac{1}{T}$, ma è sensibile ai valori anomali. Per ovviare a questo problema si può ricorrere alla mediana delle medie. Il processo prevede quanto segue:\n",
"\n",
"* Suddivisione delle istantanee di “ $T$ ” in $K$ gruppi di uguale dimensione\n",
"* Calcolo della media all’interno di ciascun gruppo\n",
"* Calcolo della mediana delle medie dei gruppi\n",
"\n",
"Questo tutorial illustra l'applicazione della tomografia locale dell'ombra di Pauli per stimare le proprietà dello stato GHZ, utilizzando innanzitutto un esempio simulato a 10 qubit e dimostrando poi un esperimento su hardware a 100 qubit.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "55b94021",
"metadata": {},
"source": [
"\n",
"\n",
"## Requisiti\n",
"\n",
"Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:\n",
"\n",
"* Qiskit SDK v2.0 o versioni successive, con supporto [alla visualizzazione](/docs/api/qiskit/visualization)\n",
"* IBM Quantum Compute v0.28 o versioni successive (`pip install qiskit-ibm-runtime`)\n",
"* Qiskit Paulice (`pip install qiskit-paulice`)\n",
"\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "7db2e559",
"metadata": {},
"source": [
"\n",
"\n",
"## Configura\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": [
"Le seguenti funzioni di supporto vengono utilizzate nell'esempio relativo all'hardware su larga scala presentato più avanti nel tutorial:\n",
"\n",
"* `bad_cz`, `bad_readout`, `bad_coherence` — Interrogare i dati di calibrazione del backend per identificare i gate, i qubit e i bordi che superano le soglie di errore o di coerenza, in modo da poterli escludere dal layout del circuito.\n",
"* `parallel_ghz` — Costruisce un circuito GHZ di una determinata dimensione utilizzando una ricerca in larghezza sulla mappa di accoppiamento, distribuendo gli strati CNOT in parallelo per ridurre al minimo la profondità del circuito.\n",
"* `weighted_coverage`, `active_wires`, `z_trace_backward` — Calcolare quale porzione della superficie di errore del circuito GHZ sia protetta da un dato insieme di controlli di parità spazio-temporali, utilizzando la tecnica della propagazione all’indietro descritta in [\\[6\\]](#references).\n",
"* `add_measurements` — Inserisce le misurazioni dei qubit GHZ e dei qubit di controllo in registri classici separati, in modo che possano essere post-selezionate in modo indipendente.\n",
"* `build_shadow_circuits` — Genera un circuito per ogni base di Pauli casuale aggiungendo le rotazioni appropriate di un singolo qubit prima della misurazione.\n",
"* `parse_hw_results_w_checks` — Suddivide l'output del Sampler in triplette `(ghz_bits, check_bits, bases)` , una per ogni shot.\n",
"* `adapt_to_flat_bits`, `estimate_zz_nn_hw` — Mappare i bit del registro GHZ nello spazio indice completo del dispositivo e calcolare la stima della mediana delle medie ( $\\langle Z_i Z_{i+1}\\rangle$ ).\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": [
"\n",
"\n",
"## Esempio di simulatore su piccola scala\n",
"\n",
"Per prima cosa, preparare uno stato GHZ a 10 qubit utilizzando la simulazione con vettori di stato e stimarne le proprietà ricorrendo alle ombre di Pauli locali. Lo stato GHZ è uno stato quantistico con entanglement massimo, costituito da tre o più qubit, che funge da punto di riferimento per l'entanglement multipartito. Si tratta di un elemento fondamentale per i protocolli di comunicazione quantistica e di crittografia. Inoltre, il circuito GHZ presenta proprietà analitiche ben note, il che lo rende ideale per la verifica in questo tutorial.\n",
"\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "988ee237",
"metadata": {},
"source": [
"\n",
"\n",
"### Fase 1: Mappare gli input classici su un problema quantistico\n",
"\n",
"Preparare uno stato GHZ ideale a 10 qubit.\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": [
"\n",
"\n",
"### Fase 2: Ottimizzazione per l'esecuzione su hardware quantistico\n",
"\n",
"Successivamente, definire una funzione che raccolga le istantanee d’ombra ruotando lo stato GHZ in basi di Pauli casuali. Questo rispecchia il funzionamento della primitiva Sampler sull’hardware nell’esempio su larga scala.\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": [
"\n",
"\n",
"### Fase 3: Eseguire il comando utilizzando Qiskit primitives\n",
"\n",
"Per questa simulazione, utilizzare direttamente il vettore di stato.\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": [
"\n",
"\n",
"### Fase 4: Elaborazione finale e restituzione del risultato nel formato classico desiderato\n",
"\n",
"Definire lo stimatore “shadow” e lo stimatore della mediana delle medie (MOM) per le basi di misurazione casuali utilizzate per acquisire le istantanee “shadow”. Utilizza queste funzioni per stimare le proprietà dello stato GHZ, tra cui la magnetizzazione di un singolo qubit, i correlatori a due qubit e l'RDM a due qubit, sulla base delle misurazioni raccolte. Quindi confronta le stime con i valori esatti.\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": [
"\n",
"\n",
"#### Visualizza i risultati\n",
"\n",
"Il primo grafico mostra i valori stimati della magnetizzazione dei singoli qubit per i diversi qubit. Il secondo grafico mostra la distribuzione dei valori stimati del correlatore a due punti per tutte le coppie di qubit. Il grafico finale mostra l'ampiezza delle voci della matrice di densità ridotta a due qubit. Per le dimensioni del sistema in questione viene utilizzata una raccolta piuttosto ampia di istantanee, pertanto i risultati sono vicini a quelli esatti.\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": [
""
]
},
"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": [
"\n",
"\n",
"#### Valutare il ridimensionamento delle immagini\n",
"\n",
"In questo tutorial vengono utilizzate 5.000 istantanee per stimare le grandezze osservabili. Il ridimensionamento degli shot nella tomografia delle ombre viene valutato tracciando l'errore medio di una delle funzioni stimate in funzione dei diversi budget di shot. L'errore nel correlatore a due qubit decresce secondo la legge $1/\\sqrt{T}$, come previsto per la convergenza statistica standard di uno stimatore della media. L'importanza della tomografia d'ombra risiede nel fatto che il prefattore costante che precede l' $1/\\sqrt{T}$ e non cresce con il numero di qubit per gli osservabili locali, come quelli esaminati in questo 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": [
""
]
},
"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": [
"\n",
"\n",
"## Esempio di hardware su larga scala\n",
"\n",
"La simulazione precedente conferma che la tomografia locale dell’ombra di Pauli ricostruisce le osservabili corrette su un vettore di stato privo di rumore. L'esecuzione dello stesso protocollo utilizzando 100 qubit su hardware reale comporta due difficoltà che non si presentano nella simulazione: gli errori di gate e di lettura si accumulano lungo l'intero circuito GHZ, e per individuare la disposizione ottimale dei qubit è necessario effettuare una ricerca nel grafico di accoppiamento del dispositivo, anziché scegliere manualmente una disposizione.\n",
"\n",
"Per affrontare la prima sfida, utilizzare un sistema di rilevamento degli errori a basso overhead con codici spazio-temporali [\\[6\\]](#references). Un numero limitato di qubit ancilla adiacenti all’albero GHZ viene utilizzato per i controlli di parità. Ogni ancilla è intrecciata con due dei suoi vicini GHZ prima della misurazione, quindi un singolo errore di inversione di bit su uno dei due vicini inverte il risultato dell'ancilla da `0` a `1`. Le istantanee in cui una qualsiasi lettura ancilla di controllo viene `1` scartata, lasciando un sottoinsieme più pulito da cui stimare l' $\\langle Z_i Z_{i+1}\\rangle$ e.\n",
"\n",
"La seconda sfida viene affrontata effettuando una ricerca tra le possibili configurazioni dell’albero GHZ e selezionando quella che massimizza la copertura del controllo di parità, ovvero la frazione di linee del circuito monitorate da almeno un controllo, entro un aumento accettabile della profondità del circuito.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25",
"metadata": {},
"source": [
"\n",
"\n",
"### Passaggi da 1 a 4 (eseguiti in un unico blocco di codice)\n",
"\n",
"Lo stesso flusso di lavoro in quattro fasi viene ora applicato a uno stato GHZ da 100 qubit su hardware quantistico reale, aggiungendo due caratteristiche specifiche dell’hardware che non si riscontrano nella simulazione:\n",
"\n",
"* **La fase 1** caratterizza il rumore del backend (errori di gate, errori di lettura e tempi di \" T2 \") in modo da escludere i qubit e i fronti inaffidabili dal layout del circuito.\n",
"* **Il passaggio 2** analizza le possibili configurazioni dell'albero GHZ per individuare quella che massimizza la copertura del controllo di parità entro una profondità del circuito accettabile, sostituendo la catena banale utilizzata nella simulazione.\n",
"* **Il passaggio 3** costruisce i circuiti ombra (uno per ogni base di Pauli casuale), li traspone in porte native dell'hardware e li invia alla primitiva Sampler con `shots=8` (ogni circuito produce otto stringhe di bit indipendenti).\n",
"* **Il passaggio 4** effettua una post-selezione delle istantanee in cui tutti i qubit di controllo danno valore 0 (nessun errore rilevato), quindi stima l' $\\langle Z_i Z_{i+1}\\rangle$ e con e senza tale post-selezione per dimostrare il miglioramento derivante dal rilevamento degli errori.\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": [
""
]
},
"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": [
"\n",
"\n",
"#### Interpretare i risultati\n",
"\n",
"I due grafici mostrano l' $\\langle Z_i Z_{i+1}\\rangle$ e stimata tramite tomografia d'ombra per ciascuna coppia di qubit adiacenti nella catena GHZ da 100 qubit, prima e dopo la post-selezione sui controlli di parità.\n",
"\n",
"In uno stato GHZ ideale, tutti i correlatori tra vicini più prossimi sono pari a $1$. Sull'hardware, il rumore riduce le stime grezze e introduce una dispersione da qubit a qubit. Dopo il rilevamento degli errori, vengono eliminate le istantanee in cui si è verificato almeno un errore rilevabile, mentre le stime rimanenti vengono avvicinate a un $1$ e con varianza ridotta, il che è coerente con uno stato effettivo più pulito.\n",
"\n",
"Il miglioramento osservato dopo il rilevamento degli errori è coerente con il fatto che i controlli abbiano individuato errori reali.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "75f48e6a-c7e4-46f3-9d39-a7a877427a04",
"metadata": {},
"source": [
"\n",
"\n",
"## Passi successivi\n",
"\n",
"Consulta le seguenti risorse per saperne di più:\n",
"\n",
"\n",
" * [Formulazione generale dell'informazione quantistica](/learning/courses/general-formulation-of-quantum-information)\n",
" * [Tutorial sul rilevamento degli errori a basso overhead tramite codici spazio-temporali](/docs/tutorials/spacetime-codes)\n",
" * [Documentazione sull'API dell'add-on Qiskit Paulice](/docs/api/qiskit-addon-paulice)\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9d876315",
"metadata": {},
"source": [
"\n",
"\n",
"## Riferimenti\n",
"\n",
"* \\[1] Huang, H. Y., Kueng, R., & Preskill, J. (2020). [Prevedere molte proprietà di un sistema quantistico a partire da pochissime misurazioni](https://www.nature.com/articles/s41567-020-0932-7). *Nature Physics*, 16, 1050–1057.\n",
"* \\[2] Zhao, A., Rubin, N. C., e Miyake, A. (2021). [Tomografia parziale fermionica tramite ombre classiche](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.127.110504). *Fis. Rev. Lett.*, 127, 110504.\n",
"* \\[3] Low, G. H. (2022). [Ombre classiche dei fermioni con simmetria del numero di particelle](https://arxiv.org/abs/2208.08964). *arXiv:2208.08964*.\n",
"* \\[4] Bertoni, C., Haferkamp, J., Hinsche, M. et al. (2024). [Ombre superficiali: stima delle aspettative mediante circuiti di Clifford casuali a bassa profondità](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.133.020602). *Fis. Rev. Lett.*, 133, 020602.\n",
"* \\[5] Hu, H.-Y., Gu, A., Majumder, S. et al. (2025). [Dimostrazione di un apprendimento robusto ed efficiente delle proprietà quantistiche mediante ombre poco profonde](https://doi.org/10.1038/s41467-025-57349-w). *Nature Communications*, 16, 2943.\n",
"* \\[6] Martiel, S. & Javadi-Abhari, A. (2025). [Rilevamento degli errori a basso overhead con codici spazio-temporali](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
}