Tomografía de sombras del estado de Greenberger-Horne-Zeilinger

Estimación de uso: menos de 10 segundos en un procesador Heron r2 (NOTA: Es solo una estimación. Tu tiempo de ejecución puede variar.)

Resultados de aprendizaje

• Cómo la tomografía de sombras construye una sombra clásica de un estado cuántico y por qué requiere exponencialmente menos mediciones que la tomografía completa del estado

• Cómo implementar la tomografía de sombras de Pauli local con Qiskit para estimar magnetizaciones de un solo qubit, correladores de dos puntos y matrices de densidad reducidas de dos qubits

• Cómo usar qiskit-paulice para la detección de errores de bajo sobrecoste mediante códigos espaciotemporales en un experimento de hardware a gran escala

Requisitos previos

• Mediciones

• Tomografía de estados cuánticos

• Detección de errores mediante códigos espaciotemporales

Antecedentes

La tomografía completa del estado cuántico de un sistema de $N$ qubits requiere un número exponencial de mediciones para estimar los $4^N -1$ parámetros independientes de la matriz de densidad. Para un sistema de 50 qubits, por ejemplo, este proceso requiere del orden de $10^{30}$ shots, que es más de mil millones de veces el número de estrellas del universo observable. La tomografía de sombras [1] es un protocolo de información cuántica que construye una “sombra clásica” del estado cuántico, a partir de la cual se pueden estimar observables objetivo usando un número de shots que escala logarítmicamente con el número de observables y es independiente del tamaño del sistema $N$, lo que la hace relevante en muchas áreas de la computación cuántica, incluidas la estimación de energía, los cálculos de entropía de entrelazamiento y funciones de correlación, la estimación de kernels en aprendizaje automático cuántico, entre otras.

Este tutorial muestra cómo usar la tomografía de sombras para estimar observables locales, como la magnetización de un solo qubit $\langle Z_i \rangle$ del qubit $i$, la función de correlación de dos puntos $\langle Z_i Z_j\rangle$ y la matriz de densidad reducida de dos qubits (2-RDM), del estado de Greenberger-Horne-Zeilinger (GHZ), altamente entrelazado. El tutorial ilustra primero un flujo de trabajo de tomografía de sombras a pequeña escala con un estado GHZ de 10 qubits mediante simulación de vector de estado. Después muestra una demostración de 100 qubits en hardware cuántico real.

Protocolo de tomografía de sombras

El protocolo de tomografía de sombras tiene los siguientes pasos:

1. Aplica una unitaria aleatoria $U$, extraída de un conjunto unitario, al estado cuántico inicial.
2. Mide el resultado en la base computacional, obteniendo una cadena de bits $\ket{b}$.
3. Almacena la instantánea clásica, $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, donde $\mathcal{M}^{-1}$ es la inversa del canal de medición $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Repite el proceso $T$ veces, donde la colección de $T$ instantáneas $\{\hat{\rho}_t\}$ forma una sombra clásica $S(\rho, T)$.
5. Estima el observable $o$ del operador $O$ como $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Variantes de la tomografía de sombras

La elección del conjunto unitario, es decir, de la base, en el paso 1 del protocolo afecta a la complejidad de muestreo, a la varianza de las estimaciones y a la completitud de la información de una tomografía de sombras dada. En la literatura se han considerado muchos grupos de conjuntos, cada uno con compromisos entre complejidad y eficiencia en la estimación de observables. Las propiedades de las variantes más usadas se resumen en la siguiente tabla, aunque esta lista no es exhaustiva.

Variante	Conjunto	Complejidad de muestreo	Profundidad del circuito	Observables adecuados
Pauli local [1]	Paulis de un solo qubit $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	locales/de pocos cuerpos
Clifford global [1]	Clifford de $N$ qubits	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	de rango bajo, fidelidad
gaussiana fermiónica [2], [3]	circuitos matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermiónicos
sombras poco profundas [4], [5]	circuitos de poca profundidad	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	no locales (adecuados para el hardware)

Los términos usados en la tabla se definen de la siguiente manera:

• "complejidad de muestreo" se refiere al número de mediciones necesarias para una variante dada de tomografía de sombras.

• "profundidad del circuito" describe lo costoso que es implementar las mediciones en hardware cuántico.

• $N$ es el número de qubits.

• $k$ es la localidad del observable (por ejemplo, $\langle Z_i \rangle$ tiene $k=1$ y $\langle Z_i Z_j\rangle$ tiene $k=2$).

• $M$ es el número de observables que se estiman simultáneamente.

• $\varepsilon$ es la tolerancia de error aditivo.

• $\mu$ es el número de fermiones.

Aunque el conjunto de Pauli local es el más adecuado para el hardware, paga una penalización de complejidad de muestreo de $3^k$, lo que lo hace poco apropiado para estimar propiedades no locales. De forma más general, la cota de complejidad de muestreo de las sombras locales escala como $4^k$, mejorando a $3^k$ en el caso especial de observables de producto tensorial.

Por el contrario, aunque el conjunto de Clifford global puede usarse para estimar observables no locales, su profundidad de circuito no es favorable. Las sombras gaussianas fermiónicas son útiles en simulaciones de física de muchos cuerpos, química y materiales, pero requieren circuitos matchgate, que pueden conllevar un sobrecoste significativo en dispositivos de propósito general basados en puertas.

Las sombras poco profundas tienen una profundidad de circuito intermedia entre la de Pauli local y la de Clifford global, lo que las convierte en una variante muy prometedora para la implementación práctica, aunque aún quedan algunas preguntas abiertas sobre la profundidad óptima del circuito para un observable dado.

A medida que la frontera de la investigación en tomografía de sombras sigue expandiéndose, algunos de los problemas más relevantes por abordar incluyen elegir de forma óptima el conjunto de mediciones para un conjunto objetivo de observables, mejorar las cotas de complejidad de muestreo y extender el protocolo a observables evolucionados en el tiempo. Aun así, la tomografía de sombras sigue siendo uno de los principales flujos de trabajo para extraer la máxima información clásica de un estado cuántico con un número mínimo de mediciones, lo que la hace ideal para caracterizar estados cuánticos a una escala en la que tanto la tomografía completa como la simulación clásica son imposibles.

El conjunto de Pauli local

Este tutorial usa la base de Pauli aleatoria, es decir, puertas de Clifford de un solo qubit, para recopilar sombras y estimar propiedades del estado GHZ. Ten en cuenta que las mediciones de Pauli son una gran opción para estimar observables locales (como la densidad de espín, las funciones de correlación locales y la densidad de probabilidad de posición), pero no globales (como la energía total, el momento y la carga).

Aquí, se elige de forma independiente una base de Pauli aleatoria $\{X, Y, Z\}$ por cada qubit y se aplica la rotación de un solo qubit correspondiente, seguida de la medición:

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

Para este conjunto, la inversa del canal de medición es:

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

El factor $3$ proviene de promediar sobre las tres bases de Pauli, lo que produce un canal despolarizante con parámetro $\frac{1}{3}$. El estimador completo se factoriza entonces como un producto de estimadores de un solo qubit. Como se muestra en la tabla, para $M$ observables, cada uno de localidad $k$, el número de sombras de Pauli necesarias para estimar todos los observables con error aditivo $\varepsilon$ y alta probabilidad es $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. El tamaño del sistema no afecta al coste de medición mientras los observables sigan siendo locales, por lo que se pueden estimar simultáneamente exponencialmente muchos observables con un coste solo logarítmico.

Mediana de medias

Una media simple sobre $T$ instantáneas reduce la varianza de una estimación de un observable con una sola shot como $\frac{1}{T}$, pero es sensible a los valores atípicos. Se puede usar la mediana de medias para remediarlo. El proceso consiste en lo siguiente:

• Dividir las $T$ instantáneas en $K$ grupos de igual tamaño
• Calcular la media dentro de cada grupo
• Devolver la mediana de las medias de los grupos

Este tutorial realiza tomografía de sombras de Pauli local para estimar propiedades del estado GHZ, primero con un ejemplo simulado de 10 qubits y después con un experimento en hardware de 100 qubits.

Requisitos

Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:

• Qiskit SDK v2.0 o posterior, con soporte de visualización

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

• Qiskit Paulice (pip install qiskit-paulice)

Configuración

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

El ejemplo de hardware a gran escala que aparece más adelante en el tutorial usa las siguientes funciones auxiliares:

• bad_cz, bad_readout, bad_coherence: consultan los datos de calibración del Backend para identificar puertas, qubits y aristas que superan los umbrales de error o coherencia, de modo que puedan excluirse del layout del circuito.

• parallel_ghz: construye un circuito GHZ de un tamaño dado mediante una búsqueda en anchura sobre el mapa de acoplamiento, repartiendo las capas de CNOT en paralelo para minimizar la profundidad del circuito.

• weighted_coverage, active_wires, z_trace_backward: calculan qué parte de la superficie de error del circuito GHZ está protegida por un conjunto dado de comprobaciones de paridad espaciotemporales, usando la técnica de propagación hacia atrás de [6].

• add_measurements: añade las mediciones de los qubits GHZ y de los qubits de comprobación en registros clásicos separados para que puedan posseleccionarse de forma independiente.

• build_shadow_circuits: produce un circuito por cada base de Pauli aleatoria añadiendo las rotaciones de un solo qubit adecuadas antes de la medición.

• parse_hw_results_w_checks: desempaqueta la salida de Sampler en tripletas (ghz_bits, check_bits, bases), una por shot.

• adapt_to_flat_bits, estimate_zz_nn_hw: reasignan los bits del registro GHZ al espacio de índices del dispositivo completo y calculan la estimación de la mediana de medias de $\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

Ejemplo a pequeña escala con simulador

Primero, prepara un estado GHZ de 10 qubits mediante simulación de vector de estado y estima sus propiedades usando sombras de Pauli locales. El estado GHZ es un estado cuántico máximamente entrelazado de tres o más qubits que sirve como referencia para el entrelazamiento multipartito. Es un bloque de construcción habitual en protocolos de comunicación y criptografía cuánticas. Además, el circuito GHZ tiene propiedades analíticas completamente conocidas, lo que lo hace ideal para la verificación en este tutorial.

Paso 1: Asignar entradas clásicas a un problema cuántico

Prepara un estado GHZ ideal de 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)

Paso 2: Optimizar para la ejecución en hardware cuántico

A continuación, define una función que recopile instantáneas de sombras rotando el estado GHZ a bases de Pauli aleatorias. Esto refleja lo que hace la primitiva Sampler en hardware en el ejemplo a gran escala.

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

Paso 3: Ejecutar usando primitivas de Qiskit

Para esta simulación, usa directamente el vector de estado.

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.

Paso 4: Posprocesar y devolver el resultado en el formato clásico deseado

Define el estimador de shadow y el estimador de mediana de medias (MOM) para las bases de medición aleatorias utilizadas para recopilar los snapshots de shadow. Usa estas funciones para estimar propiedades del estado GHZ, incluidas la magnetización de un solo qubit, los correladores de dos qubits y la RDM de dos qubits, a partir de las mediciones recopiladas. Después, compara las estimaciones con los valores exactos.

_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

Visualizar los resultados

El primer gráfico muestra los valores estimados de magnetización de un solo qubit para los distintos qubits. El segundo gráfico muestra la distribución de los valores estimados del correlador de dos puntos para todos los pares de qubits. El último gráfico muestra la magnitud de las entradas de la matriz de densidad reducida de dos qubits. Se usa una colección bastante grande de snapshots para el tamaño del sistema dado, por lo que los resultados se acercan a los exactos.

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

Evaluar el escalado de shots

En este tutorial se usan 5000 instantáneas para estimar los observables. El escalado de shots de la tomografía de sombras se evalúa graficando el error medio de una de las funciones estimadas frente a distintos presupuestos de shots. El error en el correlador de dos qubits decae como $1/\sqrt{T}$, como cabe esperar de la convergencia estadística estándar de un estimador de la media. La importancia de la tomografía de sombras es que el prefactor constante delante de $1/\sqrt{T}$ no crece con el número de qubits para observables locales, como los examinados en este 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

Ejemplo de hardware a gran escala

La simulación anterior confirma que la tomografía de sombras de Pauli local recupera los observables correctos en un vector de estado sin ruido. Ejecutar el mismo protocolo con 100 qubits en hardware real introduce dos desafíos que no aparecen en la simulación: los errores de puertas y de lectura se acumulan a lo largo del circuito GHZ profundo, y encontrar la disposición óptima de qubits requiere buscar en el grafo de acoplamiento del dispositivo en lugar de elegir una disposición a mano.

Para abordar el primer desafío, usa detección de errores de bajo sobrecoste con códigos espaciotemporales [6]. Se utiliza un pequeño número de qubits ancilla adyacentes al árbol GHZ como comprobaciones de paridad. Cada ancilla se entrelaza con dos de sus vecinos GHZ antes de la medición, de modo que un único error de inversión de bit en cualquiera de los vecinos cambia el resultado de la ancilla de 0 a 1. Las instantáneas en las que cualquier ancilla de comprobación lee 1 se descartan, lo que deja un subconjunto más limpio a partir del cual estimar $\langle Z_i Z_{i+1}\rangle$.

El segundo desafío se resuelve buscando entre disposiciones candidatas del árbol GHZ y seleccionando la que maximiza la cobertura de las comprobaciones de paridad, es decir, la fracción de líneas del circuito monitoreadas por al menos una comprobación, dentro de un aumento aceptable de la profundidad del circuito.

Pasos 1-4 (realizados en un solo bloque de código)

Ahora se aplica el mismo flujo de trabajo de cuatro pasos a un estado GHZ de 100 qubits en hardware cuántico real, añadiendo dos características específicas del hardware que no surgen en la simulación:

• El paso 1 caracteriza el ruido del backend (errores de puertas, errores de lectura y tiempos T2) para que los qubits y las aristas poco fiables se excluyan de la disposición del circuito.

• El paso 2 busca entre las posibles disposiciones del árbol GHZ para encontrar una que maximice la cobertura de las comprobaciones de paridad dentro de una profundidad de circuito aceptable, en sustitución de la cadena trivial usada en la simulación.

• El paso 3 construye los circuitos de sombras (uno por cada base de Pauli aleatoria), los transpila a puertas nativas del hardware y los envía a la primitiva Sampler con shots=8 (cada circuito produce ocho cadenas de bits independientes).

• El paso 4 posselecciona los snapshots cuyos qubits de comprobación leen todos 0 (sin error detectado) y luego estima $\langle Z_i Z_{i+1}\rangle$ con y sin esa posselección para demostrar la mejora que aporta la detección de errores.

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

Interpretar los resultados

Los dos gráficos muestran $\langle Z_i Z_{i+1}\rangle$ estimado mediante tomografía de shadow para cada par de qubits vecinos en la cadena GHZ de 100 qubits, antes y después de posseleccionar según las comprobaciones de paridad.

Para un estado GHZ ideal, todos los correladores de vecinos más cercanos son iguales a $1$. En el hardware, el ruido reduce las estimaciones sin procesar e introduce dispersión entre qubits. Tras la detección de errores, se eliminan los snapshots en los que ocurrió al menos un error detectable, y las estimaciones restantes se acercan más a $1$ con menor varianza, lo que es coherente con un estado efectivo más limpio.

La mejora tras la detección de errores es coherente con que las comprobaciones detecten errores reales.

Próximos pasos

Explora los siguientes recursos para saber más:

Recomendaciones

• Formulación general de la información cuántica
• Tutorial de detección de errores de bajo sobrecoste con códigos espaciotemporales
• Documentación de la API del complemento Qiskit Paulice

Referencias

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