Tomographie par ombres de l'état de Greenberger-Horne-Zeilinger

Estimation d'utilisation : moins de 10 secondes sur un processeur Heron r2 (REMARQUE : il ne s'agit que d'une estimation. Ton temps d'exécution peut varier.)

Objectifs d'apprentissage

• Comment la tomographie par ombres construit une ombre classique d'un état quantique et pourquoi elle nécessite exponentiellement moins de mesures que la tomographie d'état complète

• Comment implémenter la tomographie par ombres de Pauli locale avec Qiskit pour estimer des aimantations à un qubit, des corrélateurs à deux points et des matrices densité réduites à deux qubits

• Comment utiliser qiskit-paulice pour une détection d'erreurs à faible surcoût à l'aide de codes espace-temps dans une expérience matérielle à grande échelle

Prérequis

• Mesures

• Tomographie d'état quantique

• Détection d'erreurs à l'aide de codes espace-temps

Contexte

La tomographie d'état quantique complète d'un système de $N$ qubits nécessite un nombre exponentiel de mesures pour estimer les $4^N -1$ paramètres indépendants de la matrice densité. Pour un système de 50 qubits, par exemple, ce processus requiert de l'ordre de $10^{30}$ shots, soit plus d'un milliard de fois le nombre d'étoiles de l'univers observable. La tomographie par ombres [1] est un protocole d'information quantique qui construit une « ombre classique » de l'état quantique, à partir de laquelle des observables cibles peuvent être estimées avec un nombre de shots qui croît de façon logarithmique avec le nombre d'observables et est indépendant de la taille du système $N$, ce qui la rend pertinente dans de nombreux domaines du calcul quantique, notamment l'estimation d'énergie, le calcul d'entropie d'intrication et de fonctions de corrélation, l'estimation de noyaux en apprentissage automatique quantique, et bien d'autres.

Ce tutoriel montre comment utiliser la tomographie par ombres pour estimer des observables locales, comme l'aimantation à un qubit $\langle Z_i \rangle$ pour le qubit $i$, la fonction de corrélation à deux points $\langle Z_i Z_j\rangle$ et la matrice densité réduite à deux qubits (2-RDM), de l'état de Greenberger-Horne-Zeilinger (GHZ), fortement intriqué. Le tutoriel illustre d'abord un workflow de tomographie par ombres à petite échelle sur un état GHZ à 10 qubits à l'aide d'une simulation de vecteur d'état. Il présente ensuite une démonstration à 100 qubits sur du matériel quantique réel.

Protocole de tomographie par ombres

Le protocole de tomographie par ombres comporte les étapes suivantes :

1. Appliquer un unitaire aléatoire $U$, tiré d'un ensemble unitaire, à l'état quantique initial.
2. Mesurer le résultat dans la base computationnelle, ce qui donne une chaîne de bits $\ket{b}$.
3. Stocker l'instantané classique, $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, où $\mathcal{M}^{-1}$ est l'inverse du canal de mesure $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Répéter le processus $T$ fois, la collection des $T$ instantanés $\{\hat{\rho}_t\}$ formant une ombre classique $S(\rho, T)$.
5. Estimer l'observable $o$ de l'opérateur $O$ par $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Variantes de la tomographie par ombres

Le choix de l'ensemble unitaire, c'est-à-dire de la base, à l'étape 1 du protocole influence la complexité d'échantillonnage, la variance des estimations et l'exhaustivité de l'information d'une tomographie par ombres donnée. De nombreux groupes d'ensembles ont été étudiés dans la littérature, chacun avec des compromis entre complexité et efficacité d'estimation des observables. Les propriétés des variantes les plus couramment utilisées sont résumées dans le tableau suivant, bien que cette liste ne soit pas exhaustive.

Variante	Ensemble	Complexité d'échantillonnage	Profondeur du Circuit	Observables adaptées
Pauli local [1]	Paulis à un qubit $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	locales/à peu de corps
Clifford global [1]	Clifford à $N$ qubits	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	de rang faible, fidélité
gaussien fermionique [2], [3]	Circuits matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermioniques
ombres peu profondes [4], [5]	Circuits de faible profondeur	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	non locales (adaptées au matériel)

Les termes utilisés dans le tableau sont définis comme suit :

• « complexité d'échantillonnage » désigne le nombre de mesures nécessaires pour une variante de tomographie par ombres donnée.

• « profondeur du Circuit » décrit le coût d'implémentation des mesures sur du matériel quantique.

• $N$ est le nombre de qubits.

• $k$ est la localité de l'observable (par exemple, $\langle Z_i \rangle$ a $k=1$ et $\langle Z_i Z_j\rangle$ a $k=2$).

• $M$ est le nombre d'observables à estimer simultanément.

• $\varepsilon$ est la tolérance d'erreur additive.

• $\mu$ est le nombre de fermions.

Bien que l'ensemble de Pauli local soit le plus adapté au matériel, il subit une pénalité de complexité d'échantillonnage en $3^k$, ce qui le rend mal adapté à l'estimation de propriétés non locales. Plus généralement, la borne de complexité d'échantillonnage des ombres locales croît en $4^k$, et s'améliore à $3^k$ dans le cas particulier des observables en produit tensoriel.

À l'inverse, si l'ensemble de Clifford global peut servir à estimer des observables non locales, sa profondeur de Circuit n'est pas favorable. Les ombres gaussiennes fermioniques sont utiles pour les simulations de physique à N corps, de chimie et de matériaux, mais nécessitent des Circuits matchgate, qui peuvent entraîner un surcoût important sur des dispositifs généralistes à base de Gates.

Les ombres peu profondes ont une profondeur de Circuit intermédiaire entre le Pauli local et le Clifford global, ce qui en fait une variante très prometteuse pour une mise en œuvre pratique, même s'il reste des questions ouvertes sur la profondeur de Circuit optimale pour une observable donnée.

Alors que la frontière de la recherche sur la tomographie par ombres continue de s'élargir, parmi les problèmes les plus pertinents à traiter figurent le choix optimal de l'ensemble de mesure pour un ensemble cible d'observables, l'amélioration des bornes de complexité d'échantillonnage et l'extension du protocole aux observables évoluant dans le temps. Néanmoins, la tomographie par ombres reste l'un des principaux workflows pour extraire un maximum d'information classique d'un état quantique avec un minimum de mesures, ce qui la rend idéale pour caractériser des états quantiques à une échelle où la tomographie complète comme la simulation classique sont impossibles.

L'ensemble de Pauli local

Ce tutoriel utilise la base de Pauli aléatoire, c'est-à-dire des Gates de Clifford à un qubit, pour collecter des ombres et estimer des propriétés de l'état GHZ. Note que les mesures de Pauli sont un excellent choix pour estimer des observables locales (comme la densité de spin, les fonctions de corrélation locales et la densité de probabilité de position), mais pas des observables globales (comme l'énergie totale, l'impulsion et la charge).

Ici, une base de Pauli aléatoire $\{X, Y, Z\}$ est choisie indépendamment pour chaque qubit et la rotation à un qubit correspondante est appliquée, suivie de la mesure :

$$\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}$$

Pour cet ensemble, l'inverse du canal de mesure est :

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

Le facteur $3$ provient de la moyenne sur les trois bases de Pauli, qui produit un canal dépolarisant de paramètre $\frac{1}{3}$. L'estimateur complet se factorise alors en un produit d'estimateurs à un qubit. Comme le montre le tableau, pour $M$ observables, chacune de localité $k$, le nombre d'ombres de Pauli nécessaires pour estimer toutes les observables avec une erreur additive $\varepsilon$ et une forte probabilité est $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. La taille du système n'affecte pas le coût de mesure tant que les observables restent locales, de sorte qu'un nombre exponentiel d'observables peut être estimé simultanément avec seulement un coût logarithmique.

Médiane des moyennes

Une moyenne simple sur $T$ instantanés réduit la variance d'une estimation à un seul shot d'une observable en $\frac{1}{T}$, mais elle est sensible aux valeurs aberrantes. La médiane des moyennes permet d'y remédier. Le processus consiste à :

• Répartir les $T$ instantanés en $K$ groupes de taille égale
• Calculer la moyenne au sein de chaque groupe
• Renvoyer la médiane des moyennes des groupes

Ce tutoriel réalise une tomographie par ombres de Pauli locale pour estimer des propriétés de l'état GHZ, d'abord avec un exemple simulé à 10 qubits, puis en présentant une expérience matérielle à 100 qubits.

Prérequis techniques

Avant de commencer ce tutoriel, assure-toi d'avoir installé les éléments suivants :

• Qiskit SDK v2.0 ou ultérieur, avec la prise en charge de la visualisation

• IBM Quantum Compute v0.28 ou ultérieur (pip install qiskit-ibm-runtime)

• Qiskit Paulice (pip install qiskit-paulice)

Configuration

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

Les fonctions utilitaires suivantes sont utilisées par l'exemple matériel à grande échelle plus loin dans le tutoriel :

• bad_cz, bad_readout, bad_coherence — Interrogent les données de calibration du Backend pour identifier les Gates, qubits et arêtes qui dépassent des seuils d'erreur ou de cohérence, afin de pouvoir les exclure du layout du Circuit.

• parallel_ghz — Construit un Circuit GHZ d'une taille donnée à l'aide d'un parcours en largeur de la carte de couplage, en répartissant les couches de CNOT en parallèle pour minimiser la profondeur du Circuit.

• weighted_coverage, active_wires, z_trace_backward — Calculent quelle part de la surface d'erreur du Circuit GHZ est protégée par un ensemble donné de contrôles de parité espace-temps, à l'aide de la technique de rétropropagation de [6].

• add_measurements — Ajoute les mesures des qubits GHZ et des qubits de contrôle dans des registres classiques séparés afin de pouvoir les post-sélectionner indépendamment.

• build_shadow_circuits — Produit un Circuit par base de Pauli aléatoire en ajoutant les rotations à un qubit appropriées avant la mesure.

• parse_hw_results_w_checks — Décompose la sortie du Sampler en triplets (ghz_bits, check_bits, bases), un par shot.

• adapt_to_flat_bits, estimate_zz_nn_hw — Remettent les bits du registre GHZ dans l'espace d'indices complet du dispositif et calculent l'estimation de $\langle Z_i Z_{i+1}\rangle$ par médiane des moyennes.

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

Exemple à petite échelle sur simulateur

Prépare d'abord un état GHZ à 10 qubits à l'aide d'une simulation de vecteur d'état et estime ses propriétés avec des ombres de Pauli locales. L'état GHZ est un état quantique maximalement intriqué d'au moins trois qubits, qui sert de référence pour l'intrication multipartite. C'est un élément de base courant des protocoles de communication et de cryptographie quantiques. De plus, le Circuit GHZ a des propriétés analytiques entièrement connues, ce qui le rend idéal pour la vérification dans ce tutoriel.

Étape 1 : Traduire des entrées classiques en un problème quantique

Prépare un état GHZ idéal à 10 qubits.

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)

Étape 2 : Optimiser pour l'exécution sur matériel quantique

Définis ensuite une fonction qui collecte des instantanés d'ombre en faisant pivoter l'état GHZ vers des bases de Pauli aléatoires. Cela reproduit ce que fait la primitive Sampler sur le matériel dans l'exemple à grande échelle.

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

Étape 3 : Exécuter avec les primitives Qiskit

Pour cette simulation, utilise directement le vecteur d'état.

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.

Étape 4 : Post-traiter et renvoyer le résultat dans le format classique souhaité

Définis l'estimateur de shadow et l'estimateur médiane des moyennes (MOM) pour les bases de mesure aléatoires utilisées afin de collecter les snapshots de shadow. Utilise ces fonctions pour estimer des propriétés de l'état GHZ, notamment la magnétisation à un seul Qubit, les corrélateurs à deux Qubits et la RDM à deux Qubits, à partir des mesures collectées. Compare ensuite les estimations avec les valeurs exactes.

_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

Visualiser les résultats

Le premier graphique montre les valeurs estimées de la magnétisation à un seul Qubit pour les différents Qubits. Le deuxième graphique montre la distribution des valeurs estimées des corrélateurs à deux points pour toutes les paires de Qubits. Le dernier graphique montre la magnitude des éléments de la matrice densité réduite à deux Qubits. Une collection de snapshots assez grande est utilisée pour la taille de système donnée, de sorte que les résultats sont proches des valeurs exactes.

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

Évaluer le scaling du nombre de shots

Dans ce tutoriel, 5 000 instantanés sont utilisés pour estimer les observables. L'évolution de l'erreur avec le nombre de shots en tomographie par ombres est évaluée en traçant l'erreur moyenne de l'une des fonctions estimées en fonction de différents budgets de shots. L'erreur du corrélateur à deux Qubits décroît comme $1/\sqrt{T}$, comme on peut s'y attendre pour la convergence statistique standard d'un estimateur de moyenne. L'intérêt de la tomographie par ombres est que le préfacteur constant devant $1/\sqrt{T}$ ne croît pas avec le nombre de Qubits pour les observables locales, comme celles examinées dans ce tutoriel.

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

Exemple matériel à grande échelle

La simulation précédente confirme que la tomographie par ombres de Pauli locale retrouve les bonnes observables sur un vecteur d'état sans bruit. Exécuter le même protocole avec 100 Qubits sur du matériel réel introduit deux défis qui n'apparaissent pas en simulation : les erreurs de Gates et de lecture s'accumulent tout au long du Circuit GHZ profond, et trouver la disposition optimale des Qubits nécessite de parcourir le graphe de couplage de l'appareil plutôt que de choisir une disposition à la main.

Pour relever le premier défi, utilise la détection d'erreurs à faible surcoût avec des codes espace-temps [6]. Un petit nombre de Qubits ancillaires adjacents à l'arbre GHZ sont utilisés comme contrôles de parité. Chaque ancilla est intriquée avec deux de ses voisins GHZ avant la mesure, de sorte qu'une seule erreur de bit-flip sur l'un ou l'autre voisin fait passer le résultat de l'ancilla de 0 à 1. Les instantanés dans lesquels une ancilla de contrôle donne 1 sont écartés, ce qui laisse un sous-ensemble plus propre à partir duquel estimer $\langle Z_i Z_{i+1}\rangle$.

Le second défi est traité en parcourant des dispositions candidates de l'arbre GHZ et en sélectionnant celle qui maximise la couverture des contrôles de parité, c'est-à-dire la fraction des fils du Circuit surveillés par au moins un contrôle, pour une augmentation acceptable de la profondeur du Circuit.

Étapes 1 à 4 (réalisées dans un seul bloc de code)

Le même workflow en quatre étapes est maintenant appliqué à un état GHZ de 100 Qubits sur du matériel quantique réel, en ajoutant deux fonctionnalités spécifiques au matériel qui n'apparaissent pas en simulation :

• L'étape 1 caractérise le bruit du Backend (erreurs de Gates, erreurs de lecture et temps T2) afin d'exclure de la disposition du Circuit les Qubits et les arêtes peu fiables.

• L'étape 2 parcourt les dispositions possibles de l'arbre GHZ pour en trouver une qui maximise la couverture des contrôles de parité pour une profondeur de Circuit acceptable, en remplacement de la chaîne triviale utilisée en simulation.

• L'étape 3 construit les Circuits de shadow (un par base de Pauli aléatoire), les transpile en Gates natives du matériel et les soumet à la primitive Sampler avec shots=8 (chaque Circuit produit huit chaînes de bits indépendantes).

• L'étape 4 postsélectionne les snapshots dont les Qubits de contrôle donnent tous 0 (aucune erreur détectée), puis estime $\langle Z_i Z_{i+1}\rangle$ avec et sans cette postsélection pour démontrer l'amélioration apportée par la détection d'erreurs.

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

Interpréter les résultats

Les deux graphiques montrent $\langle Z_i Z_{i+1}\rangle$ estimé par tomographie de shadow pour chaque paire de Qubits voisins de la chaîne GHZ de 100 Qubits, avant et après la postsélection sur les contrôles de parité.

Pour un état GHZ idéal, tous les corrélateurs entre plus proches voisins valent $1$. Sur le matériel, le bruit abaisse les estimations brutes et introduit une dispersion d'un Qubit à l'autre. Après la détection d'erreurs, les snapshots dans lesquels au moins une erreur détectable s'est produite sont supprimés, et les estimations restantes se rapprochent de $1$ avec une variance réduite, ce qui est cohérent avec un état effectif plus propre.

L'amélioration après la détection d'erreurs est cohérente avec le fait que les contrôles détectent de véritables erreurs.

Étapes suivantes

Explore les ressources suivantes pour en savoir plus :

Recommandations

• Formulation générale de l'information quantique
• Tutoriel sur la détection d'erreurs à faible surcoût avec des codes spacetime
• Documentation de l'API de l'add-on Qiskit Paulice

Références

• [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