Shadow-Tomographie des Greenberger-Horne-Zeilinger-Zustands

Nutzungsschätzung: Unter 10 Sekunden auf einem Heron-r2-Prozessor (HINWEIS: Dies ist nur eine Schätzung. Deine Laufzeit kann abweichen.)

Lernergebnisse

• Wie die Shadow-Tomographie einen klassischen Schatten eines Quantenzustands konstruiert und warum sie exponentiell weniger Messungen benötigt als die vollständige Zustandstomographie

• Wie du lokale Pauli-Shadow-Tomographie mit Qiskit implementierst, um Ein-Qubit-Magnetisierungen, Zweipunkt-Korrelatoren und reduzierte Zwei-Qubit-Dichtematrizen zu schätzen

• Wie du qiskit-paulice für Fehlererkennung mit geringem Overhead mittels Spacetime-Codes in einem großskaligen Hardware-Experiment nutzt

Voraussetzungen

• Messungen

• Quantenzustandstomographie

• Fehlererkennung mit Spacetime-Codes

Hintergrund

Die vollständige Quantenzustandstomographie eines $N$-Qubit-Systems erfordert exponentiell viele Messungen, um alle $4^N -1$ unabhängigen Parameter der Dichtematrix zu schätzen. Für ein System mit 50 Qubits benötigt dieser Prozess zum Beispiel in der Größenordnung von $10^{30}$ Shots, mehr als das Milliardenfache der Anzahl der Sterne im beobachtbaren Universum. Die Shadow-Tomographie [1] ist ein quanteninformatisches Protokoll, das einen „klassischen Schatten“ des Quantenzustands konstruiert, aus dem Ziel-Observablen mit einer Shot-Anzahl geschätzt werden können, die logarithmisch mit der Anzahl der Observablen skaliert und unabhängig von der Systemgröße $N$ ist. Dadurch ist sie in vielen Bereichen des Quantencomputings relevant, darunter Energieschätzung, Berechnungen von Verschränkungsentropie und Korrelationsfunktionen, Kernel-Schätzung im Quanten-Machine-Learning und anderen.

Dieses Tutorial zeigt, wie du mit der Shadow-Tomographie lokale Observablen schätzt, etwa die Ein-Qubit-Magnetisierung $\langle Z_i \rangle$ für Qubit $i$, die Zweipunkt-Korrelationsfunktion $\langle Z_i Z_j\rangle$ und die reduzierte Zwei-Qubit-Dichtematrix (2-RDM) des stark verschränkten Greenberger-Horne-Zeilinger-Zustands (GHZ). Das Tutorial veranschaulicht zunächst einen kleinen Shadow-Tomographie-Workflow an einem 10-Qubit-GHZ-Zustand mittels Statevector-Simulation. Danach zeigt es eine 100-Qubit-Demonstration auf echter Quantenhardware.

Protokoll der Shadow-Tomographie

Das Protokoll der Shadow-Tomographie besteht aus den folgenden Schritten:

1. Wende eine zufällige Unitäre $U$, gezogen aus einem unitären Ensemble, auf den Anfangszustand an.
2. Miss das Ergebnis in der Rechenbasis und erhalte einen Bitstring $\ket{b}$.
3. Speichere den klassischen Snapshot $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, wobei $\mathcal{M}^{-1}$ die Inverse des Messkanals $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$ ist.
4. Wiederhole den Prozess $T$-mal, wobei die Sammlung der $T$ Snapshots $\{\hat{\rho}_t\}$ einen klassischen Schatten $S(\rho, T)$ bildet.
5. Schätze die Observable $o$ des Operators $O$ als $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Varianten der Shadow-Tomographie

Die Wahl des unitären Ensembles, also der Basis, in Schritt 1 des Protokolls beeinflusst die Stichprobenkomplexität, die Varianz der Schätzungen und die Informationsvollständigkeit einer gegebenen Shadow-Tomographie. In der Literatur wurden viele Ensemble-Gruppen betrachtet, jeweils mit Kompromissen zwischen Komplexität und Effizienz der Observablen-Schätzung. Die Eigenschaften der am häufigsten verwendeten Varianten sind in der folgenden Tabelle zusammengefasst, wobei diese Liste nicht erschöpfend ist.

Variante	Ensemble	Stichprobenkomplexität	Circuit-Tiefe	Geeignete Observablen
lokale Pauli [1]	Ein-Qubit-Paulis $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	lokal/wenige Körper
globale Clifford [1]	$N$-Qubit-Clifford	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	niedriger Rang, Fidelity
fermionische Gauß [2], [3]	Matchgate-Circuits	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermionisch
Shallow Shadows [4], [5]	Circuits geringer Tiefe	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	nicht-lokal (hardwarefreundlich)

Die in der Tabelle verwendeten Begriffe sind wie folgt definiert:

• „Stichprobenkomplexität“ bezeichnet die Anzahl der Messungen, die für eine gegebene Variante der Shadow-Tomographie benötigt werden.

• „Circuit-Tiefe“ beschreibt, wie aufwendig die Messungen auf Quantenhardware zu implementieren sind.

• $N$ ist die Anzahl der Qubits.

• $k$ ist die Lokalität der Observablen (zum Beispiel hat $\langle Z_i \rangle$ den Wert $k=1$ und $\langle Z_i Z_j\rangle$ den Wert $k=2$).

• $M$ ist die Anzahl der gleichzeitig zu schätzenden Observablen.

• $\varepsilon$ ist die additive Fehlertoleranz.

• $\mu$ ist die Anzahl der Fermionen.

Obwohl das lokale Pauli-Ensemble am hardwarefreundlichsten ist, zahlt es einen Aufschlag von $3^k$ bei der Stichprobenkomplexität, wodurch es für die Schätzung nicht-lokaler Eigenschaften schlecht geeignet ist. Allgemeiner skaliert die Schranke der Stichprobenkomplexität für lokale Shadows wie $4^k$ und verbessert sich im Spezialfall von Tensorprodukt-Observablen auf $3^k$.

Umgekehrt kann das globale Clifford-Ensemble zwar zur Schätzung nicht-lokaler Observablen verwendet werden, aber seine Circuit-Tiefe ist ungünstig. Fermionische Gaußsche Shadows sind nützlich bei Simulationen von Vielteilchenphysik, Chemie und Materialien, erfordern aber Matchgate-Circuits, die auf universellen gatterbasierten Geräten erheblichen Overhead verursachen können.

Shallow Shadows haben eine Circuit-Tiefe zwischen lokalem Pauli und globalem Clifford, was sie zu einer sehr vielversprechenden Variante für die praktische Umsetzung macht, auch wenn es noch offene Fragen zur optimalen Circuit-Tiefe für eine gegebene Observable gibt.

Da die Forschungsfront der Shadow-Tomographie weiter wächst, gehören zu den relevantesten zu lösenden Problemen die optimale Wahl des Messensembles für eine Zielmenge von Observablen, die Verbesserung der Schranken der Stichprobenkomplexität und die Erweiterung des Protokolls auf zeitentwickelte Observablen. Dennoch bietet die Shadow-Tomographie weiterhin einen der führenden Workflows, um mit einer minimalen Anzahl von Messungen maximale klassische Information aus einem Quantenzustand zu gewinnen. Das macht sie ideal zur Charakterisierung von Quantenzuständen in Größenordnungen, in denen sowohl die vollständige Tomographie als auch die klassische Simulation unmöglich sind.

Das lokale Pauli-Ensemble

Dieses Tutorial verwendet die zufällige Pauli-Basis, also Ein-Qubit-Clifford-Gates, um Shadows zu sammeln und Eigenschaften des GHZ-Zustands zu schätzen. Beachte, dass Pauli-Messungen eine hervorragende Wahl zur Schätzung lokaler Observablen sind (wie Spindichte, lokale Korrelationsfunktionen und Ortswahrscheinlichkeitsdichte), aber nicht globaler (wie Gesamtenergie, Impuls und Ladung).

Hier wird pro Qubit unabhängig eine zufällige Pauli-Basis $\{X, Y, Z\}$ gewählt und die entsprechende Ein-Qubit-Rotation angewendet, gefolgt von der Messung:

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

Für dieses Ensemble ist die Inverse des Messkanals:

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

Der Faktor $3$ ergibt sich aus der Mittelung über die drei Pauli-Basen, die einen depolarisierenden Kanal mit Parameter $\frac{1}{3}$ erzeugt. Der vollständige Schätzer zerfällt dann in ein Produkt von Ein-Qubit-Schätzern. Wie in der Tabelle gezeigt, beträgt für $M$ Observablen mit jeweils Lokalität $k$ die Anzahl der Pauli-Shadows, die nötig ist, um alle Observablen mit hoher Wahrscheinlichkeit bis auf den additiven Fehler $\varepsilon$ zu schätzen, $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. Die Systemgröße beeinflusst die Messkosten nicht, solange die Observablen lokal bleiben; exponentiell viele Observablen können also gleichzeitig mit nur logarithmischen Kosten geschätzt werden.

Median of Means

Ein einfacher Mittelwert über $T$ Snapshots reduziert die Varianz einer Einzel-Shot-Schätzung einer Observablen wie $\frac{1}{T}$, ist aber empfindlich gegenüber Ausreißern. Abhilfe kann Median of Means schaffen. Der Prozess umfasst Folgendes:

• Die $T$ Snapshots in $K$ gleich große Gruppen aufteilen
• Den Mittelwert innerhalb jeder Gruppe berechnen
• Den Median der Gruppenmittelwerte zurückgeben

Dieses Tutorial führt lokale Pauli-Shadow-Tomographie durch, um Eigenschaften des GHZ-Zustands zu schätzen, zunächst an einem simulierten 10-Qubit-Beispiel und anschließend mit einem 100-Qubit-Hardware-Experiment.

Anforderungen

Stelle vor dem Start dieses Tutorials sicher, dass Folgendes installiert ist:

• Qiskit SDK v2.0 oder höher, mit Unterstützung für Visualisierung

• IBM Quantum Compute v0.28 oder höher (pip install qiskit-ibm-runtime)

• Qiskit Paulice (pip install qiskit-paulice)

Einrichtung

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

Die folgenden Hilfsfunktionen werden vom großskaligen Hardware-Beispiel später im Tutorial verwendet:

• bad_cz, bad_readout, bad_coherence — Fragen die Kalibrierdaten des Backends ab, um Gates, Qubits und Kanten zu identifizieren, die Fehler- oder Kohärenzschwellen überschreiten, damit sie aus dem Circuit-Layout ausgeschlossen werden können.

• parallel_ghz — Baut einen GHZ-Circuit gegebener Größe mittels Breitensuche über die Kopplungskarte und verteilt CNOT-Schichten parallel, um die Circuit-Tiefe zu minimieren.

• weighted_coverage, active_wires, z_trace_backward — Berechnen, wie viel der Fehleroberfläche des GHZ-Circuits durch eine gegebene Menge von Spacetime-Paritätsprüfungen geschützt ist, mit der Rückwärtspropagations-Technik aus [6].

• add_measurements — Hängt Messungen der GHZ-Qubits und der Prüf-Qubits an separate klassische Register an, damit sie unabhängig voneinander nachselektiert werden können.

• build_shadow_circuits — Erzeugt einen Circuit pro zufälliger Pauli-Basis, indem vor der Messung die passenden Ein-Qubit-Rotationen angehängt werden.

• parse_hw_results_w_checks — Entpackt die Sampler-Ausgabe in (ghz_bits, check_bits, bases)-Tripel, eines pro Shot.

• adapt_to_flat_bits, estimate_zz_nn_hw — Bilden die Bits des GHZ-Registers zurück in den Indexraum des gesamten Geräts ab und berechnen die Median-of-Means-Schätzung von $\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

Kleinskaliges Simulator-Beispiel

Bereite zunächst mit Statevector-Simulation einen 10-Qubit-GHZ-Zustand vor und schätze seine Eigenschaften mit lokalen Pauli-Shadows. Der GHZ-Zustand ist ein maximal verschränkter Quantenzustand mit drei oder mehr Qubits, der als Benchmark für Mehrparteien-Verschränkung dient. Er ist ein gängiger Baustein für Quantenkommunikations- und Kryptographieprotokolle. Außerdem hat der GHZ-Circuit vollständig bekannte analytische Eigenschaften, was ihn ideal für die Verifikation in diesem Tutorial macht.

Schritt 1: Klassische Eingaben auf ein Quantenproblem abbilden

Bereite einen idealen 10-Qubit-GHZ-Zustand vor.

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)

Schritt 2: Für die Ausführung auf Quantenhardware optimieren

Definiere als Nächstes eine Funktion, die Shadow-Snapshots sammelt, indem sie den GHZ-Zustand in zufällige Pauli-Basen rotiert. Das spiegelt wider, was die Sampler-Primitive im großskaligen Beispiel auf der Hardware tut.

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

Schritt 3: Mit Qiskit-Primitives ausführen

Für diese Simulation verwenden wir direkt den Statevector.

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.

Schritt 4: Nachbearbeiten und das Ergebnis im gewünschten klassischen Format zurückgeben

Definiere den Shadow-Estimator und den Median-of-Means-Estimator (MOM) für die zufälligen Messbasen, die zum Sammeln der Shadow-Snapshots verwendet werden. Verwende diese Funktionen, um aus den gesammelten Messungen Eigenschaften des GHZ-Zustands zu schätzen, darunter die Magnetisierung einzelner Qubits, Zwei-Qubit-Korrelatoren und die Zwei-Qubit-RDM. Vergleiche anschließend die Schätzwerte mit den exakten Werten.

_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

Ergebnisse visualisieren

Die erste Grafik zeigt die geschätzten Magnetisierungswerte einzelner Qubits für die verschiedenen Qubits. Die zweite Grafik zeigt die Verteilung der geschätzten Zwei-Punkt-Korrelatorwerte für alle Qubit-Paare. Die letzte Grafik zeigt den Betrag der Einträge der reduzierten Zwei-Qubit-Dichtematrix. Für die gegebene Systemgröße wird eine recht große Menge an Snapshots verwendet, daher liegen die Ergebnisse nahe an den exakten Werten.

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

Shot-Skalierung bewerten

In diesem Tutorial werden 5.000 Snapshots verwendet, um die Observablen zu schätzen. Die Shot-Skalierung der Shadow-Tomographie wird bewertet, indem der mittlere Fehler einer der geschätzten Funktionen gegen verschiedene Shot-Budgets aufgetragen wird. Der Fehler im Zwei-Qubit-Korrelator nimmt wie $1/\sqrt{T}$ ab, wie es für die übliche statistische Konvergenz eines Mittelwertschätzers zu erwarten ist. Die Bedeutung der Shadow-Tomographie liegt darin, dass der konstante Vorfaktor vor $1/\sqrt{T}$ für lokale Observablen, wie sie in diesem Tutorial untersucht werden, nicht mit der Anzahl der Qubits wächst.

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

Beispiel für Hardware im großen Maßstab

Die vorangegangene Simulation bestätigt, dass die lokale Pauli-Shadow-Tomographie auf einem rauschfreien Statevector die korrekten Observablen liefert. Wird dasselbe Protokoll mit 100 Qubits auf echter Hardware ausgeführt, entstehen zwei Herausforderungen, die in der Simulation nicht auftreten: Gate- und Auslesefehler summieren sich über den tiefen GHZ-Circuit auf, und das Finden des optimalen Qubit-Layouts erfordert eine Suche im Kopplungsgraphen des Geräts, anstatt ein Layout von Hand zu wählen.

Um die erste Herausforderung zu bewältigen, verwende eine Fehlererkennung mit geringem Overhead durch Spacetime-Codes [6]. Eine kleine Anzahl von Ancilla-Qubits neben dem GHZ-Baum dient als Paritätsprüfungen. Jede Ancilla wird vor der Messung mit zwei ihrer GHZ-Nachbarn verschränkt, sodass ein einzelner Bitflip-Fehler bei einem der Nachbarn das Ergebnis der Ancilla von 0 auf 1 kippt. Snapshots, bei denen eine Prüf-Ancilla 1 liefert, werden verworfen. Übrig bleibt eine sauberere Teilmenge, aus der $\langle Z_i Z_{i+1}\rangle$ geschätzt wird.

Die zweite Herausforderung wird gelöst, indem mögliche GHZ-Baum-Layouts durchsucht und dasjenige ausgewählt wird, das die Abdeckung der Paritätsprüfungen maximiert, also den Anteil der Circuit-Leitungen, die von mindestens einer Prüfung überwacht werden, bei einem akzeptablen Anstieg der Circuit-Tiefe.

Schritte 1-4 (in einem einzigen Codeblock ausgeführt)

Derselbe vierstufige Workflow wird nun auf einen 100-Qubit-GHZ-Zustand auf echter Quantenhardware angewendet. Dabei kommen zwei hardwarespezifische Merkmale hinzu, die in der Simulation nicht auftreten:

• Schritt 1 charakterisiert das Rauschen des Backends (Gate-Fehler, Auslesefehler und T2-Zeiten), damit unzuverlässige Qubits und Kanten aus dem Circuit-Layout ausgeschlossen werden.

• Schritt 2 durchsucht mögliche GHZ-Baum-Layouts, um eines zu finden, das die Abdeckung der Paritätsprüfungen bei akzeptabler Circuit-Tiefe maximiert – und ersetzt damit die in der Simulation verwendete einfache Kette.

• Schritt 3 erstellt die Shadow-Circuits (eine pro zufälliger Pauli-Basis), transpiliert sie in hardwarenative Gates und übergibt sie mit shots=8 an das Sampler-Primitive (jeder Circuit erzeugt acht unabhängige Bitstrings).

• Schritt 4 wählt per Postselektion die Snapshots aus, bei denen alle Prüf-Qubits 0 lesen (kein erkannter Fehler), und schätzt dann $\langle Z_i Z_{i+1}\rangle$ mit und ohne diese Postselektion, um die Verbesserung durch die Fehlererkennung zu zeigen.

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

Ergebnisse interpretieren

Die beiden Grafiken zeigen $\langle Z_i Z_{i+1}\rangle$, geschätzt aus der Shadow-Tomographie für jedes benachbarte Qubit-Paar in der 100-Qubit-GHZ-Kette, vor und nach der Postselektion auf die Paritätsprüfungen.

Für einen idealen GHZ-Zustand sind alle Nächste-Nachbar-Korrelatoren gleich $1$. Auf der Hardware drückt Rauschen die Rohschätzwerte nach unten und führt zu Streuung von Qubit zu Qubit. Nach der Fehlererkennung werden Snapshots entfernt, in denen mindestens ein erkennbarer Fehler aufgetreten ist, und die verbleibenden Schätzwerte rücken mit geringerer Varianz näher an $1$, was mit einem saubereren effektiven Zustand übereinstimmt.

Die Verbesserung nach der Fehlererkennung stimmt damit überein, dass die Prüfungen echte Fehler erkennen.

Nächste Schritte

Weitere Ressourcen zum Vertiefen:

Empfehlungen

• General formulation of quantum information
• Low-overhead error detection using spacetime codes tutorial
• Qiskit Paulice add-on API documentation

Referenzen

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