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

Estimación de tiempo de ejecución: menos de 10 segundos en un procesador Heron r2 (NOTA: Se trata únicamente de una estimación). (El tiempo de ejecución puede variar.)

Resultados del aprendizaje

Cómo la tomografía de sombra construye una sombra clásica de un estado cuántico y por qué requiere un número exponencialmente menor de mediciones que la tomografía de estado completa
Cómo implementar la tomografía de sombra de Pauli local utilizando Qiskit para estimar magnetizaciones de un solo qubit, correladores de dos puntos y matrices de densidad reducida de dos qubits
Cómo utilizar qiskit-paulice para la detección de errores con baja sobrecarga mediante códigos espaciotemporales en un experimento de hardware a gran escala

Requisitos previos

En segundo plano

La tomografía completa del estado cuántico de un sistema de qubits « $N$ » requiere un número exponencial de mediciones para estimar todos los parámetros independientes de la matriz de densidad $4^N -1$ . En el caso de un sistema de 50 qubits, por ejemplo, este proceso requiere del orden de $10^{30}$ s, lo que supone más de mil millones de veces el número de estrellas del universo observable. La tomografía de sombra [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 los observables de interés utilizando un número de disparos que varía 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, entre ellas la estimación de energía, los cálculos de la entropía de entrelazamiento y de la función de correlación, la estimación de núcleos en el aprendizaje automático cuántico y otras.

Este tutorial muestra cómo utilizar la tomografía de sombra para estimar observables locales, como la magnetización de un solo qubit $\langle Z_i \rangle$ para el 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), que presenta un alto grado de entrelazamiento. En primer lugar, el tutorial ilustra un flujo de trabajo de tomografía de sombra a pequeña escala en un estado GHZ de 10 qubits mediante la simulación de vectores de estado. A continuación, muestra una demostración de 100 qubits en hardware cuántico real.

Protocolo de tomografía de sombra

El protocolo de tomografía de sombra consta de los siguientes pasos:

Aplica al estado cuántico inicial un « $U$ » unitario aleatorio, extraído de un conjunto unitario.
Mide el resultado en la base computacional, obteniendo una cadena de bits $\ket{b}$ .
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]$ .
Repite el proceso $T$ veces, de modo que el conjunto de instantáneas $T$ $\{\hat{\rho}_t\}$ forme una sombra clásica $S(\rho, T)$ .
Estimar la variable 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 sombra

La elección del conjunto unitario —es decir, la base— en el paso 1 del protocolo influye en la complejidad de la muestra, la varianza de las estimaciones y la exhaustividad de la información de una tomografía de sombra determinada. En la bibliografía se han analizado numerosos métodos de ensamblaje, cada uno de los cuales presenta un equilibrio entre la complejidad y la eficiencia en la estimación de las variables observables. En la siguiente tabla se resumen las propiedades de las variantes más utilizadas, 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	local/pocos cuerpos
Clifford global [1]	$N$ -qubit de Clifford	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	rango bajo, fidelidad
Gaussiano fermiónico [2], [3]	circuitos Matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermiónico
sombras superficiales [4], [5]	circuitos de baja profundidad	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	no local (compatible con el hardware)

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

El término «complejidad de la muestra» se refiere al número de mediciones necesarias para una variante determinada de la tomografía de sombra.
El «grosor del circuito» describe el coste que supone llevar a cabo las mediciones en hardware cuántico.
$N$ es el número de qubits.
$k$ es la ubicación de la variable 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 variables observables que se van a estimar simultáneamente.
$\varepsilon$ es la tolerancia al error aditivo.
$\mu$ es el número de fermiones.

Aunque el conjunto de Pauli local es el que mejor se adapta al hardware, conlleva una penalización en la complejidad de muestreo $3^k$ , lo que lo hace poco adecuado para la estimación de propiedades no locales. En términos más generales, el límite de complejidad de muestreo de la sombra local varía proporcionalmente a $4^k$ , y mejora hasta $3^k$ en el caso especial de los observables de producto tensorial.

Por el contrario, aunque el conjunto de Clifford global puede utilizarse para estimar observables no locales, su profundidad de circuito no es favorable. Las sombras gaussianas fermiónicas resultan útiles en simulaciones de física de muchos cuerpos, química y materiales, pero requieren circuitos «matchgate», lo que podría suponer una carga significativa en dispositivos de uso general basados en puertas lógicas.

«Shallow shadows» presenta una profundidad de circuito situada entre la del circuito local de Pauli y la del circuito global de Clifford, lo que lo convierte en una variante muy prometedora para su aplicación práctica, aunque aún quedan algunas cuestiones por resolver en cuanto a la profundidad óptima del circuito para un observable determinado.

A medida que la investigación en tomografía de sombra sigue avanzando, algunos de los problemas más relevantes que hay que abordar son la elección óptima del conjunto de mediciones para un conjunto objetivo de observables, la mejora de los límites de complejidad de la muestra y la ampliación del protocolo a observables que evolucionan en el tiempo. No obstante, la tomografía de sombra sigue siendo uno de los principales métodos para extraer la máxima información clásica de un estado cuántico utilizando un número mínimo de mediciones, lo que la convierte en la técnica ideal para caracterizar estados cuánticos a una escala en la que tanto la tomografía completa como la simulación clásica resultan imposibles.

El conjunto local de Pauli

Este tutorial utiliza la base aleatoria de Pauli, es decir, puertas de Clifford de un solo qubit, para recopilar sombras y estimar las propiedades del estado GHZ. Hay que tener en cuenta que las medidas de Pauli son una excelente 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 las globales (como la energía total, el momento y la carga).

En este caso, se elige de forma independiente una base de Pauli aleatoria $\{X, Y, Z\}$ por qubit y se aplica la rotación correspondiente a un solo qubit, 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$ se obtiene al promediar las tres bases de Pauli, lo que da lugar a un canal despolarizador con el parámetro $\frac{1}{3}$ . El estimador completo se factoriza entonces como un producto de estimadores de un solo qubit. Tal y como se muestra en la tabla, para observables de tipo « $M$ », en cada una de las localidades $k$ , el número de sombras de Pauli necesarias para estimar todas las observables con un error aditivo $\varepsilon$ con alta probabilidad es $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$ . El tamaño del sistema no afecta al coste de la medición siempre que las observables sigan siendo locales, por lo que se puede estimar simultáneamente un número exponencial de observables con un coste meramente logarítmico.

Mediana de las medias

Una media simple calculada a partir de instantáneas de $T$ reduce la varianza de una estimación de una sola toma de una variable observable según la fórmula $\frac{1}{T}$ , pero es sensible a los valores atípicos. La mediana de las medias puede servir para solucionar este problema. El proceso consiste en lo siguiente:

Dividir las instantáneas de « $T$ » en $K$ grupos de igual tamaño
Cálculo de la media dentro de cada grupo
Devuelve la mediana de las medias de los grupos

En este tutorial se lleva a cabo una tomografía de sombra de Pauli local para estimar las propiedades del estado GHZ; primero, mediante un ejemplo simulado de 10 qubits y, a continuación, mediante un experimento con hardware de 100 qubits.

Requisitos

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

Qiskit SDK v2.0 o posterior, con soporte para visualización
IBM Quantum Compute v0.28 o posterior (pip install qiskit-ibm-runtime)
Qiskit Paulice (pip install qiskit-paulice)

Configuración

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)

Las siguientes funciones auxiliares se utilizan en el ejemplo de hardware a gran escala que aparece más adelante en el tutorial:

bad_cz, bad_readout, bad_coherence — Consultar los datos de calibración del backend para identificar puertas, qubits y flancos que superen los umbrales de error o de coherencia, de modo que puedan excluirse del diseño del circuito.
parallel_ghz — Construye un circuito GHZ de un tamaño determinado mediante una búsqueda en anchura sobre el mapa de acoplamientos, distribuyendo las capas CNOT en paralelo para minimizar la profundidad del circuito.
weighted_coverage, active_wires, z_trace_backward — Calcula qué parte de la superficie de error del circuito GHZ queda protegida por un conjunto dado de comprobaciones de paridad espaciotemporal, utilizando la técnica de propagación hacia atrás descrita en [6].
add_measurements — Añade las mediciones de los qubits GHZ y de los qubits de comprobación a registros clásicos independientes, de modo que puedan seleccionarse posteriormente de forma independiente.
build_shadow_circuits — Genera 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 — Descompone la salida del Sampler en triples (ghz_bits, check_bits, bases) , uno por nota.
adapt_to_flat_bits, estimate_zz_nn_hw — Asignar los bits del registro GHZ al espacio de índices completo del dispositivo y calcular la estimación de la mediana de las medias $\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 de simulador a pequeña escala

En primer lugar, prepara un estado GHZ de 10 qubits mediante simulación con vectores de estado y estima sus propiedades utilizando las sombras de Pauli locales. El estado GHZ es un estado cuántico con el máximo nivel de entrelazamiento, formado por tres o más qubits, que sirve de referencia para el entrelazamiento multipartito. Es un componente básico habitual de los protocolos de comunicación cuántica y criptografía. Además, el circuito GHZ tiene propiedades analíticas totalmente conocidas, lo que lo convierte en ideal para su 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})"
)

Output:

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 sombra rotando el estado GHZ a bases de Pauli aleatorias. Esto refleja lo que hace la primitiva «Sampler» en el 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 con el comando « Qiskit primitives »

Para esta simulación, utiliza 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.")

Output:

Collected 5000 snapshots.

Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado

Define el estimador de sombra y el estimador de la mediana de las medias (MOM) para las bases de medición aleatorias utilizadas para recopilar las instantáneas de sombra. Utiliza estas funciones para estimar las propiedades del estado GHZ —entre ellas, la magnetización de un solo qubit, los correladores de dos qubits y el RDM de dos qubits— a partir de las mediciones recopiladas. A continuación, 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}"
)

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

Visualiza 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 gráfico final muestra la magnitud de las entradas de la matriz de densidad reducida de dos qubits. Para el tamaño del sistema en cuestión se utiliza una colección bastante amplia de instantáneas, 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"])

Output:

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

Evaluar el escalado de las tomas

En este tutorial se utilizan 5.000 instantáneas para estimar las variables observables. El escalado de disparos de la tomografía de sombras se evalúa representando gráficamente el error medio de una de las funciones estimadas en función de diferentes presupuestos de disparos. El error en el correlador de dos qubits decae según la relación $1/\sqrt{T}$ , tal y como cabría esperar de la convergencia estadística estándar de un estimador de la media. La importancia de la tomografía de sombra radica en que el prefactor constante que precede a $1/\sqrt{T}$ no aumenta con el número de qubits en el caso de observables locales, como los que se analizan 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()

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

Ejemplo de hardware a gran escala

La simulación anterior confirma que la tomografía de sombra de Pauli local recupera los observables correctos a partir de un vector de estado sin ruido. La ejecución del mismo protocolo utilizando 100 qubits en hardware real plantea dos retos que no se dan en la simulación: los errores de puerta y de lectura se acumulan a lo largo del circuito GHZ profundo, y para encontrar la disposición óptima de los qubits es necesario realizar una búsqueda en el grafo de acoplamiento del dispositivo, en lugar de elegir una disposición manualmente.

Para hacer frente al primer reto, utiliza la detección de errores con baja sobrecarga mediante códigos espaciotemporales [6]. Se utiliza un pequeño número de qubits auxiliares adyacentes al árbol GHZ como controles de paridad. Cada ancilla se entrelaza con dos de sus vecinas GHZ antes de la medición, por lo que un solo error de inversión de bit en cualquiera de las vecinas invierte el resultado de la ancilla de 0 a 1. Se 1 descartan las instantáneas en las que aparece cualquier lectura «ancilla», lo que deja un subconjunto más limpio a partir del cual estimar el valor de « $\langle Z_i Z_{i+1}\rangle$ ».

El segundo reto se aborda buscando entre las posibles disposiciones de árboles GHZ y seleccionando aquella que maximice la cobertura de la comprobación de paridad, es decir, la fracción de líneas del circuito supervisadas por al menos una comprobación, dentro de un aumento aceptable de la profundidad del circuito.

Pasos 1-4 (ejecutados en un único 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 se dan en la simulación:

El paso 1 caracteriza el ruido del backend (errores de puerta, errores de lectura y tiempos de « T2 ») para que los qubits y los flancos poco fiables queden excluidos del diseño del circuito.
El paso 2 analiza las posibles disposiciones del árbol GHZ para encontrar aquella que maximice la cobertura de la comprobación de paridad dentro de una profundidad de circuito aceptable, sustituyendo así la cadena trivial utilizada en la simulación.
En el paso 3 se construyen los circuitos de sombra (uno por cada base de Pauli aleatoria), se transpilan a puertas nativas del hardware y se envían a la primitiva «Sampler» ( shots=8 cada circuito genera ocho cadenas de bits independientes).
En el paso 4 se realizan postselecciones de instantáneas cuyos qubits de comprobación dan todos un valor de 0 (sin error detectado); a continuación, se estima un « $\langle Z_i Z_{i+1}\rangle$ » con y sin dicha postselección para demostrar la mejora derivada de 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)))

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

Interpretar los resultados

Los dos gráficos muestran un valor de « $\langle Z_i Z_{i+1}\rangle$ » estimado a partir de la tomografía de sombra para cada par de qubits adyacentes en la cadena GHZ de 100 qubits, antes y después de la postselección en las comprobaciones de paridad.

Para que se dé un estado GHZ ideal, todos los correladores de vecinos más cercanos deben ser iguales a $1$ . En el hardware, el ruido reduce las estimaciones brutas e introduce dispersión entre qubits. Tras la detección de errores, se eliminan las instantáneas en las que se ha producido al menos un error detectable, y las estimaciones restantes se acercan a un $1$ e con varianza reducida, lo cual es coherente con un estado efectivo más limpio.

La mejora observada tras la detección de errores concuerda con el hecho de que las comprobaciones detectan errores reales.

Próximos pasos

Echa un vistazo a los siguientes recursos para obtener más información:

Recomendaciones

Referencias

[1] Huang, H. Y., Kueng, R., y Preskill, J. (2020). Predecir muchas propiedades de un sistema cuántico a partir de muy pocas mediciones. Nature Physics, 16, 1050-1057.
[2] Zhao, A., Rubin, N. C. y Miyake, A. (2021). Tomografía parcial de fermiones mediante sombras clásicas. Fís. Rev. Lett., 127, 110504.
[3] Low, G. H. (2022). Sombras clásicas de fermiones con simetría de número de partículas. arXiv:2208.08964.
[4] Bertoni, C., Haferkamp, J., Hinsche, M. y otros (2024). Sombras superficiales: estimación de expectativas mediante circuitos de Clifford aleatorios de baja profundidad. Fís. Rev. Lett., 133, 020602.
[5] Hu, H.-Y., Gu, A., Majumder, S. y otros (2025). Demostración de un aprendizaje robusto y eficiente de propiedades cuánticas mediante «sombras» poco profundas. Nature Communications, 16, 2943.
[6] Martiel, S. y Javadi-Abhari, A. (2025). Detección de errores con baja sobrecarga mediante códigos espaciotemporales. arXiv:2504.15725.

¿Le ha resultado útil esta página?

Informe de un error, de una errata o solicite contenido en GitHub.