Shadow tomography do estado de Greenberger-Horne-Zeilinger

Estimativa de uso: menos de 10 segundos em um processador Heron r2 (OBSERVAÇÃO: esta é apenas uma estimativa. Seu tempo de execução pode variar.)

Resultados de aprendizado

• Como a shadow tomography constrói uma sombra clássica de um estado quântico e por que ela exige exponencialmente menos medições do que a tomografia completa de estado

• Como implementar a shadow tomography de Pauli local usando o Qiskit para estimar magnetizações de um único qubit, correlatores de dois pontos e matrizes densidade reduzidas de dois qubits

• Como usar o qiskit-paulice para detecção de erros de baixo overhead usando códigos espaço-temporais em um experimento de hardware em larga escala

Pré-requisitos

• Medições

• Tomografia de estado quântico

• Detecção de erros usando códigos espaço-temporais

Contexto

A tomografia completa de estado quântico de um sistema de $N$ qubits exige um número exponencial de medições para estimar todos os $4^N -1$ parâmetros independentes da matriz densidade. Para um sistema de 50 qubits, por exemplo, esse processo exige da ordem de $10^{30}$ shots, o que é mais de um bilhão de vezes o número de estrelas no universo observável. A shadow tomography [1] é um protocolo de informação quântica que constrói uma “sombra clássica” do estado quântico, a partir da qual observáveis-alvo podem ser estimados usando um número de shots que cresce logaritmicamente com o número de observáveis e é independente do tamanho do sistema $N$, o que a torna relevante em muitas áreas da computação quântica, incluindo estimativa de energia, cálculos de entropia de emaranhamento e funções de correlação, estimativa de kernels em aprendizado de máquina quântico, entre outras.

Este tutorial mostra como usar a shadow tomography para estimar observáveis locais, como a magnetização de um único qubit $\langle Z_i \rangle$ para o qubit $i$, a função de correlação de dois pontos $\langle Z_i Z_j\rangle$ e a matriz densidade reduzida de dois qubits (2-RDM), do estado altamente emaranhado de Greenberger-Horne-Zeilinger (GHZ). O tutorial primeiro ilustra um fluxo de trabalho de shadow tomography em pequena escala com um estado GHZ de 10 qubits usando simulação de vetor de estado. Em seguida, mostra uma demonstração com 100 qubits em hardware quântico real.

Protocolo de shadow tomography

O protocolo de shadow tomography tem as seguintes etapas:

1. Aplique um unitário aleatório $U$, sorteado de um ensemble unitário, ao estado quântico inicial.
2. Meça o resultado na base computacional, obtendo uma bitstring $\ket{b}$.
3. Armazene o snapshot clássico, $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, em que $\mathcal{M}^{-1}$ é o inverso do canal de medição $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Repita o processo $T$ vezes, em que a coleção de $T$ snapshots $\{\hat{\rho}_t\}$ forma uma sombra clássica $S(\rho, T)$.
5. Estime o observável $o$ do operador $O$ como $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Variantes da shadow tomography

A escolha do ensemble unitário, isto é, da base, na etapa 1 do protocolo afeta a complexidade amostral, a variância das estimativas e a completude de informação de uma dada shadow tomography. Muitos grupos de ensembles já foram considerados na literatura, cada um com compromissos entre complexidade e eficiência de estimativa de observáveis. As propriedades das variantes mais usadas estão resumidas na tabela a seguir, embora esta lista não seja exaustiva.

Variante	Ensemble	Complexidade amostral	Profundidade do circuito	Observáveis adequados
Pauli local [1]	Paulis de um qubit $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	locais/de poucos corpos
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 posto baixo, fidelidade
gaussiana fermiônica [2], [3]	circuitos matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermiônicos
shallow shadows [4], [5]	circuitos de baixa profundidade	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	não locais (adequados ao hardware)

Os termos usados na tabela são definidos a seguir:

• "complexidade amostral" refere-se ao número de medições necessárias para uma dada variante de shadow tomography.

• "profundidade do circuito" descreve quão caro é implementar as medições em hardware quântico.

• $N$ é o número de qubits.

• $k$ é a localidade do observável (por exemplo, $\langle Z_i \rangle$ tem $k=1$ e $\langle Z_i Z_j\rangle$ tem $k=2$).

• $M$ é o número de observáveis a serem estimados simultaneamente.

• $\varepsilon$ é a tolerância de erro aditivo.

• $\mu$ é o número de férmions.

Embora o ensemble de Pauli local seja o mais adequado ao hardware, ele paga uma penalidade de complexidade amostral de $3^k$, o que o torna pouco indicado para a estimativa de propriedades não locais. De forma mais geral, o limite de complexidade amostral das sombras locais cresce como $4^k$, melhorando para $3^k$ no caso especial de observáveis de produto tensorial.

Por outro lado, embora o ensemble de Clifford global possa ser usado para estimar observáveis não locais, sua profundidade de circuito não é favorável. As sombras gaussianas fermiônicas são úteis em simulações de física de muitos corpos, química e materiais, mas exigem circuitos matchgate, que podem ter um overhead significativo em dispositivos de uso geral baseados em portas.

As shallow shadows têm uma profundidade de circuito entre a de Pauli local e a de Clifford global, o que as torna uma variante muito promissora para implementação prática, embora ainda haja algumas questões em aberto sobre a profundidade de circuito ideal para um dado observável.

À medida que a fronteira da pesquisa em shadow tomography continua se expandindo, alguns dos problemas mais relevantes a serem abordados incluem escolher de forma ótima o ensemble de medição para um conjunto-alvo de observáveis, melhorar os limites de complexidade amostral e estender o protocolo a observáveis com evolução temporal. Ainda assim, a shadow tomography continua sendo um dos principais fluxos de trabalho para extrair o máximo de informação clássica de um estado quântico usando um número mínimo de medições, o que a torna ideal para caracterizar estados quânticos em escalas nas quais tanto a tomografia completa quanto a simulação clássica são impossíveis.

O ensemble de Pauli local

Este tutorial usa a base de Pauli aleatória, isto é, portas de Clifford de um qubit, para coletar sombras e estimar propriedades do estado GHZ. Observe que as medições de Pauli são uma ótima escolha para estimar observáveis locais (como densidade de spin, funções de correlação locais e densidade de probabilidade de posição), mas não observáveis globais (como energia total, momento e carga).

Aqui, uma base de Pauli aleatória $\{X, Y, Z\}$ por qubit é escolhida de forma independente e a rotação de um qubit correspondente é aplicada, seguida da medição:

$$\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 esse ensemble, o inverso do canal de medição é:

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

O fator $3$ vem da média sobre as três bases de Pauli, que produz um canal despolarizante com parâmetro $\frac{1}{3}$. O estimador completo então se fatora como um produto de estimadores de um qubit. Como mostra a tabela, para $M$ observáveis, cada um de localidade $k$, o número de sombras de Pauli necessárias para estimar todos os observáveis com erro aditivo $\varepsilon$ e alta probabilidade é $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. O tamanho do sistema não afeta o custo de medição, desde que os observáveis permaneçam locais, de modo que um número exponencialmente grande de observáveis pode ser estimado simultaneamente com custo apenas logarítmico.

Mediana das médias

Uma média simples sobre $T$ snapshots reduz a variância de uma estimativa de shot único de um observável como $\frac{1}{T}$, mas é sensível a outliers. A mediana das médias pode ser usada para remediar isso. O processo envolve o seguinte:

• Dividir os $T$ snapshots em $K$ grupos de mesmo tamanho
• Calcular a média dentro de cada grupo
• Retornar a mediana das médias dos grupos

Este tutorial executa a shadow tomography de Pauli local para estimar propriedades do estado GHZ, primeiro com um exemplo simulado de 10 qubits e depois demonstrando um experimento de hardware com 100 qubits.

Requisitos

Antes de começar este tutorial, certifique-se de ter o seguinte instalado:

• Qiskit SDK v2.0 ou posterior, com suporte a visualização

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

• Qiskit Paulice (pip install qiskit-paulice)

Configuração

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

As funções auxiliares a seguir são usadas pelo exemplo de hardware em larga escala mais adiante no tutorial:

• bad_cz, bad_readout, bad_coherence — Consultam os dados de calibração do backend para identificar portas, qubits e arestas que excedem limites de erro ou de coerência, para que possam ser excluídos do layout do circuito.

• parallel_ghz — Constrói um circuito GHZ de um dado tamanho usando uma busca em largura sobre o mapa de acoplamento, distribuindo camadas de CNOT em paralelo para minimizar a profundidade do circuito.

• weighted_coverage, active_wires, z_trace_backward — Calculam quanto da superfície de erro do circuito GHZ é protegida por um dado conjunto de verificações de paridade espaço-temporais, usando a técnica de propagação reversa de [6].

• add_measurements — Acrescenta medições dos qubits GHZ e dos qubits de verificação em registradores clássicos separados, para que possam ser pós-selecionados de forma independente.

• build_shadow_circuits — Produz um circuito por base de Pauli aleatória, acrescentando as rotações de um qubit apropriadas antes da medição.

• parse_hw_results_w_checks — Desempacota a saída do Sampler em triplas (ghz_bits, check_bits, bases), uma por shot.

• adapt_to_flat_bits, estimate_zz_nn_hw — Mapeiam os bits do registrador GHZ de volta para o espaço de índices do dispositivo completo e calculam a estimativa de mediana das médias 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

Exemplo em pequena escala com simulador

Primeiro, prepare um estado GHZ de 10 qubits usando simulação de vetor de estado e estime suas propriedades usando sombras de Pauli locais. O estado GHZ é um estado quântico maximamente emaranhado, com três ou mais qubits, que serve como benchmark de emaranhamento multipartite. Ele é um bloco de construção comum em protocolos de comunicação e criptografia quânticas. Além disso, o circuito GHZ tem propriedades analíticas totalmente conhecidas, o que o torna ideal para verificação neste tutorial.

Etapa 1: Mapear entradas clássicas para um problema quântico

Prepare um 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)

Etapa 2: Otimizar para execução em hardware quântico

Em seguida, defina uma função que coleta snapshots de sombra girando o estado GHZ para bases de Pauli aleatórias. Isso espelha o que a primitiva Sampler faz no hardware no exemplo em larga 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

Etapa 3: Executar usando primitives do Qiskit

Para esta simulação, use diretamente o vetor 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.

Etapa 4: Pós-processar e retornar o resultado no formato clássico desejado

Defina o estimador de shadow e o estimador median-of-means (MOM) para as bases de medição aleatórias usadas para coletar os snapshots de shadow. Use essas funções para estimar propriedades do estado GHZ, incluindo a magnetização de um único qubit, correlatores de dois qubits e a RDM de dois qubits, a partir das medições coletadas. Depois, compare as estimativas com os valores exatos.

_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

Visualize os resultados

O primeiro gráfico mostra os valores estimados de magnetização de um único qubit para os diferentes qubits. O segundo gráfico mostra a distribuição dos valores estimados do correlator de dois pontos para todos os pares de qubits. O último gráfico mostra a magnitude das entradas da matriz densidade reduzida de dois qubits. Uma coleção bastante grande de snapshots é usada para o tamanho de sistema dado, então os resultados ficam próximos dos exatos.

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

Avalie o escalonamento de shots

Neste tutorial, 5.000 snapshots são usados para estimar os observáveis. O escalonamento de shots da shadow tomography é avaliado plotando o erro médio de uma das funções estimadas em relação a diferentes orçamentos de shots. O erro no correlator de dois qubits decai como $1/\sqrt{T}$, como esperado para a convergência estatística padrão de um estimador de média. A importância da shadow tomography é que o prefator constante diante de $1/\sqrt{T}$ não cresce com o número de qubits para observáveis locais, como os examinados neste 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

Exemplo de hardware em grande escala

A simulação anterior confirma que a shadow tomography de Pauli local recupera os observáveis corretos em um vetor de estado sem ruído. Executar o mesmo protocolo com 100 qubits em hardware real introduz dois desafios que não aparecem na simulação: os erros de porta e de leitura se acumulam ao longo do circuito GHZ profundo, e encontrar o layout ótimo de qubits exige pesquisar o grafo de acoplamento do dispositivo em vez de escolher um layout manualmente.

Para lidar com o primeiro desafio, use detecção de erros de baixo overhead com códigos espaço-temporais [6]. Um pequeno número de qubits ancilla adjacentes à árvore GHZ é usado como verificações de paridade. Cada ancilla é emaranhada com dois de seus vizinhos GHZ antes da medição, de modo que um único erro de bit-flip em qualquer um dos vizinhos inverte o resultado da ancilla de 0 para 1. Snapshots em que qualquer ancilla de verificação lê 1 são descartados, deixando um subconjunto mais limpo a partir do qual estimar $\langle Z_i Z_{i+1}\rangle$.

O segundo desafio é tratado pesquisando layouts candidatos da árvore GHZ e selecionando aquele que maximiza a cobertura de verificações de paridade, ou seja, a fração de fios do circuito monitorados por pelo menos uma verificação, dentro de um aumento aceitável na profundidade do circuito.

Etapas 1-4 (executadas em um único bloco de código)

O mesmo fluxo de trabalho de quatro etapas é aplicado agora a um estado GHZ de 100 qubits em hardware quântico real, adicionando dois recursos específicos de hardware que não surgem na simulação:

• A Etapa 1 caracteriza o ruído do Backend (erros de porta, erros de leitura e tempos T2) para que qubits e arestas não confiáveis sejam excluídos do layout do circuito.

• A Etapa 2 pesquisa possíveis layouts de árvore GHZ para encontrar um que maximize a cobertura de verificações de paridade dentro de uma profundidade de circuito aceitável — substituindo a cadeia trivial usada na simulação.

• A Etapa 3 constrói os circuitos de shadow (um por base de Pauli aleatória), os transpila para portas nativas do hardware e os envia à primitiva Sampler com shots=8 (cada circuito produz oito bitstrings independentes).

• A Etapa 4 pós-seleciona os snapshots cujos qubits de verificação leem todos 0 (nenhum erro detectado) e, em seguida, estima $\langle Z_i Z_{i+1}\rangle$ com e sem essa pós-seleção para demonstrar a melhoria obtida com a detecção de erros.

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

Interprete os resultados

Os dois gráficos mostram $\langle Z_i Z_{i+1}\rangle$ estimado pela tomografia de shadow para cada par de qubits vizinhos na cadeia GHZ de 100 qubits, antes e depois da pós-seleção nas verificações de paridade.

Para um estado GHZ ideal, todos os correlatores entre vizinhos mais próximos são iguais a $1$. No hardware, o ruído reduz as estimativas brutas e introduz dispersão entre qubits. Após a detecção de erros, os snapshots em que ocorreu pelo menos um erro detectável são removidos, e as estimativas restantes se aproximam de $1$ com variância reduzida, o que é consistente com um estado efetivo mais limpo.

A melhoria após a detecção de erros é consistente com as verificações detectando erros reais.

Próximos passos

Explore os recursos a seguir para saber mais:

Recomendações

• Formulação geral da informação quântica
• Tutorial de detecção de erros de baixo overhead usando códigos espaço-temporais
• Documentação da API do add-on Qiskit Paulice

Referências

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