Shadow tomography dello stato di Greenberger-Horne-Zeilinger

Stima di utilizzo: meno di 10 secondi su un processore Heron r2 (NOTA: è solo una stima. Il tuo tempo di esecuzione potrebbe variare.)

Obiettivi di apprendimento

• Come la shadow tomography costruisce una shadow classica di uno stato quantistico e perché richiede un numero di misure esponenzialmente inferiore rispetto alla tomografia completa dello stato

• Come implementare la shadow tomography di Pauli locale usando Qiskit per stimare magnetizzazioni a singolo qubit, correlatori a due punti e matrici densità ridotte a due qubit

• Come usare qiskit-paulice per il rilevamento degli errori a basso overhead usando codici spazio-temporali in un esperimento hardware su larga scala

Prerequisiti

• Misure

• Tomografia dello stato quantistico

• Rilevamento degli errori con codici spazio-temporali

Contesto

La tomografia completa dello stato quantistico di un sistema di $N$ qubit richiede un numero esponenziale di misure per stimare tutti i $4^N -1$ parametri indipendenti della matrice densità. Per un sistema di 50 qubit, ad esempio, questo processo richiede nell'ordine di $10^{30}$ shot, ovvero più di un miliardo di volte il numero di stelle nell'universo osservabile. La shadow tomography [1] è un protocollo di informazione quantistica che costruisce una “shadow classica” dello stato quantistico, dalla quale le osservabili di interesse possono essere stimate usando un numero di shot che scala in modo logaritmico con il numero di osservabili ed è indipendente dalla dimensione del sistema $N$. Questo la rende rilevante in molti ambiti del calcolo quantistico, tra cui la stima dell'energia, i calcoli di entropia di entanglement e di funzioni di correlazione, la stima di kernel nel machine learning quantistico e altri.

Questo tutorial mostra come usare la shadow tomography per stimare osservabili locali, come la magnetizzazione a singolo qubit $\langle Z_i \rangle$ per il qubit $i$, la funzione di correlazione a due punti $\langle Z_i Z_j\rangle$ e la matrice densità ridotta a due qubit (2-RDM), dello stato di Greenberger-Horne-Zeilinger (GHZ) altamente entangled. Il tutorial illustra prima un workflow di shadow tomography su piccola scala con uno stato GHZ a 10 qubit usando la simulazione dello statevector. Poi mostra una dimostrazione a 100 qubit su hardware quantistico reale.

Protocollo di shadow tomography

Il protocollo di shadow tomography ha i seguenti passaggi:

1. Applica allo stato quantistico iniziale un'unitaria casuale $U$, estratta da un ensemble unitario.
2. Misura il risultato nella base computazionale, ottenendo una stringa di bit $\ket{b}$.
3. Memorizza lo snapshot classico, $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, dove $\mathcal{M}^{-1}$ è l'inversa del canale di misura $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Ripeti il processo $T$ volte, dove la raccolta di $T$ snapshot $\{\hat{\rho}_t\}$ forma una shadow classica $S(\rho, T)$.
5. Stima l'osservabile $o$ dell'operatore $O$ come $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Varianti della shadow tomography

La scelta dell'ensemble unitario, cioè della base, nel passaggio 1 del protocollo influisce sulla complessità di campionamento, sulla varianza delle stime e sulla completezza dell'informazione di una data shadow tomography. In letteratura sono stati considerati molti gruppi di ensemble, ciascuno con compromessi tra complessità ed efficienza di stima delle osservabili. Le proprietà delle varianti più comuni sono riassunte nella tabella seguente, anche se questo elenco non è esaustivo.

Variante	Ensemble	Complessità di campionamento	Profondità del circuito	Osservabili adatte
Pauli locale [1]	Pauli a singolo qubit $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	locali/a pochi corpi
Clifford globale [1]	Clifford a $N$ qubit	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	a basso rango, fedeltà
gaussiana fermionica [2], [3]	circuiti matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermioniche
shallow shadows [4], [5]	circuiti a bassa profondità	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	non locali (adatte all'hardware)

I termini usati nella tabella sono definiti come segue:

• la "complessità di campionamento" si riferisce al numero di misure necessarie per una data variante di shadow tomography.

• la "profondità del circuito" descrive quanto è costoso implementare le misure su hardware quantistico.

• $N$ è il numero di qubit.

• $k$ è la località dell'osservabile (ad esempio, $\langle Z_i \rangle$ ha $k=1$ e $\langle Z_i Z_j\rangle$ ha $k=2$).

• $M$ è il numero di osservabili da stimare simultaneamente.

• $\varepsilon$ è la tolleranza additiva dell'errore.

• $\mu$ è il numero di fermioni.

Anche se l'ensemble di Pauli locale è il più adatto all'hardware, paga una penalità di complessità di campionamento pari a $3^k$, che lo rende poco adatto alla stima di proprietà non locali. Più in generale, il limite di complessità di campionamento delle shadow locali scala come $4^k$, migliorando a $3^k$ nel caso speciale di osservabili a prodotto tensoriale.

Al contrario, mentre l'ensemble di Clifford globale può essere usato per stimare osservabili non locali, la sua profondità del circuito non è favorevole. Le shadow gaussiane fermioniche sono utili nelle simulazioni di fisica a molti corpi, chimica e materiali, ma richiedono circuiti matchgate, che possono comportare un overhead significativo sui dispositivi generici basati su gate.

Le shallow shadows hanno una profondità del circuito intermedia tra quella di Pauli locale e di Clifford globale, il che le rende una variante molto promettente per l'implementazione pratica, anche se restano alcune questioni aperte sulla profondità ottimale del circuito per una data osservabile.

Mentre la frontiera della ricerca sulla shadow tomography continua ad ampliarsi, alcuni dei problemi più rilevanti da affrontare includono la scelta ottimale dell'ensemble di misura per un insieme di osservabili di interesse, il miglioramento dei limiti di complessità di campionamento e l'estensione del protocollo alle osservabili evolute nel tempo. Ciononostante, la shadow tomography resta uno dei workflow principali per estrarre la massima informazione classica da uno stato quantistico usando un numero minimo di misure, il che la rende ideale per caratterizzare stati quantistici su scale in cui sono impossibili sia la tomografia completa sia la simulazione classica.

L'ensemble di Pauli locale

Questo tutorial usa la base di Pauli casuale, cioè gate di Clifford a singolo qubit, per raccogliere le shadow e stimare le proprietà dello stato GHZ. Nota che le misure di Pauli sono un'ottima scelta per stimare osservabili locali (come la densità di spin, le funzioni di correlazione locali e la densità di probabilità di posizione), ma non quelle globali (come energia totale, quantità di moto e carica).

Qui una base di Pauli casuale $\{X, Y, Z\}$ per ogni qubit viene scelta in modo indipendente e viene applicata la corrispondente rotazione a singolo qubit, seguita dalla misura:

$$\begin{equation} U = U_0 \otimes U_1 \otimes \cdots \otimes U_{n-1}, U_q \in \{I, H, S^{\dagger}H\} \end{equation} \tag{1}$$

Per questo ensemble, l'inversa del canale di misura è:

$$\begin{equation} \hat{\rho}_q=3U_q^{\dagger}\ket{b_q}\bra{b_q}U_q - I \end{equation} \tag{2}$$

Il fattore $3$ deriva dalla media sulle tre basi di Pauli, che produce un canale depolarizzante con parametro $\frac{1}{3}$. Lo stimatore completo si fattorizza quindi come prodotto di stimatori a singolo qubit. Come mostrato nella tabella, per $M$ osservabili, ciascuna di località $k$, il numero di shadow di Pauli necessarie per stimare tutte le osservabili con errore additivo $\varepsilon$ con alta probabilità è $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. La dimensione del sistema non influisce sul costo di misura finché le osservabili restano locali, quindi un numero esponenziale di osservabili può essere stimato simultaneamente con un costo solo logaritmico.

Mediana delle medie

Una semplice media su $T$ snapshot riduce la varianza della stima di un'osservabile da singolo shot come $\frac{1}{T}$, ma è sensibile agli outlier. La mediana delle medie può essere usata per rimediare. Il processo prevede quanto segue:

• Dividere i $T$ snapshot in $K$ gruppi di uguale dimensione
• Calcolare la media all'interno di ciascun gruppo
• Restituire la mediana delle medie dei gruppi

Questo tutorial esegue la shadow tomography di Pauli locale per stimare le proprietà dello stato GHZ, prima con un esempio simulato a 10 qubit e poi dimostrando un esperimento hardware a 100 qubit.

Requisiti

Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:

• Qiskit SDK v2.0 o successivo, con supporto per la visualizzazione

• IBM Quantum Compute v0.28 o successivo (pip install qiskit-ibm-runtime)

• Qiskit Paulice (pip install qiskit-paulice)

Configurazione

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "qiskit_paulice": "qiskit-paulice"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import matplotlib.pyplot as plt
from collections import deque, defaultdict
from qiskit import QuantumCircuit, ClassicalRegister
from qiskit.quantum_info import (
    Statevector,
    partial_trace,
    Operator,
    Pauli,
    SparsePauliOp,
)
from qiskit.transpiler import generate_preset_pass_manager
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit_paulice.layout import get_low_overhead_ancillas

seed = 42  # for reproducibility
rng = np.random.default_rng(seed)

Le seguenti funzioni di supporto sono usate dall'esempio hardware su larga scala più avanti nel tutorial:

• bad_cz, bad_readout, bad_coherence — Interrogano i dati di calibrazione del backend per identificare gate, qubit e connessioni che superano le soglie di errore o di coerenza, così da poterli escludere dal layout del circuito.

• parallel_ghz — Costruisce un circuito GHZ di una data dimensione usando una ricerca in ampiezza sulla coupling map, distribuendo in parallelo i layer di CNOT per minimizzare la profondità del circuito.

• weighted_coverage, active_wires, z_trace_backward — Calcolano quanta parte della superficie di errore del circuito GHZ è protetta da un dato insieme di controlli di parità spazio-temporali, usando la tecnica di propagazione all'indietro da [6].

• add_measurements — Aggiunge le misure dei qubit GHZ e dei qubit di controllo in registri classici separati, in modo che possano essere postselezionate in modo indipendente.

• build_shadow_circuits — Produce un circuito per ogni base di Pauli casuale aggiungendo le opportune rotazioni a singolo qubit prima della misura.

• parse_hw_results_w_checks — Scompone l'output del Sampler in triple (ghz_bits, check_bits, bases), una per shot.

• adapt_to_flat_bits, estimate_zz_nn_hw — Rimappano i bit del registro GHZ nello spazio degli indici dell'intero dispositivo e calcolano la stima della mediana delle medie di $\langle Z_i Z_{i+1}\rangle$.

def bad_cz(target, threshold: float = 0.01) -> list:
    "Return edges whose CZ error exceeds threshold."
    edges = []
    for edge in backend.target.build_coupling_map().get_edges():
        if (edge[1], edge[0]) not in edges:
            edges.append(edge)
    cz_errors = {edge: target["cz"][edge].error for edge in edges}
    worst = sorted(cz_errors.items(), key=lambda x: x[1], reverse=True)
    return [list(edge) for edge, error in worst if error > threshold]

def bad_readout(target, threshold: float = 0.01) -> list:
    "Return nodes whose measurement error exceeds threshold."
    meas_errors = {
        node: target["measure"][(node,)].error
        for node in range(backend.num_qubits)
    }
    worst = sorted(meas_errors.items(), key=lambda x: x[1], reverse=True)
    return [node for node, error in worst if error > threshold]

def bad_coherence(target, threshold: float = 60) -> list:
    "Return nodes whose T2 value is below threshold."
    t2s = {
        node: (target.qubit_properties[node].t2 or 0) * 1e6
        for node in range(backend.num_qubits)
    }
    return [
        node
        for node, val in sorted(t2s.items(), key=lambda x: x[1])
        if val < threshold
    ]

def parallel_ghz(
    root: int, num_qubits: int, backend, bad_edges: list, skip: list
) -> tuple:
    "Build a GHZ circuit of size num_qubits using BFS from root, one neighbor per layer."
    edges = [
        e
        for e in [list(e) for e in backend.coupling_map.get_edges()]
        if e not in bad_edges
        and [e[1], e[0]] not in bad_edges
        and e[0] not in skip
        and e[1] not in skip
    ]
    adj = defaultdict(list)
    for u, v in edges:
        adj[u].append(v)
        adj[v].append(u)

    qc = QuantumCircuit(backend.num_qubits)
    visited = [root]
    queue = deque([root])
    explored = defaultdict(set)
    layers = []
    qc.h(root)

    while queue and len(visited) < num_qubits:
        layer = []
        for node in list(queue):
            queue.popleft()
            unvisited = [
                nb
                for nb in adj[node]
                if nb not in visited and nb not in explored[node]
            ]
            if unvisited:
                nb = unvisited[0]
                visited.append(nb)
                queue.append(node)
                queue.append(nb)
                explored[node].add(nb)
                layer.append((node, nb))
                if len(visited) == num_qubits:
                    break
        if layer:
            layers.append(layer)
        else:
            break

    for layer in layers:
        for q1, q2 in layer:
            qc.cx(q1, q2)
        qc.barrier()

    return qc, visited, layers

def weighted_coverage(
    layers: list, parities: list, w_idle: float = 0.2, w_gate: float = 0.8
) -> float:
    "Compute weighted fraction of wires covered by at least one parity check."
    wires = active_wires(layers)
    covered_by_any = {n_layer: set() for n_layer in range(len(layers))}
    for parity in parities:
        for n_layer, qs in z_trace_backward(layers, parity).items():
            covered_by_any[n_layer] |= qs
    covered_weight = total_weight = 0
    for n_layer in range(len(layers)):
        idle = wires[n_layer]["idle"]
        gate = wires[n_layer]["gate"]
        total_weight += w_idle * len(idle) + w_gate * len(gate)
        covered_weight += w_idle * len(
            covered_by_any[n_layer] & idle
        ) + w_gate * len(covered_by_any[n_layer] & gate)
    return covered_weight / total_weight if total_weight > 0 else 0

def active_wires(layers: list) -> dict:
    """
    Returns per-layer dict with two sets:
    - 'idle': activated wires that are idle in this layer
    - 'gate': activated wires that are control/target of a CNOT at this layer
    """
    first_activation = {}
    for n_layer, layer in enumerate(layers):
        for c, t in layer:
            first_activation.setdefault(c, n_layer)
            first_activation.setdefault(t, n_layer)
    result = {}
    for n_layer in range(len(layers)):
        active = {
            q
            for q, n_layer0 in first_activation.items()
            if n_layer >= n_layer0
        }
        gate = {q for c, t in layers[n_layer] for q in (c, t)}
        result[n_layer] = {"idle": active - gate, "gate": gate}
    return result

def z_trace_backward(layers: list, initial_Zs: list) -> dict:
    """
    Backward propagate Zs with parity cancellation.
    Returns {layer: set of qubits with odd parity Z at that layer}.
    """
    wires = active_wires(layers)
    support = set(initial_Zs)
    trace = {}
    for n_layer in range(len(layers) - 1, -1, -1):
        active = wires[n_layer]["idle"] | wires[n_layer]["gate"]
        trace[n_layer] = support & active
        new_support = set()
        for q in support:
            hit = False
            for c, t in layers[n_layer]:
                if q == t:
                    new_support ^= {t, c}
                    hit = True
                    break
                elif q == c:
                    new_support ^= {c}
                    hit = True
                    break
            if not hit:
                new_support ^= {q}
        support = new_support
    return trace

def add_measurements(
    qc: QuantumCircuit, ghz_qubits: list, checks: list
) -> QuantumCircuit:
    "Add GHZ and check qubit measurements into separate classical registers."
    c1 = ClassicalRegister(len(ghz_qubits), "c1")
    qc.add_register(c1)
    for q, c in zip(ghz_qubits, c1):
        qc.measure(q, c)
    c2 = ClassicalRegister(len(checks), "c2")
    qc.add_register(c2)
    for q, c in zip(checks, c2):
        qc.measure(q, c)
    return qc

def build_shadow_circuits(
    base_circuit: QuantumCircuit,
    n_measure: int,
    ghz_qubits: list,
    checks: list,
) -> tuple:
    "Build n_measure circuits by appending random Pauli basis rotations to each GHZ qubit."
    n = base_circuit.num_qubits
    all_bases = rng.integers(0, 3, size=(n_measure, n))
    circuits = []
    for m in range(n_measure):
        qc = base_circuit.copy()
        qc.barrier()
        for q in ghz_qubits:
            b = all_bases[m, q]
            if b == 1:
                qc.h(q)
            elif b == 2:
                qc.sdg(q)
                qc.h(q)
        add_measurements(qc, ghz_qubits, checks)
        circuits.append(qc)
    return circuits, all_bases

def parse_hw_results_w_checks(
    result, bases: np.ndarray, ghz_qubits: list, checks: list
) -> list:
    "Parse hardware bitstrings into (ghz_bits, check_bits, bases) triples."
    snapshots = []
    n_ghz = len(ghz_qubits)
    n_checks = len(checks)
    for t, pub_result in enumerate(result):
        ghz_strings = (
            pub_result.data.c1.get_bitstrings()
        )  # one entry per shot
        chk_strings = pub_result.data.c2.get_bitstrings()
        for ghz_bitstring, check_bitstring in zip(ghz_strings, chk_strings):
            ghz_bits = [
                int(ghz_bitstring[n_ghz - 1 - i]) for i in range(n_ghz)
            ]
            check_bits = [
                int(check_bitstring[n_checks - 1 - i])
                for i in range(n_checks)
            ]
            snapshots.append((ghz_bits, check_bits, bases[t]))
    return snapshots

def adapt_to_flat_bits(
    ghz_bits: list, ghz_qubits: list, n_total: int
) -> list:
    "Map GHZ-register bits back to a full n_total-length bit vector."
    bits = [0] * n_total
    for i, q in enumerate(ghz_qubits):
        bits[q] = ghz_bits[i]
    return bits

def estimate_zz_nn_hw(
    snapshots: list, ghz_qubits: list, n_total: int, n_groups: int
) -> np.ndarray:
    "Estimate <ZiZi+1> for each neighboring GHZ qubit pair from flat-bit snapshots."
    pairs = list(zip(ghz_qubits[:-1], ghz_qubits[1:]))
    all_est = np.column_stack(
        [
            snapshot_estimate(
                ["Z" if q in pair else "I" for q in range(n_total)], snapshots
            )
            for pair in pairs
        ]
    )
    return mom(all_est, n_groups).real

Esempio su simulatore su piccola scala

Per prima cosa, prepara uno stato GHZ a 10 qubit usando la simulazione dello statevector e stima le sue proprietà usando shadow di Pauli locali. Lo stato GHZ è uno stato quantistico massimamente entangled con tre o più qubit che funge da benchmark per l'entanglement multipartito. È un elemento costruttivo comune per i protocolli di comunicazione e crittografia quantistica. Inoltre, il circuito GHZ ha proprietà analitiche completamente note, il che lo rende ideale per la verifica in questo tutorial.

Passaggio 1: mappa gli input classici su un problema quantistico

Prepara uno stato GHZ ideale a 10 qubit.

def ghz_statevector(num_qubits):
    qc = QuantumCircuit(num_qubits)
    qc.h(0)
    for q in range(num_qubits - 1):
        qc.cx(q, q + 1)
    return Statevector(qc)

num_qubits = 10  # number of qubits
psi = ghz_statevector(num_qubits)
psi.seed(seed)
print(
    f"{num_qubits}-qubit GHZ state ready (norm = {np.linalg.norm(psi):.6f})"
)

10-qubit GHZ state ready (norm = 1.000000)

Passaggio 2: ottimizza per l'esecuzione su hardware quantistico

Poi, definisci una funzione che raccoglie snapshot di shadow ruotando lo stato GHZ in basi di Pauli casuali. Questo rispecchia ciò che fa il primitive Sampler sull'hardware nell'esempio su larga scala.

def basis_circuit(b):
    """0 -> Z basis (I), 1 -> X basis (H), 2 -> Y basis (Sdg;H)."""
    qc = QuantumCircuit(1)
    if b == 1:
        qc.h(0)
    elif b == 2:
        qc.sdg(0)
        qc.h(0)
    # b == 0: Z basis, no rotation
    return qc

_basis = {
    b: basis_circuit(b) for b in (0, 1, 2)
}  # built once, shared with the estimator's _U

def collect_shadow_snapshots(psi, n_snapshots, num_qubits):
    snapshots = []
    for _ in range(n_snapshots):
        bases = rng.integers(0, 3, size=num_qubits)
        rot = QuantumCircuit(num_qubits)
        for q, b in enumerate(bases):
            rot.compose(_basis[b], qubits=[q], inplace=True)
        evolved = psi.evolve(rot)
        evolved.seed(int(rng.integers(0, 2**31)))
        bitstr = next(
            iter(evolved.sample_counts(shots=1))
        )  # Born rule internal
        bits = [int(x) for x in bitstr[::-1]]  # little-endian -> qubit
        snapshots.append((bits, bases))
    return snapshots

Passaggio 3: esegui usando i primitive di Qiskit

Per questa simulazione, usa direttamente lo statevector.

n_snapshots = 5000  # number of shadow snapshots
n_groups = 20  # groups/batches for median-of-means

snapshots = collect_shadow_snapshots(psi, n_snapshots, num_qubits)
print(f"Collected {n_snapshots} snapshots.")

Collected 5000 snapshots.

Passaggio 4: post-elabora e restituisci il risultato nel formato classico desiderato

Definisci lo stimatore shadow e lo stimatore della mediana delle medie (MOM) per le basi di misura casuali usate per raccogliere gli snapshot di shadow. Usa queste funzioni per stimare le proprietà dello stato GHZ, tra cui la magnetizzazione del singolo qubit, i correlatori a due qubit e la RDM a due qubit, a partire dalle misure raccolte. Poi confronta le stime con i valori esatti.

_U = {
    b: Operator(basis_circuit(b)) for b in (0, 1, 2)
}  # rotation Uᵢ per basis
_P = {p: Operator(Pauli(p)) for p in "IXYZ"}  # single-qubit Paulis

def shadow_estimator(P, b, bit):  # P is a label char, e.g. "X"
    UPU = _U[b] @ _P[P] @ _U[b].adjoint()  # Uᵢ P Uᵢ†
    return (3 * UPU.data[bit, bit] - np.trace(_P[P].data)).real

def mom(estimates: np.ndarray, n_groups: int) -> np.ndarray:
    "Median-of-means outlier suppression."
    T = len(estimates)
    spg = T // n_groups  # shots per group
    trimmed = estimates[
        : spg * n_groups
    ]  # drop any leftover shots that don't cleanly fit into a group
    groups = trimmed.reshape(n_groups, spg, *estimates.shape[1:])
    means = groups.mean(axis=1)  # take the mean of each group
    return np.median(means.real, axis=0) + 1j * np.median(
        means.imag, axis=0
    )  # take the median of means for real and imaginary parts separately

def snapshot_estimate(
    pauli_string, snapshots
):  # full-width label, len == num_qubits
    out = np.empty(len(snapshots), dtype=complex)
    for t, (bits, bases) in enumerate(snapshots):
        e = 1.0 + 0j
        for q, P in enumerate(pauli_string):
            e *= shadow_estimator(P, bases[q], bits[q])
        out[t] = e
    return out

Z_est = np.column_stack(
    [
        snapshot_estimate(
            ["Z" if q == i else "I" for q in range(num_qubits)], snapshots
        )
        for i in range(num_qubits)
    ]
)
Z_vals = mom(Z_est, n_groups).real

ZZ_est = np.column_stack(
    [
        snapshot_estimate(
            ["Z" if q in (i, i + 1) else "I" for q in range(num_qubits)],
            snapshots,
        )
        for i in range(num_qubits - 1)
    ]
)
ZZ_nn = mom(ZZ_est, n_groups).real

ZZ_full = np.zeros((n_snapshots, num_qubits, num_qubits), dtype=complex)
for i in range(num_qubits):
    ZZ_full[:, i, i] = 1.0  # diagonal is always 1
    for j in range(i + 1, num_qubits):
        col = snapshot_estimate(
            ["Z" if q in (i, j) else "I" for q in range(num_qubits)],
            snapshots,
        )
        ZZ_full[:, i, j] = col
        ZZ_full[:, j, i] = col  # symmetric
ZZ_mat = mom(ZZ_full, n_groups).real

# 2-qubit reduced density matrix (RDM) for pair (qi, qj)
# the elementwise complex median (real/imag split) for the RDM coefficients generalizes the scalar median in Ref.[1]
qi, qj = 0, num_qubits - 1
terms = []
for a in "IXYZ":
    for b in "IXYZ":
        label = "".join({qi: a, qj: b}.get(q, "I") for q in range(num_qubits))
        c = mom(snapshot_estimate(label, snapshots).reshape(-1, 1), n_groups)[
            0
        ]
        terms.append((b + a, c / 4))
rho_shadow = SparsePauliOp.from_list(terms).to_matrix()

traced_out = [q for q in range(num_qubits) if q not in (qi, qj)]
rho_exact = partial_trace(
    psi, traced_out
).data  # diag(0.5, 0, 0, 0.5) for GHZ

print(f"{'Observable':<22} {'Estimated':>12} {'Ideal':>8} {'Max error':>12}")
print("-" * 57)
print(
    f"{'<Zi> (all qubits)':<22} {'-':>12} {'0':>8} {np.max(np.abs(Z_vals)):>12.4f}"
)
print(
    f"{'<ZiZi+1> mean':<22} {np.mean(ZZ_nn):>12.4f} {'1':>8} {np.max(np.abs(ZZ_nn - 1)):>12.4f}"
)
print(
    f"{'2-RDM Frobenius err':<22} {'-':>12} {'0':>8} {np.linalg.norm(rho_shadow - rho_exact):>12.4f}"
)

Observable                Estimated    Ideal    Max error
---------------------------------------------------------
<Zi> (all qubits)                 -        0       0.0540
<ZiZi+1> mean                1.0040        1       0.1520
2-RDM Frobenius err               -        0       0.1080

Visualizza i risultati

Il primo grafico mostra i valori stimati di magnetizzazione del singolo 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. L'ultimo grafico mostra il modulo degli elementi della matrice densità ridotta a due qubit. Per la dimensione del sistema considerata si usa una raccolta di snapshot piuttosto ampia, quindi i risultati sono vicini a quelli esatti.

fig, axes = plt.subplots(
    1, 3, figsize=(14, 4), gridspec_kw={"width_ratios": [1, 1, 1]}
)

z_exact = psi.expectation_value(
    SparsePauliOp.from_sparse_list([("Z", [0], 1.0)], num_qubits)
).real  # 0.0
zz_exact = psi.expectation_value(
    SparsePauliOp.from_sparse_list([("ZZ", [0, 1], 1.0)], num_qubits)
).real  # 1.0

# <Zi>
axes[0].bar(range(num_qubits), Z_vals, color="#4a90d9", alpha=0.85)
axes[0].axhline(0, color="black", lw=0.8)
axes[0].set_ylim(-0.5, 0.5)
axes[0].set_title(r"$\langle Z_i \rangle$ (ideal = %0.2f)" % z_exact)
axes[0].set_xlabel("Qubit index")
axes[0].grid(axis="y", alpha=0.3)

# <ZiZj>
im = axes[1].imshow(ZZ_mat, vmin=-1, vmax=1, cmap="magma")
plt.colorbar(im, ax=axes[1])
axes[1].set_title(r"$\langle Z_i Z_j \rangle$ (ideal = %0.2f)" % zz_exact)
axes[1].set_xlabel("Qubit j")
axes[1].set_ylabel("Qubit i")

# 2-RDM magnitude
im2 = axes[2].imshow(np.abs(rho_shadow), vmin=0, vmax=0.5, cmap="Blues")
plt.colorbar(im2, ax=axes[2])
axes[2].set_title(r"$|\rho_2|$ for qubits (0, N-1)")
axes[2].set_xticks(range(4))
axes[2].set_yticks(range(4))
axes[2].set_xticklabels(["00", "01", "10", "11"])
axes[2].set_yticklabels(["00", "01", "10", "11"])

[Text(0, 0, '00'), Text(0, 1, '01'), Text(0, 2, '10'), Text(0, 3, '11')]

Valuta la scalabilità rispetto agli shot

In questo tutorial vengono usati 5.000 snapshot per stimare le osservabili. La scalabilità rispetto agli shot della shadow tomography viene valutata tracciando l'errore medio di una delle funzioni stimate in funzione di diversi budget di shot. L'errore nel correlatore a due qubit decresce come $1/\sqrt{T}$, come previsto per la convergenza statistica standard di uno stimatore della media. L'importanza della shadow tomography sta nel fatto che il prefattore costante davanti a $1/\sqrt{T}$ non cresce con il numero di qubit per le osservabili locali, come quelle esaminate in questo tutorial.

shot_counts = [100, 250, 500, 1000, 2000, 5000, 10000]
errors = []

for N_test in shot_counts:
    snaps = collect_shadow_snapshots(psi, N_test, num_qubits)
    zz_vals = []
    for i in range(num_qubits - 1):
        ests = []
        for bits, bases in snaps:
            e = 1.0 + 0j
            for q in range(num_qubits):
                op = "Z" if q in (i, i + 1) else "I"
                e *= shadow_estimator(op, bases[q], bits[q])
            ests.append(e.real)
        zz_vals.append(np.mean(ests))
    errors.append(float(np.mean(np.abs(np.array(zz_vals) - 1.0))))
    print(f"T={N_test:5d} error = {errors[-1]:.4f}")

fig, ax = plt.subplots(figsize=(6, 4))
T_arr = np.array(shot_counts, dtype=float)
ax.loglog(T_arr, errors, "o-", color="#4a90d9", lw=2, label="Shadow estimate")
ax.loglog(
    T_arr,
    errors[0] * np.sqrt(shot_counts[0] / T_arr),
    "k--",
    alpha=0.5,
    label=r"$1/\sqrt{T}$ reference",
)
ax.set_xlabel("Shots T")
ax.set_ylabel(r"Mean $|\langle Z_i Z_{i+1} \rangle - 1|$")
ax.set_title(f"Convergence for {num_qubits}-qubit GHZ")
ax.legend()
plt.tight_layout()
plt.show()

T=  100 error = 0.1856
T=  250 error = 0.1613
T=  500 error = 0.1027
T= 1000 error = 0.0827
T= 2000 error = 0.0513
T= 5000 error = 0.0311
T=10000 error = 0.0162

Esempio su hardware su larga scala

La simulazione precedente conferma che la shadow tomography di Pauli locale recupera le osservabili corrette su uno statevector privo di rumore. Eseguire lo stesso protocollo con 100 qubit su hardware reale introduce due sfide che non compaiono in simulazione: gli errori dei gate e di lettura si accumulano lungo il circuito GHZ profondo, e trovare il layout ottimale dei qubit richiede di cercare nel grafo di accoppiamento del dispositivo invece di scegliere a mano un layout.

Per affrontare la prima sfida, usa il rilevamento degli errori a basso overhead con i codici spazio-temporali [6]. Un piccolo numero di qubit ancilla adiacenti all'albero GHZ viene usato come controlli di parità. Ogni ancilla viene entangled con due dei suoi vicini GHZ prima della misura, quindi un singolo errore di bit-flip su uno dei due vicini inverte l'esito dell'ancilla da 0 a 1. Gli snapshot in cui un qualsiasi ancilla di controllo legge 1 vengono scartati, lasciando un sottoinsieme più pulito da cui stimare $\langle Z_i Z_{i+1}\rangle$.

La seconda sfida viene gestita cercando tra i possibili layout dell'albero GHZ e scegliendo quello che massimizza la copertura dei controlli di parità, cioè la frazione di fili del circuito monitorati da almeno un controllo, entro un aumento accettabile della profondità del circuito.

Passi 1-4 (eseguiti in un unico blocco di codice)

Lo stesso workflow in quattro passi viene ora applicato a uno stato GHZ a 100 qubit su hardware quantistico reale, aggiungendo due funzionalità specifiche dell'hardware che non emergono in simulazione:

• Il passaggio 1 caratterizza il rumore del Backend (errori dei gate, errori di lettura e tempi T2) in modo che qubit e collegamenti inaffidabili siano esclusi dal layout del circuito.

• Il passaggio 2 cerca tra i possibili layout dell'albero GHZ per trovarne uno che massimizzi la copertura dei controlli di parità entro una profondità del circuito accettabile, sostituendo la semplice catena usata in simulazione.

• Il passaggio 3 costruisce i circuiti shadow (uno per ogni base di Pauli casuale), li transpila verso i gate nativi dell'hardware e li invia alla primitiva Sampler con shots=8 (ogni circuito produce otto bitstring indipendenti).

• Il passaggio 4 applica la postselezione agli snapshot i cui qubit di controllo leggono tutti 0 (nessun errore rilevato), quindi stima $\langle Z_i Z_{i+1}\rangle$ con e senza tale postselezione per mostrare il miglioramento dovuto al rilevamento degli errori.

# ── Step 1: Map classical inputs to a quantum problem ─────────────────────────
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)
coupling_map = backend.target.build_coupling_map()
n_total = backend.num_qubits

thresh_cz, thresh_meas, thresh_t2 = 0.025, 0.15, 10
bad_edges = bad_cz(backend.target, threshold=thresh_cz)
bad_nodes_readout = bad_readout(backend.target, threshold=thresh_meas)
dead_qubits = bad_readout(backend.target, threshold=0.4)
bad_nodes = list(
    set(bad_nodes_readout)
    | set(bad_coherence(backend.target, threshold=thresh_t2))
)

# ── Step 2: Optimize for quantum hardware execution ───────────────────────────
# Search over starting qubits and random node-skip sets to find the GHZ layout
# with the best parity-check coverage within an acceptable circuit depth.
ghz_size = 100
max_skips, shuffles, max_depth_increase = 10, 200, 10
w_idle, w_gate = 0.2, 0.8

# Find the root that yields the shallowest GHZ — exclude bad_nodes from the start
best_root = None
base_depth = float("inf")
for root in range(backend.num_qubits):
    if root in bad_nodes:
        continue
    qc, ghz_qubits, _ = parallel_ghz(
        root, ghz_size, backend, bad_edges, bad_nodes
    )
    if len(ghz_qubits) != ghz_size:
        continue
    depth = qc.depth(lambda x: x.operation.num_qubits == 2)
    if depth < base_depth:
        best_root, base_depth = root, depth
if best_root is None:
    raise Exception(
        "No qualifying GHZ root layout found. Relax error thresholds."
    )

qc, ghz_qubits, layers = parallel_ghz(
    best_root, ghz_size, backend, bad_edges, bad_nodes
)

degree_two_nodes = [
    i
    for i in ghz_qubits
    if all(n in ghz_qubits for n in coupling_map.neighbors(i))
    and len(coupling_map.neighbors(i)) >= 2
]

best_covered_fraction = -1
best_qc = best_checks = best_parities = best_layers = best_ghz_qubits = None

for num_skips in range(max_skips):
    for _ in range(shuffles):
        skip_trial = list(
            rng.choice(degree_two_nodes, num_skips, replace=False)
        )
        qc, ghz_qubits, layers = parallel_ghz(
            best_root, ghz_size, backend, bad_edges, skip_trial + bad_nodes
        )
        if len(ghz_qubits) != ghz_size:
            continue
        depth = qc.depth(lambda x: x.operation.num_qubits == 2)
        checks, parities = [], []
        for ancilla, neighbors in get_low_overhead_ancillas(
            coupling_map, ghz_qubits
        ).items():
            if ancilla in dead_qubits:
                continue
            good_neighbors = [
                n
                for n in neighbors
                if [ancilla, n] not in bad_edges
                and [n, ancilla] not in bad_edges
            ]
            if len(good_neighbors) < 2:
                continue
            checks.append(ancilla)
            parities.append((good_neighbors[0], good_neighbors[1]))
            qc.cx(good_neighbors[0], ancilla)
            qc.cx(good_neighbors[1], ancilla)
        covered_fraction = weighted_coverage(
            layers=layers, parities=parities, w_idle=w_idle, w_gate=w_gate
        )
        if (
            covered_fraction > best_covered_fraction
            and depth <= base_depth + max_depth_increase
        ):
            best_covered_fraction = covered_fraction
            (
                best_qc,
                best_ghz_qubits,
                best_checks,
                best_parities,
                best_layers,
            ) = qc, ghz_qubits, checks, parities, layers

qc, ghz_qubits, checks, parities, layers = (
    best_qc,
    best_ghz_qubits,
    best_checks,
    best_parities,
    best_layers,
)
if len(ghz_qubits) != ghz_size:
    raise Exception("No GHZ found. Relax error thresholds.")

print(f"GHZ qubits: {ghz_qubits} ({len(ghz_qubits)})")
print(f"Check qubits: {checks} ({len(checks)})")
print(
    "Covered fraction (no idle): ",
    weighted_coverage(
        layers=layers, parities=parities, w_idle=0.0, w_gate=1.0
    ),
)

# ── Step 3: Execute using Qiskit primitives ───────────────────────────────────
# Build one shadow circuit per random Pauli basis, each repeated for shots_per_basis
# independent bitstrings. Transpile to hardware-native gates, then submit using the
# Sampler primitive. Each circuit yields shots_per_basis snapshots in that basis.
n_hw = 1000
shots_per_basis = 8  # multiple shots per basis circuit
shadow_circs, hw_bases = build_shadow_circuits(qc, n_hw, ghz_qubits, checks)
sampler = Sampler(backend)
sampler.options.environment.job_tags = ["TUT_STGHZS"]
pm = generate_preset_pass_manager(optimization_level=1, backend=backend)
isa_circuit = pm.run(shadow_circs)
job = sampler.run(isa_circuit, shots=shots_per_basis)
result = job.result()

# ── Step 4: Post-process and return result in desired classical format ─────────
# Filter snapshots whose check qubits all read 0 (no detected error), then
# estimate <ZiZi+1> with and without the error-detection postselection.
snapshots_hw = parse_hw_results_w_checks(result, hw_bases, ghz_qubits, checks)
snapshots_kept = [
    (ghz_bits, check_bits, bases)
    for (ghz_bits, check_bits, bases) in snapshots_hw
    if all(b == 0 for b in check_bits)
]
snapshots_kept_adapted = [
    (adapt_to_flat_bits(ghz_bits, ghz_qubits, n_total), bases)
    for (ghz_bits, check_bits, bases) in snapshots_kept
]
snapshots_hw_adapted = [
    (adapt_to_flat_bits(ghz_bits, ghz_qubits, n_total), bases)
    for (ghz_bits, check_bits, bases) in snapshots_hw
]

print("Shadow snapshots collected:", len(snapshots_hw))
print("Shadow snapshots kept post error-detection:", len(snapshots_kept))

n_groups = min(5, len(snapshots_kept))  # groups for median-of-means
ZZ_hw = estimate_zz_nn_hw(snapshots_hw_adapted, ghz_qubits, n_total, n_groups)
ZZ_kept = estimate_zz_nn_hw(
    snapshots_kept_adapted, ghz_qubits, n_total, n_groups
)

# Plot <ZiZi+1> before and after error detection
fig, axes = plt.subplots(1, 2, figsize=(13, 4))

axes[0].plot(ZZ_hw, "o-", color="#e05c5c", markersize=4)
axes[0].axhline(1.0, color="gray", ls="--", alpha=0.7, label="Ideal")
axes[0].set_ylim(-0.6, 3.0)
axes[0].set_title(r"$\langle Z_i Z_{i+1} \rangle$ (without error detection)")
axes[0].set_xlabel("GHZ qubit pair")
axes[0].legend()
axes[0].grid(alpha=0.3)

axes[1].plot(ZZ_kept, "o-", color="#e05c5c", markersize=4)
axes[1].axhline(1.0, color="gray", ls="--", alpha=0.7, label="Ideal")
axes[1].set_ylim(-0.6, 3.0)
axes[1].set_title(r"$\langle Z_i Z_{i+1} \rangle$ (with error detection)")
axes[1].set_xlabel("GHZ qubit pair")
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("Before detection: <ZiZi+1>= ", float(np.average(ZZ_hw)))
print("After detection: <ZiZi+1>= ", float(np.average(ZZ_kept)))

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)
Check qubits: [68, 37, 50, 63, 83, 17, 26, 18, 117, 111, 75] (11)
Covered fraction (no idle):  0.5454545454545454
Shadow snapshots collected: 8000
Shadow snapshots kept post error-detection: 566

Before detection: <ZiZi+1>=  0.573749999999951
After detection: <ZiZi+1>=  0.8817377312951776

Interpreta i risultati

I due grafici mostrano $\langle Z_i Z_{i+1}\rangle$ stimato con la shadow tomography per ogni coppia di qubit vicini nella catena GHZ a 100 qubit, prima e dopo la postselezione sui controlli di parità.

Per uno stato GHZ ideale, tutti i correlatori tra primi vicini sono pari a $1$. Su hardware, il rumore abbassa le stime grezze e introduce dispersione da qubit a qubit. Dopo il rilevamento degli errori, gli snapshot in cui si è verificato almeno un errore rilevabile vengono rimossi e le stime rimanenti si avvicinano a $1$ con varianza ridotta, il che è coerente con uno stato effettivo più pulito.

Il miglioramento dopo il rilevamento degli errori è coerente con il fatto che i controlli rilevino errori reali.

Prossimi passi

Esplora le seguenti risorse per saperne di più:

Raccomandazioni

• General formulation of quantum information
• Low-overhead error detection using spacetime codes tutorial
• Qiskit Paulice add-on API documentation

Riferimenti

• [1] Huang, H. Y., Kueng, R., & Preskill, J. (2020). Predicting many properties of a quantum system from very few measurements. Nature Physics, 16, 1050–1057.

• [2] Zhao, A., Rubin, N. C., & Miyake, A. (2021). Fermionic partial tomography via classical shadows. Phys. Rev. Lett., 127, 110504.

• [3] Low, G. H. (2022). Classical shadows of fermions with particle number symmetry. arXiv:2208.08964.

• [4] Bertoni, C., Haferkamp, J., Hinsche, M. et al. (2024). Shallow Shadows: Expectation Estimation Using Low-Depth Random Clifford Circuits. Phys. Rev. Lett., 133, 020602.

• [5] Hu, H.-Y., Gu, A., Majumder, S. et al. (2025). Demonstration of robust and efficient quantum property learning with shallow shadows. Nature Communications, 16, 2943.

• [6] Martiel, S. & Javadi-Abhari, A. (2025). Low-overhead error detection with spacetime codes. arXiv:2504.15725. © IBM Corp. 2026