Tomografi bayangan bagi keadaan Greenberger-Horne-Zeilinger

Anggaran penggunaan: Kurang daripada 10 saat pada pemproses Heron r2 (NOTA: Ini hanyalah anggaran. Masa jalan anda mungkin berbeza.)

Hasil pembelajaran

• Bagaimana tomografi bayangan membina bayangan klasik bagi keadaan kuantum dan mengapa ia memerlukan pengukuran yang jauh lebih sedikit secara eksponen berbanding tomografi keadaan penuh

• Cara melaksanakan tomografi bayangan Pauli tempatan menggunakan Qiskit untuk menganggar magnetisasi qubit tunggal, korelator dua titik, dan matriks ketumpatan terturun dua-qubit

• Cara menggunakan qiskit-paulice untuk pengesanan ralat beroverhead rendah menggunakan kod ruang-masa dalam eksperimen perkakasan berskala besar

Prasyarat

• Pengukuran

• Tomografi keadaan kuantum

• Pengesanan ralat menggunakan kod ruang-masa

Latar belakang

Tomografi keadaan kuantum penuh bagi sistem $N$-qubit memerlukan pengukuran yang banyak secara eksponen untuk menganggar semua $4^N -1$ parameter bebas bagi matriks ketumpatan. Bagi sistem 50 qubit, sebagai contoh, proses ini memerlukan kira-kira $10^{30}$ shot, iaitu lebih daripada satu bilion kali bilangan bintang di alam semesta yang boleh diperhatikan. Tomografi bayangan [1] ialah protokol maklumat kuantum yang membina “bayangan klasik” bagi keadaan kuantum, yang daripadanya observable sasaran boleh dianggar menggunakan bilangan shot yang berskala secara logaritma dengan bilangan observable dan tidak bergantung pada saiz sistem $N$, menjadikannya relevan merentasi banyak bidang dalam pengkomputeran kuantum, termasuk anggaran tenaga, pengiraan entropi keterkaitan dan fungsi korelasi, anggaran kernel dalam pembelajaran mesin kuantum, dan lain-lain.

Tutorial ini menunjukkan cara menggunakan tomografi bayangan untuk menganggar observable tempatan, seperti magnetisasi qubit tunggal $\langle Z_i \rangle$ bagi qubit $i$, fungsi korelasi dua titik $\langle Z_i Z_j\rangle$, dan matriks ketumpatan terturun dua-qubit (2-RDM), bagi keadaan Greenberger-Horne-Zeilinger (GHZ) yang sangat terkait. Tutorial ini mula-mula menggambarkan aliran kerja tomografi bayangan berskala kecil pada keadaan GHZ 10-qubit menggunakan simulasi statevector. Kemudian ia menunjukkan demonstrasi 100-qubit pada perkakasan kuantum sebenar.

Protokol tomografi bayangan

Protokol tomografi bayangan mempunyai langkah berikut:

1. Gunakan unitari rawak $U$, yang diambil daripada ensemble unitari, pada keadaan kuantum awal.
2. Ukur hasilnya dalam asas pengiraan, untuk mendapatkan rentetan bit $\ket{b}$.
3. Simpan snapshot klasik, $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, dengan $\mathcal{M}^{-1}$ ialah songsangan saluran pengukuran $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Ulang proses ini sebanyak $T$ kali, dengan koleksi $T$ snapshot $\{\hat{\rho}_t\}$ membentuk bayangan klasik $S(\rho, T)$.
5. Anggarkan observable $o$ bagi pengendali $O$ sebagai $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Varian tomografi bayangan

Pilihan ensemble unitari, iaitu asas, dalam langkah 1 protokol, mempengaruhi kompleksiti sampel, varians anggaran, dan kelengkapan maklumat bagi sesuatu tomografi bayangan. Banyak kumpulan ensemble telah dipertimbangkan dalam literatur, masing-masing dengan pertukaran antara kompleksiti dan kecekapan penganggaran observable. Sifat varian yang paling biasa digunakan diringkaskan dalam jadual berikut, walaupun senarai ini tidak lengkap.

Varian	Ensemble	Kompleksiti sampel	Kedalaman litar	Observable yang sesuai
Pauli tempatan [1]	Pauli qubit tunggal $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	tempatan/beberapa-jasad
Clifford global [1]	Clifford $N$-qubit	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	pangkat rendah, fidelity
Gaussian fermion [2], [3]	litar matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	fermion
bayangan cetek [4], [5]	litar kedalaman rendah	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	bukan tempatan (mesra perkakasan)

Istilah yang digunakan dalam jadual ditakrifkan seperti berikut:

• "kompleksiti sampel" merujuk kepada bilangan pengukuran yang diperlukan untuk sesuatu varian tomografi bayangan.

• "kedalaman litar" menggambarkan betapa mahalnya pengukuran untuk dilaksanakan pada perkakasan kuantum.

• $N$ ialah bilangan qubit.

• $k$ ialah kelokalan observable (sebagai contoh, $\langle Z_i \rangle$ mempunyai $k=1$ dan $\langle Z_i Z_j\rangle$ mempunyai $k=2$).

• $M$ ialah bilangan observable yang hendak dianggar serentak.

• $\varepsilon$ ialah toleransi ralat aditif.

• $\mu$ ialah bilangan fermion.

Walaupun ensemble Pauli tempatan paling mesra perkakasan, ia menanggung penalti kompleksiti sampel $3^k$, yang menjadikannya kurang sesuai untuk penganggaran sifat bukan tempatan. Secara lebih umum, had kompleksiti sampel bayangan tempatan berskala sebagai $4^k$, bertambah baik kepada $3^k$ dalam kes khas observable hasil darab tensor.

Sebaliknya, walaupun ensemble Clifford global boleh digunakan untuk menganggar observable bukan tempatan, kedalaman litarnya tidak menguntungkan. Bayangan Gaussian fermion menawarkan kegunaan dalam simulasi fizik banyak-jasad, kimia, dan bahan, tetapi memerlukan litar matchgate, yang mungkin membawa overhead yang ketara pada peranti berasaskan gate tujuan am.

Bayangan cetek mempunyai kedalaman litar di antara Pauli tempatan dan Clifford global, menjadikannya varian yang sangat menjanjikan untuk pelaksanaan praktikal, walaupun masih terdapat beberapa persoalan terbuka tentang kedalaman litar optimum bagi sesuatu observable.

Apabila sempadan penyelidikan tomografi bayangan terus berkembang, antara masalah yang paling relevan untuk ditangani termasuk memilih ensemble pengukuran secara optimum bagi set observable sasaran, menambah baik had kompleksiti sampel, dan memperluas protokol kepada observable yang berevolusi masa. Walau bagaimanapun, tomografi bayangan masih merupakan salah satu aliran kerja terkemuka untuk mengekstrak maklumat klasik maksimum daripada keadaan kuantum menggunakan bilangan pengukuran minimum, menjadikannya ideal untuk mencirikan keadaan kuantum pada skala di mana tomografi penuh dan simulasi klasik kedua-duanya mustahil.

Ensemble Pauli tempatan

Tutorial ini menggunakan asas Pauli rawak, iaitu gate Clifford qubit tunggal, untuk mengumpul bayangan dan menganggar sifat keadaan GHZ. Perhatikan bahawa pengukuran Pauli ialah pilihan yang sangat baik untuk menganggar observable tempatan (seperti ketumpatan spin, fungsi korelasi tempatan, dan ketumpatan kebarangkalian kedudukan), tetapi tidak untuk yang global (seperti jumlah tenaga, momentum, dan cas).

Di sini, asas Pauli rawak $\{X, Y, Z\}$ bagi setiap qubit dipilih secara bebas dan putaran qubit tunggal yang sepadan digunakan, diikuti dengan pengukuran:

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

Bagi ensemble ini, songsangan saluran pengukuran ialah:

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

Faktor $3$ datang daripada purata ke atas tiga asas Pauli, yang menghasilkan saluran penyahkutuban dengan parameter $\frac{1}{3}$. Penganggar penuh kemudian difaktorkan sebagai hasil darab penganggar qubit tunggal. Seperti yang ditunjukkan dalam jadual, bagi $M$ observable, masing-masing dengan kelokalan $k$, bilangan bayangan Pauli yang diperlukan untuk menganggar semua observable kepada ralat aditif $\varepsilon$ dengan kebarangkalian tinggi ialah $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. Saiz sistem tidak menjejaskan kos pengukuran selagi observable kekal tempatan, jadi observable yang banyak secara eksponen boleh dianggar serentak dengan kos logaritma sahaja.

Median bagi min

Min biasa ke atas $T$ snapshot mengurangkan varians anggaran satu-shot bagi sesuatu observable sebagai $\frac{1}{T}$, tetapi sensitif terhadap outlier. Median bagi min boleh digunakan untuk mengatasi perkara ini. Prosesnya melibatkan perkara berikut:

• Memisahkan $T$ snapshot kepada $K$ kumpulan yang sama saiz
• Mengira min dalam setiap kumpulan
• Mengembalikan median merentasi min kumpulan

Tutorial ini melaksanakan tomografi bayangan Pauli tempatan untuk menganggar sifat keadaan GHZ, mula-mula menggunakan contoh simulasi 10-qubit dan kemudian menunjukkan eksperimen perkakasan 100-qubit.

Keperluan

Sebelum memulakan tutorial ini, pastikan anda telah memasang perkara berikut:

• Qiskit SDK v2.0 atau lebih baharu, dengan sokongan visualisasi

• IBM Quantum Compute v0.28 atau lebih baharu (pip install qiskit-ibm-runtime)

• Qiskit Paulice (pip install qiskit-paulice)

Persediaan

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

Fungsi pembantu berikut digunakan oleh contoh perkakasan berskala besar kemudian dalam tutorial:

• bad_cz, bad_readout, bad_coherence — Menanya data penentukuran backend untuk mengenal pasti gate, qubit, dan tepi yang melebihi ambang ralat atau koheren, supaya ia boleh dikecualikan daripada susun atur litar.

• parallel_ghz — Membina litar GHZ bagi saiz tertentu menggunakan carian pertama-lebar ke atas peta gandingan, menyebarkan lapisan CNOT secara selari untuk meminimumkan kedalaman litar.

• weighted_coverage, active_wires, z_trace_backward — Mengira sejauh mana permukaan ralat litar GHZ dilindungi oleh set semakan pariti ruang-masa tertentu, menggunakan teknik perambatan ke belakang daripada [6].

• add_measurements — Menambah pengukuran qubit GHZ dan qubit semakan ke dalam daftar klasik yang berasingan supaya ia boleh dipilih pasca secara bebas.

• build_shadow_circuits — Menghasilkan satu litar bagi setiap asas Pauli rawak dengan menambah putaran qubit tunggal yang sesuai sebelum pengukuran.

• parse_hw_results_w_checks — Menyahpek output Sampler kepada tiga serangkai (ghz_bits, check_bits, bases), satu bagi setiap shot.

• adapt_to_flat_bits, estimate_zz_nn_hw — Memetakan bit daftar GHZ kembali ke dalam ruang indeks peranti penuh dan mengira anggaran median bagi min $\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

Contoh simulator berskala kecil

Pertama, sediakan keadaan GHZ 10-qubit menggunakan simulasi statevector dan anggarkan sifatnya menggunakan bayangan Pauli tempatan. Keadaan GHZ ialah keadaan kuantum terkait secara maksimum dengan tiga atau lebih qubit yang berfungsi sebagai penanda aras untuk keterkaitan multipartit. Ia ialah blok binaan biasa untuk protokol komunikasi kuantum dan kriptografi. Selain itu, litar GHZ mempunyai sifat analisis yang diketahui sepenuhnya, menjadikannya ideal untuk pengesahan dalam tutorial ini.

Langkah 1: Petakan input klasik kepada masalah kuantum

Sediakan keadaan GHZ 10-qubit yang ideal.

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)

Langkah 2: Optimumkan untuk pelaksanaan pada perkakasan kuantum

Seterusnya, takrifkan fungsi yang mengumpul snapshot bayangan dengan memutar keadaan GHZ ke dalam asas Pauli rawak. Ini mencerminkan apa yang dilakukan oleh primitif Sampler pada perkakasan dalam contoh berskala besar.

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

Langkah 3: Laksanakan menggunakan primitif Qiskit

Untuk simulasi ini, gunakan statevector secara langsung.

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.

Langkah 4: Pasca-proses dan kembalikan hasil dalam format klasik yang dikehendaki

Takrifkan penganggar bayangan dan penganggar median bagi min (MOM) untuk asas pengukuran rawak yang digunakan bagi mengumpul snapshot bayangan. Gunakan fungsi ini untuk menganggar sifat keadaan GHZ, termasuk magnetisasi qubit tunggal, korelator dua qubit, dan RDM dua qubit, daripada pengukuran yang dikumpul. Kemudian bandingkan anggaran tersebut dengan nilai tepat.

_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

Visualkan keputusan

Plot pertama menunjukkan nilai magnetisasi qubit tunggal yang dianggar bagi qubit yang berbeza. Plot kedua menunjukkan taburan nilai korelator dua titik yang dianggar bagi semua pasangan qubit. Plot terakhir menunjukkan magnitud entri matriks ketumpatan terturun dua qubit. Koleksi snapshot yang agak besar digunakan untuk saiz sistem yang diberikan, jadi keputusan hampir sama dengan nilai tepat.

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

Menilai penskalaan shot

Dalam tutorial ini, 5,000 snapshot digunakan untuk menganggar observable. Penskalaan shot bagi tomografi bayangan dinilai dengan memplot ralat min bagi salah satu fungsi yang dianggar terhadap bajet shot yang berbeza. Ralat dalam korelator dua qubit menyusut sebagai $1/\sqrt{T}$ seperti yang dijangkakan bagi penumpuan statistik piawai penganggar min. Kepentingan tomografi bayangan ialah faktor pekali malar di hadapan $1/\sqrt{T}$ tidak bertambah dengan bilangan qubit bagi observable tempatan, seperti yang diperiksa dalam tutorial ini.

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

Contoh perkakasan berskala besar

Simulasi sebelum ini mengesahkan bahawa tomografi bayangan Pauli tempatan memulihkan observable yang betul pada vektor keadaan tanpa hingar. Menjalankan protokol yang sama menggunakan 100 qubit pada perkakasan sebenar memperkenalkan dua cabaran yang tidak wujud dalam simulasi: ralat gate dan bacaan terkumpul merentasi litar GHZ yang dalam, dan mencari susun atur qubit optimum memerlukan pencarian graf gandingan peranti dan bukannya memilih susun atur secara manual.

Untuk menangani cabaran pertama, gunakan pengesanan ralat overhead rendah dengan kod ruang-masa [6]. Sebilangan kecil qubit ancilla bersebelahan dengan pokok GHZ digunakan sebagai semakan pariti. Setiap ancilla dibuat terkait dengan dua jiran GHZ sebelum pengukuran, supaya satu ralat balikan bit pada mana-mana jiran akan membalikkan hasil ancilla daripada 0 kepada 1. Snapshot yang mana-mana ancilla semakan membaca 1 dibuang, meninggalkan subset yang lebih bersih untuk menganggar $\langle Z_i Z_{i+1}\rangle$.

Cabaran kedua ditangani dengan mencari melalui susun atur pokok GHZ calon dan memilih yang memaksimumkan liputan semakan pariti, iaitu pecahan wayar litar yang dipantau oleh sekurang-kurangnya satu semakan, dalam lingkungan peningkatan kedalaman litar yang boleh diterima.

Langkah 1-4 (dilaksanakan dalam satu blok kod)

Aliran kerja empat langkah yang sama kini diaplikasikan pada keadaan GHZ 100 qubit pada perkakasan kuantum sebenar, dengan menambah dua ciri khusus perkakasan yang tidak wujud dalam simulasi:

• Langkah 1 mencirikan hingar backend (ralat gate, ralat bacaan, dan masa T2) supaya qubit dan tepi yang tidak boleh dipercayai dikecualikan daripada susun atur litar.

• Langkah 2 mencari melalui susun atur pokok GHZ yang mungkin untuk mencari satu yang memaksimumkan liputan semakan pariti dalam kedalaman litar yang boleh diterima — menggantikan rantai remeh yang digunakan dalam simulasi.

• Langkah 3 membina litar bayangan (satu bagi setiap asas Pauli rawak), mentranspilkannya kepada gate asli perkakasan, dan menghantarnya kepada primitif Sampler dengan shots=8 (setiap litar menghasilkan lapan rentetan bit bebas).

• Langkah 4 melakukan pascapemilihan snapshot yang semua qubit semakannya membaca 0 (tiada ralat dikesan), kemudian menganggar $\langle Z_i Z_{i+1}\rangle$ dengan dan tanpa pascapemilihan itu untuk menunjukkan peningkatan daripada pengesanan ralat.

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

Tafsirkan keputusan

Kedua-dua plot menunjukkan $\langle Z_i Z_{i+1}\rangle$ yang dianggar daripada tomografi bayangan bagi setiap pasangan qubit berjiran dalam rantai GHZ 100 qubit, sebelum dan selepas pascapemilihan pada semakan pariti.

Bagi keadaan GHZ ideal, semua korelator jiran terdekat bersamaan $1$. Pada perkakasan, hingar merendahkan anggaran mentah dan memperkenalkan serakan antara qubit. Selepas pengesanan ralat, snapshot yang mengalami sekurang-kurangnya satu ralat boleh dikesan dibuang, dan anggaran yang tinggal ditarik lebih dekat kepada $1$ dengan varians yang berkurang, yang konsisten dengan keadaan berkesan yang lebih bersih.

Peningkatan selepas pengesanan ralat konsisten dengan semakan yang mengesan ralat sebenar.

Langkah seterusnya

Terokai sumber berikut untuk mengetahui lebih lanjut:

Cadangan

• Perumusan am maklumat kuantum
• Tutorial pengesanan ralat overhead rendah menggunakan kod ruang-masa
• Dokumentasi API alat tambah Qiskit Paulice

Rujukan

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