Tomografia d'ombra dello stato di Greenberger-Horne-Zeilinger

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.)

Risultati di apprendimento

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
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
Come utilizzare qiskit-paulice per il rilevamento degli errori a basso overhead tramite codici spazio-temporali in un esperimento hardware su larga scala

Prerequisiti

Sfondo

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] è 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.

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.

Protocollo di tomografia delle ombre

Il protocollo di tomografia dell'ombra prevede le seguenti fasi:

Applicare allo stato quantistico iniziale un’ $U$ e unitaria casuale, estratta da un insieme unitario.
Misurare il risultato nella base computazionale, ottenendo una stringa di bit $\ket{b}$ .
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]$ .
Ripetere il processo $T$ volte, dove l’insieme delle istantanee $T$ $\{\hat{\rho}_t\}$ forma una classica “ombra” $S(\rho, T)$ .
Stima l' $o$ e osservabile dell'operatore $O$ come $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$ .

Varianti della tomografia delle ombre

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.

Variante	Insieme	Complessità di campionamento	Profondità del circuito	Grandezze osservabili idonee
Pauli locale [1]	Paulis a singolo qubit $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	locale/a pochi corpi
Clifford globale [1]	$N$ -qubit di Clifford	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	basso rango, fedeltà
Gaussiano fermionico [2], [3]	circuiti matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermionico
ombre superficiali [4], [5]	circuiti a bassa profondità	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	non locale (compatibile con l'hardware)

I termini utilizzati nella tabella sono definiti come segue:

Per "complessità del campione" si intende il numero di misurazioni necessarie per una determinata variante della tomografia d'ombra.
Il termine "profondità del circuito" descrive quanto sia oneroso implementare le misurazioni sull'hardware quantistico.
$N$ è il numero di qubit.
$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$ ).
$M$ è il numero di grandezze osservabili da stimare contemporaneamente.
$\varepsilon$ è la tolleranza sull'errore additivo.
$\mu$ è il numero di fermioni.

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.

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.

"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.

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.

L'ensemble locale Pauli

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).

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:

\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 insieme, l'inverso 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 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.

Mediana delle medie

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:

Suddivisione delle istantanee di “ $T$ ” in $K$ gruppi di uguale dimensione
Calcolo della media all’interno di ciascun gruppo
Calcolo della mediana delle medie dei gruppi

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.

Requisiti

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

Qiskit SDK v2.0 o versioni successive, con supporto alla visualizzazione
IBM Quantum Compute v0.28 o versioni successive (pip install qiskit-ibm-runtime)
Qiskit Paulice (pip install qiskit-paulice)

Configura

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 vengono utilizzate nell'esempio relativo all'hardware su larga scala presentato più avanti nel tutorial:

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.
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.
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].
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.
build_shadow_circuits — Genera un circuito per ogni base di Pauli casuale aggiungendo le rotazioni appropriate di un singolo qubit prima della misurazione.
parse_hw_results_w_checks — Suddivide l'output del Sampler in triplette (ghz_bits, check_bits, bases) , una per ogni shot.
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$ ).

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 di simulatore su piccola scala

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.

Fase 1: Mappare gli input classici su un problema quantistico

Preparare 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})"
)

Output:

10-qubit GHZ state ready (norm = 1.000000)

Fase 2: Ottimizzazione per l'esecuzione su hardware quantistico

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.

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

Fase 3: Eseguire il comando utilizzando Qiskit primitives

Per questa simulazione, utilizzare direttamente il vettore di stato.

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.")

Output:

Collected 5000 snapshots.

Fase 4: Elaborazione finale e restituzione del risultato nel formato classico desiderato

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.

_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}"
)

Output:

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 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.

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"])

Output:

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

Valutare il ridimensionamento delle immagini

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.

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()

Output:

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 di hardware su larga scala

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.

Per affrontare la prima sfida, utilizzare un sistema di rilevamento degli errori a basso overhead con codici spazio-temporali [6]. 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.

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.

Passaggi da 1 a 4 (eseguiti in un unico blocco di codice)

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:

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.
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.
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).
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.

# ── 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)))

Output:

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

Interpretare i risultati

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à.

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.

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

Passi successivi

Consulta le seguenti risorse per saperne di più:

Suggerimenti

Riferimenti

[1] Huang, H. Y., Kueng, R., & Preskill, J. (2020). Prevedere molte proprietà di un sistema quantistico a partire da pochissime misurazioni. Nature Physics, 16, 1050–1057.
[2] Zhao, A., Rubin, N. C., e Miyake, A. (2021). Tomografia parziale fermionica tramite ombre classiche. Fis. Rev. Lett., 127, 110504.
[3] Low, G. H. (2022). Ombre classiche dei fermioni con simmetria del numero di particelle. arXiv:2208.08964.
[4] Bertoni, C., Haferkamp, J., Hinsche, M. et al. (2024). Ombre superficiali: stima delle aspettative mediante circuiti di Clifford casuali a bassa profondità. Fis. Rev. Lett., 133, 020602.
[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. Nature Communications, 16, 2943.
[6] Martiel, S. & Javadi-Abhari, A. (2025). Rilevamento degli errori a basso overhead con codici spazio-temporali. arXiv:2504.15725.

Questa pagina è stata utile?

Segnala un bug, un errore di battitura o richiedi contenuti su GitHub.