Tomografia cieni stanu Greenbergera-Horne'a-Zeilingera

Szacowane zużycie: poniżej 10 sekund na procesorze Heron r2 (UWAGA: To tylko szacunek. Twój czas działania może się różnić.)

Cele nauczania

• Jak tomografia cieni konstruuje klasyczny cień stanu kwantowego i dlaczego wymaga wykładniczo mniej pomiarów niż pełna tomografia stanu

• Jak zaimplementować lokalną tomografię cieni Pauliego z użyciem Qiskit, aby oszacować jednokubitowe magnetyzacje, korelatory dwupunktowe i dwukubitowe zredukowane macierze gęstości

• Jak używać qiskit-paulice do wykrywania błędów o niskim narzucie z użyciem kodów czasoprzestrzennych w eksperymencie sprzętowym na dużą skalę

Wymagania wstępne

• Pomiary

• Tomografia stanu kwantowego

• Wykrywanie błędów z użyciem kodów czasoprzestrzennych

Tło

Pełna tomografia stanu kwantowego układu $N$ kubitów wymaga wykładniczo wielu pomiarów, aby oszacować wszystkie $4^N -1$ niezależnych parametrów macierzy gęstości. Dla układu 50 kubitów, na przykład, proces ten wymaga rzędu $10^{30}$ shotów, czyli ponad miliard razy więcej niż liczba gwiazd w obserwowalnym wszechświecie. Tomografia cieni [1] to protokół informacji kwantowej, który konstruuje „klasyczny cień” stanu kwantowego, z którego można szacować docelowe obserwable, używając liczby shotów skalującej się logarytmicznie z liczbą obserwabli i niezależnej od rozmiaru układu $N$, co czyni go istotnym w wielu dziedzinach obliczeń kwantowych, w tym w szacowaniu energii, obliczeniach entropii splątania i funkcji korelacji, estymacji jąder w kwantowym uczeniu maszynowym i innych.

Ten samouczek pokazuje, jak użyć tomografii cieni do oszacowania lokalnych obserwabli, takich jak jednokubitowa magnetyzacja $\langle Z_i \rangle$ dla kubitu $i$, dwupunktowa funkcja korelacji $\langle Z_i Z_j\rangle$ oraz dwukubitowa zredukowana macierz gęstości (2-RDM) silnie splątanego stanu Greenbergera-Horne'a-Zeilingera (GHZ). Samouczek najpierw ilustruje małoskalowy przepływ pracy tomografii cieni na 10-kubitowym stanie GHZ z użyciem symulacji wektora stanu. Następnie pokazuje 100-kubitową demonstrację na prawdziwym sprzęcie kwantowym.

Protokół tomografii cieni

Protokół tomografii cieni składa się z następujących kroków:

1. Zastosuj losową unitarną $U$, wylosowaną z zespołu unitarnego, do początkowego stanu kwantowego.
2. Zmierz wynik w bazie obliczeniowej, uzyskując ciąg bitów $\ket{b}$.
3. Zapisz klasyczną migawkę, $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, gdzie $\mathcal{M}^{-1}$ jest odwrotnością kanału pomiaru $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Powtórz proces $T$ razy, gdzie zbiór $T$ migawek $\{\hat{\rho}_t\}$ tworzy klasyczny cień $S(\rho, T)$.
5. Oszacuj obserwablę $o$ operatora $O$ jako $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Warianty tomografii cieni

Wybór zespołu unitarnego, czyli bazy, w kroku 1 protokołu wpływa na złożoność próbkowania, wariancję oszacowań i kompletność informacji danej tomografii cieni. W literaturze rozważano wiele grup zespołów, każda z kompromisami między złożonością a efektywnością szacowania obserwabli. Właściwości najczęściej używanych wariantów podsumowano w poniższej tabeli, choć lista ta nie jest wyczerpująca.

Wariant	Zespół	Złożoność próbkowania	Głębokość obwodu	Odpowiednie obserwable
lokalny Pauli [1]	jednokubitowe operatory Pauliego $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	lokalne/kilkociałowe
globalny Clifford [1]	$N$-kubitowy Clifford	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	niskiego rzędu, wierność
fermionowy gaussowski [2], [3]	obwody matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermionowe
płytkie cienie [4], [5]	obwody o małej głębokości	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	nielokalne (przyjazne sprzętowi)

Terminy użyte w tabeli są zdefiniowane następująco:

• „złożoność próbkowania” oznacza liczbę pomiarów potrzebnych dla danego wariantu tomografii cieni.

• „głębokość obwodu” opisuje, jak kosztowna jest implementacja pomiarów na sprzęcie kwantowym.

• $N$ to liczba kubitów.

• $k$ to lokalność obserwabli (na przykład $\langle Z_i \rangle$ ma $k=1$, a $\langle Z_i Z_j\rangle$ ma $k=2$).

• $M$ to liczba obserwabli szacowanych jednocześnie.

• $\varepsilon$ to addytywna tolerancja błędu.

• $\mu$ to liczba fermionów.

Choć lokalny zespół Pauliego jest najbardziej przyjazny sprzętowi, wiąże się z karą złożoności próbkowania $3^k$, co czyni go niezbyt odpowiednim do szacowania właściwości nielokalnych. Ogólniej, granica złożoności próbkowania dla cieni lokalnych skaluje się jak $4^k$, poprawiając się do $3^k$ w szczególnym przypadku obserwabli o strukturze iloczynu tensorowego.

Z kolei choć globalny zespół Clifforda może służyć do szacowania obserwabli nielokalnych, jego głębokość obwodu jest niekorzystna. Fermionowe cienie gaussowskie są użyteczne w symulacjach fizyki wielu ciał, chemii i materiałów, ale wymagają obwodów matchgate, które mogą wiązać się ze znacznym narzutem na urządzeniach bramkowych ogólnego przeznaczenia.

Płytkie cienie mają głębokość obwodu pomiędzy lokalnym Pauliego a globalnym Clifforda, co czyni je bardzo obiecującym wariantem do praktycznej implementacji, choć nadal istnieją otwarte pytania dotyczące optymalnej głębokości obwodu dla danej obserwabli.

W miarę jak front badań nad tomografią cieni nadal się rozszerza, do najistotniejszych problemów do rozwiązania należą optymalny dobór zespołu pomiarowego dla docelowego zbioru obserwabli, poprawa granic złożoności próbkowania oraz rozszerzenie protokołu na obserwable ewoluujące w czasie. Mimo to tomografia cieni nadal stanowi jeden z wiodących przepływów pracy do wydobywania maksymalnej klasycznej informacji ze stanu kwantowego przy minimalnej liczbie pomiarów, co czyni ją idealną do charakteryzowania stanów kwantowych w skali, w której niemożliwa jest zarówno pełna tomografia, jak i symulacja klasyczna.

Lokalny zespół Pauliego

Ten samouczek używa losowej bazy Pauliego, czyli jednokubitowych bramek Clifforda, do zbierania cieni i szacowania właściwości stanu GHZ. Zauważ, że pomiary Pauliego są świetnym wyborem do szacowania obserwabli lokalnych (takich jak gęstość spinu, lokalne funkcje korelacji i gęstość prawdopodobieństwa położenia), ale nie globalnych (takich jak energia całkowita, pęd i ładunek).

Tutaj niezależnie wybierana jest losowa baza Pauliego $\{X, Y, Z\}$ dla każdego kubitu i stosowany jest odpowiadający jej obrót jednokubitowy, po którym następuje pomiar:

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

Dla tego zespołu odwrotność kanału pomiaru wynosi:

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

Czynnik $3$ pochodzi z uśrednienia po trzech bazach Pauliego, co daje kanał depolaryzujący z parametrem $\frac{1}{3}$. Pełny estymator faktoryzuje się wówczas na iloczyn estymatorów jednokubitowych. Jak pokazano w tabeli, dla $M$ obserwabli, każdej o lokalności $k$, liczba cieni Pauliego potrzebnych do oszacowania wszystkich obserwabli z addytywnym błędem $\varepsilon$ z wysokim prawdopodobieństwem wynosi $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. Rozmiar układu nie wpływa na koszt pomiarowy, dopóki obserwable pozostają lokalne, więc wykładniczo wiele obserwabli można szacować jednocześnie przy jedynie logarytmicznym koszcie.

Mediana ze średnich

Zwykła średnia z $T$ migawek zmniejsza wariancję jednostrzałowego oszacowania obserwabli jak $\frac{1}{T}$, ale jest wrażliwa na wartości odstające. Można temu zaradzić, stosując medianę ze średnich. Proces obejmuje następujące kroki:

• Podział $T$ migawek na $K$ grup o równej liczebności
• Obliczenie średniej w każdej grupie
• Zwrócenie mediany ze średnich grupowych

Ten samouczek wykonuje lokalną tomografię cieni Pauliego, aby oszacować właściwości stanu GHZ, najpierw na symulowanym przykładzie 10-kubitowym, a następnie demonstrując 100-kubitowy eksperyment sprzętowy.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane:

• Qiskit SDK w wersji 2.0 lub nowszej, z obsługą wizualizacji

• IBM Quantum Compute w wersji 0.28 lub nowszej (pip install qiskit-ibm-runtime)

• Qiskit Paulice (pip install qiskit-paulice)

Konfiguracja

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

Poniższe funkcje pomocnicze są używane w przykładzie sprzętowym na dużą skalę w dalszej części samouczka:

• bad_cz, bad_readout, bad_coherence — Odpytują dane kalibracyjne backendu, aby zidentyfikować bramki, kubity i krawędzie przekraczające progi błędu lub koherencji, tak aby można je było wykluczyć z układu obwodu.

• parallel_ghz — Buduje obwód GHZ o zadanym rozmiarze, używając przeszukiwania wszerz po mapie sprzężeń i rozkładając warstwy CNOT równolegle, aby zminimalizować głębokość obwodu.

• weighted_coverage, active_wires, z_trace_backward — Obliczają, jaka część powierzchni błędów obwodu GHZ jest chroniona przez dany zbiór czasoprzestrzennych kontroli parzystości, z użyciem techniki propagacji wstecznej z [6].

• add_measurements — Dołącza pomiary kubitów GHZ i kubitów kontrolnych do oddzielnych rejestrów klasycznych, tak aby można je było postselekcjonować niezależnie.

• build_shadow_circuits — Tworzy jeden obwód dla każdej losowej bazy Pauliego, dołączając odpowiednie obroty jednokubitowe przed pomiarem.

• parse_hw_results_w_checks — Rozpakowuje wyjście Samplera do trójek (ghz_bits, check_bits, bases), po jednej na shot.

• adapt_to_flat_bits, estimate_zz_nn_hw — Odwzorowują bity rejestru GHZ z powrotem na pełną przestrzeń indeksów urządzenia i obliczają oszacowanie $\langle Z_i Z_{i+1}\rangle$ metodą mediany ze średnich.

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

Małoskalowy przykład na symulatorze

Najpierw przygotuj 10-kubitowy stan GHZ z użyciem symulacji wektora stanu i oszacuj jego właściwości za pomocą lokalnych cieni Pauliego. Stan GHZ to maksymalnie splątany stan kwantowy z trzema lub więcej kubitami, służący jako benchmark splątania wielocząstkowego. Jest powszechnym elementem składowym protokołów komunikacji kwantowej i kryptografii. Dodatkowo obwód GHZ ma w pełni znane właściwości analityczne, co czyni go idealnym do weryfikacji w tym samouczku.

Krok 1: Odwzoruj klasyczne dane wejściowe na problem kwantowy

Przygotuj idealny 10-kubitowy stan GHZ.

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)

Krok 2: Zoptymalizuj pod kątem wykonania na sprzęcie kwantowym

Następnie zdefiniuj funkcję zbierającą migawki cieni przez obracanie stanu GHZ do losowych baz Pauliego. Odzwierciedla to to, co robi prymityw Sampler na sprzęcie w przykładzie na dużą skalę.

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

Krok 3: Wykonaj za pomocą prymitywów Qiskit

W tej symulacji użyj bezpośrednio wektora stanu.

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.

Krok 4: Przetwórz końcowo i zwróć wynik w żądanym formacie klasycznym

Zdefiniuj estymator cieni oraz estymator mediany ze średnich (MOM) dla losowych baz pomiarowych użytych do zebrania migawek cieni. Użyj tych funkcji do oszacowania właściwości stanu GHZ, w tym magnetyzacji pojedynczego kubitu, korelatorów dwukubitowych i dwukubitowej zredukowanej macierzy gęstości (RDM) na podstawie zebranych pomiarów. Następnie porównaj oszacowania z dokładnymi wartościami.

_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

Wizualizacja wyników

Pierwszy wykres przedstawia oszacowane wartości magnetyzacji pojedynczego kubitu dla poszczególnych kubitów. Drugi wykres pokazuje rozkład oszacowanych wartości korelatorów dwupunktowych dla wszystkich par kubitów. Ostatni wykres pokazuje wielkość elementów dwukubitowej zredukowanej macierzy gęstości. Dla danego rozmiaru układu użyto dość dużego zbioru migawek, więc wyniki są bliskie dokładnym.

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

Ocena skalowania liczby strzałów

W tym samouczku do oszacowania obserwabli użyto 5000 migawek. Skalowanie liczby strzałów w tomografii cieni oceniamy, wykreślając średni błąd jednej z oszacowanych funkcji względem różnych budżetów strzałów. Błąd korelatora dwukubitowego maleje jak $1/\sqrt{T}$, zgodnie z oczekiwaniami dla standardowej zbieżności statystycznej estymatora średniej. Znaczenie tomografii cieni polega na tym, że stały współczynnik przy $1/\sqrt{T}$ nie rośnie wraz z liczbą kubitów dla obserwabli lokalnych, takich jak te badane w tym samouczku.

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

Przykład na dużą skalę na sprzęcie

Powyższa symulacja potwierdza, że lokalna tomografia cieni Pauliego odtwarza poprawne obserwable na bezszumowym wektorze stanu. Uruchomienie tego samego protokołu na 100 kubitach na prawdziwym sprzęcie wprowadza dwa wyzwania, które nie występują w symulacji: błędy bramek i odczytu kumulują się w głębokim obwodzie GHZ, a znalezienie optymalnego układu kubitów wymaga przeszukania grafu połączeń urządzenia zamiast ręcznego wyboru układu.

Aby sprostać pierwszemu wyzwaniu, użyj detekcji błędów o niskim narzucie z kodami czasoprzestrzennymi [6]. Niewielka liczba kubitów pomocniczych sąsiadujących z drzewem GHZ służy jako kontrole parzystości. Każdy kubit pomocniczy jest splątywany z dwoma sąsiadującymi kubitami GHZ przed pomiarem, więc pojedynczy błąd odwrócenia bitu na którymkolwiek z sąsiadów zmienia wynik kubitu pomocniczego z 0 na 1. Migawki, w których dowolny pomocniczy kubit kontrolny daje 1, są odrzucane, co pozostawia czystszy podzbiór do oszacowania $\langle Z_i Z_{i+1}\rangle$.

Drugie wyzwanie rozwiązujemy, przeszukując kandydackie układy drzewa GHZ i wybierając ten, który maksymalizuje pokrycie kontrolami parzystości, czyli odsetek linii obwodu monitorowanych przez co najmniej jedną kontrolę, przy akceptowalnym wzroście głębokości obwodu.

Kroki 1-4 (wykonywane w jednym bloku kodu)

Ten sam czterokrokowy przepływ pracy stosujemy teraz do 100-kubitowego stanu GHZ na prawdziwym sprzęcie kwantowym, dodając dwie cechy specyficzne dla sprzętu, które nie występują w symulacji:

• Krok 1 charakteryzuje szum Backendu (błędy bramek, błędy odczytu i czasy T2), tak aby zawodne kubity i krawędzie zostały wykluczone z układu obwodu.

• Krok 2 przeszukuje możliwe układy drzewa GHZ, aby znaleźć taki, który maksymalizuje pokrycie kontrolami parzystości przy akceptowalnej głębokości obwodu — zastępując trywialny łańcuch użyty w symulacji.

• Krok 3 buduje obwody cieni (po jednym na każdą losową bazę Pauliego), transpiluje je do bramek natywnych dla sprzętu i przesyła do prymitywu Sampler z shots=8 (każdy obwód daje osiem niezależnych ciągów bitów).

• Krok 4 postselekcjonuje migawki, w których wszystkie kubity kontrolne odczytują 0 (nie wykryto błędu), a następnie szacuje $\langle Z_i Z_{i+1}\rangle$ z postselekcją i bez niej, aby pokazać poprawę wynikającą z detekcji błędów.

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

Interpretacja wyników

Dwa wykresy pokazują $\langle Z_i Z_{i+1}\rangle$ oszacowane z tomografii cieni dla każdej pary sąsiednich kubitów w 100-kubitowym łańcuchu GHZ, przed postselekcją według kontroli parzystości i po niej.

Dla idealnego stanu GHZ wszystkie korelatory najbliższych sąsiadów są równe $1$. Na sprzęcie szum obniża surowe oszacowania i wprowadza rozrzut między kubitami. Po detekcji błędów migawki, w których wystąpił co najmniej jeden wykrywalny błąd, są usuwane, a pozostałe oszacowania zbliżają się do $1$ ze zmniejszoną wariancją, co jest zgodne z czystszym stanem efektywnym.

Poprawa po detekcji błędów jest zgodna z tym, że kontrole wykrywają rzeczywiste błędy.

Następne kroki

Aby dowiedzieć się więcej, zapoznaj się z poniższymi zasobami:

Zalecenia

• Ogólne sformułowanie informacji kwantowej
• Samouczek: detekcja błędów o niskim narzucie z użyciem kodów czasoprzestrzennych
• Dokumentacja API dodatku Qiskit Paulice

Bibliografia

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