Тіньова томографія стану Грінберґера—Горна—Цайлінґера

Оцінка використання: менше 10 секунд на процесорі Heron r2 (ПРИМІТКА: це лише оцінка. Ваш час виконання може відрізнятися.)

Результати навчання

• Як тіньова томографія будує класичну тінь квантового стану та чому вона вимагає експоненціально менше вимірювань, ніж повна томографія стану

• Як реалізувати локальну паулі-тіньову томографію за допомогою Qiskit для оцінювання одно-кубітових намагнеченостей, двоточкових кореляторів і двокубітових зведених матриць густини

• Як використовувати qiskit-paulice для виявлення помилок з низькими накладними витратами за допомогою просторово-часових кодів у масштабному апаратному експерименті

Передумови

• Вимірювання

• Томографія квантового стану

• Виявлення помилок за допомогою просторово-часових кодів

Передісторія

Повна томографія квантового стану системи з $N$ кубітів вимагає експоненціально багато вимірювань, щоб оцінити всі $4^N -1$ незалежних параметрів матриці густини. Наприклад, для системи з 50 кубітів цей процес вимагає порядку $10^{30}$ шотів, що більш ніж у мільярд разів перевищує кількість зірок у видимому Всесвіті. Тіньова томографія [1] — це протокол квантової інформації, який будує «класичну тінь» квантового стану, з якої цільові observable можна оцінювати, використовуючи кількість шотів, що масштабується логарифмічно від кількості observable і не залежить від розміру системи $N$. Це робить її значущою в багатьох галузях квантових обчислень, зокрема в оцінюванні енергії, обчисленнях ентропії заплутаності та кореляційних функцій, оцінюванні ядер у квантовому машинному навчанні тощо.

Цей навчальний посібник показує, як використовувати тіньову томографію для оцінювання локальних observable, як-от одно-кубітова намагніченість $\langle Z_i \rangle$ для кубіта $i$, двоточкова кореляційна функція $\langle Z_i Z_j\rangle$ та двокубітова зведена матриця густини (2-RDM) сильно заплутаного стану Грінберґера—Горна—Цайлінґера (GHZ). Спочатку посібник ілюструє невеликий робочий процес тіньової томографії на 10-кубітному стані GHZ за допомогою симуляції вектора стану. Потім він показує демонстрацію зі 100 кубітами на реальному квантовому обладнанні.

Протокол тіньової томографії

Протокол тіньової томографії має такі кроки:

1. Застосуйте до початкового квантового стану випадковий унітарний оператор $U$, взятий з унітарного ансамблю.
2. Виміряйте результат у обчислювальному базисі, отримавши бітовий рядок $\ket{b}$.
3. Збережіть класичний знімок $\hat{\rho} = \mathcal{M}^{-1}(U^{\dagger}\ket{b}\bra{b}U)$, де $\mathcal{M}^{-1}$ — обернення каналу вимірювання $\mathcal{M}(\rho)=\mathbb{E}[U^{\dagger}\ket{b}\bra{b}U]$.
4. Повторіть процес $T$ разів, де сукупність $T$ знімків $\{\hat{\rho}_t\}$ утворює класичну тінь $S(\rho, T)$.
5. Оцініть observable $o$ оператора $O$ як $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

Варіанти тіньової томографії

Вибір унітарного ансамблю, тобто базису, на кроці 1 протоколу впливає на вибіркову складність, дисперсію оцінок і повноту інформації певної тіньової томографії. У літературі розглядалося багато ансамблевих груп, кожна з компромісами між складністю та ефективністю оцінювання observable. Властивості найпоширеніших варіантів узагальнено в таблиці нижче, хоча цей перелік не вичерпний.

Варіант	Ансамбль	Вибіркова складність	Глибина circuit	Придатні observable
локальний Паулі [1]	однокубітові Паулі $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	локальні/кількачастинкові
глобальний Кліффорд [1]	$N$-кубітний Кліффорд	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	малоранґові, вірність
фермionний гауссів [2], [3]	circuit matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	ферміонні
мілкі тіні [4], [5]	circuit малої глибини	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	нелокальні (зручні для обладнання)

Терміни, використані в таблиці, визначені так:

• «вибіркова складність» означає кількість вимірювань, потрібних для певного варіанта тіньової томографії.

• «глибина circuit» описує, наскільки дорого реалізувати вимірювання на квантовому обладнанні.

• $N$ — кількість кубітів.

• $k$ — локальність observable (наприклад, $\langle Z_i \rangle$ має $k=1$, а $\langle Z_i Z_j\rangle$ має $k=2$).

• $M$ — кількість observable, які оцінюються одночасно.

• $\varepsilon$ — допустима адитивна похибка.

• $\mu$ — кількість ферміонів.

Хоча локальний паулівський ансамбль є найзручнішим для обладнання, він платить штраф вибіркової складності $3^k$, що робить його непридатним для оцінювання нелокальних властивостей. Загалом межа вибіркової складності для локальних тіней масштабується як $4^k$, покращуючись до $3^k$ у особливому випадку observable у вигляді тензорних добутків.

Натомість, хоча глобальний кліффордів ансамбль можна використовувати для оцінювання нелокальних observable, його глибина circuit несприятлива. Ферміонні гауссові тіні корисні в симуляціях фізики багатьох тіл, хімії та матеріалів, але вимагають circuit matchgate, які можуть мати значні накладні витрати на універсальних пристроях на основі гейтів.

Мілкі тіні мають глибину circuit між локальним паулівським і глобальним кліффордовим ансамблями, що робить їх дуже перспективним варіантом для практичної реалізації, хоча досі є відкриті питання щодо оптимальної глибини circuit для заданого observable.

Оскільки фронтир досліджень тіньової томографії продовжує розширюватися, серед найактуальніших задач — оптимальний вибір ансамблю вимірювань для цільового набору observable, покращення меж вибіркової складності та розширення протоколу на observable, що еволюціонують у часі. Проте тіньова томографія залишається одним із провідних робочих процесів для добування максимуму класичної інформації з квантового стану за мінімальної кількості вимірювань, що робить її ідеальною для характеризації квантових станів у масштабі, де неможливі і повна томографія, і класична симуляція.

Локальний паулівський ансамбль

Цей навчальний посібник використовує випадковий паулівський базис, тобто однокубітові гейти Кліффорда, щоб збирати тіні й оцінювати властивості стану GHZ. Зауважте, що вимірювання Паулі — чудовий вибір для оцінювання локальних observable (як-от густина спіну, локальні кореляційні функції та густина ймовірності положення), але не глобальних (як-от повна енергія, імпульс і заряд).

Тут для кожного кубіта незалежно обирається випадковий базис Паулі $\{X, Y, Z\}$ і застосовується відповідне однокубітове обертання, після чого виконується вимірювання:

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

Для цього ансамблю обернення каналу вимірювання таке:

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

Множник $3$ виникає від усереднення за трьома базисами Паулі, що дає деполяризаційний канал із параметром $\frac{1}{3}$. Повний оцінювач тоді факторизується як добуток однокубітових оцінювачів. Як показано в таблиці, для $M$ observable, кожне з локальністю $k$, кількість паулівських тіней, потрібних для оцінювання всіх observable з адитивною похибкою $\varepsilon$ з високою ймовірністю, дорівнює $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. Розмір системи не впливає на вартість вимірювань, доки observable залишаються локальними, тож експоненціально багато observable можна оцінювати одночасно лише з логарифмічною вартістю.

Медіана середніх

Просте середнє за $T$ знімками зменшує дисперсію одношотової оцінки observable як $\frac{1}{T}$, але чутливе до викидів. Щоб це виправити, можна використати медіану середніх. Процес включає таке:

• Розбиття $T$ знімків на $K$ груп однакового розміру
• Обчислення середнього в кожній групі
• Повернення медіани середніх за групами

Цей навчальний посібник виконує локальну паулівську тіньову томографію, щоб оцінити властивості стану GHZ: спочатку на 10-кубітному змодельованому прикладі, а потім демонструє 100-кубітний експеримент на обладнанні.

Вимоги

Перш ніж почати цей навчальний посібник, переконайтеся, що у вас встановлено таке:

• Qiskit SDK v2.0 або новіший з підтримкою візуалізації

• IBM Quantum Compute v0.28 або новіший (pip install qiskit-ibm-runtime)

• Qiskit Paulice (pip install qiskit-paulice)

Налаштування

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

Наведені нижче допоміжні функції використовуються в масштабному апаратному прикладі далі в посібнику:

• bad_cz, bad_readout, bad_coherence — запитують дані калібрування backend, щоб визначити гейти, кубіти та ребра, які перевищують пороги помилок або когерентності, щоб їх можна було виключити з layout circuit.

• parallel_ghz — будує circuit GHZ заданого розміру за допомогою пошуку вшир по coupling map, паралельно розподіляючи шари CNOT, щоб мінімізувати глибину circuit.

• weighted_coverage, active_wires, z_trace_backward — обчислюють, яка частина поверхні помилок circuit GHZ захищена заданим набором просторово-часових перевірок парності, використовуючи техніку зворотного поширення з [6].

• add_measurements — додає вимірювання кубітів GHZ і кубітів перевірки в окремі класичні регістри, щоб їх можна було постселектувати незалежно.

• build_shadow_circuits — створює по одному circuit для кожного випадкового базису Паулі, додаючи відповідні однокубітові обертання перед вимірюванням.

• parse_hw_results_w_checks — розпаковує вивід Sampler у трійки (ghz_bits, check_bits, bases), по одній на шот.

• adapt_to_flat_bits, estimate_zz_nn_hw — відображають біти регістра GHZ назад у повний індексний простір пристрою та обчислюють оцінку $\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

Невеликий приклад із симулятором

Спочатку підготуйте 10-кубітний стан GHZ за допомогою симуляції вектора стану й оцініть його властивості за допомогою локальних паулівських тіней. Стан GHZ — це максимально заплутаний квантовий стан із трьома чи більше кубітами, що слугує еталоном багаточастинкової заплутаності. Він є поширеним будівельним блоком протоколів квантового зв'язку та криптографії. Крім того, circuit GHZ має повністю відомі аналітичні властивості, що робить його ідеальним для перевірки в цьому посібнику.

Крок 1: Відобразіть класичні вхідні дані на квантову задачу

Підготуйте ідеальний 10-кубітний стан 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)

Крок 2: Оптимізуйте для виконання на квантовому обладнанні

Далі визначте функцію, яка збирає тіньові знімки, обертаючи стан GHZ у випадкові базиси Паулі. Це віддзеркалює те, що примітив Sampler робить на обладнанні в масштабному прикладі.

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

Крок 3: Виконайте за допомогою примітивів Qiskit

Для цієї симуляції використовуй вектор стану безпосередньо.

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.

Крок 4: Виконайте постобробку та поверніть результат у потрібному класичному форматі

Визнач оцінювач тіней (shadow estimator) та оцінювач медіани середніх (MOM) для випадкових базисів вимірювання, використаних для збору знімків тіней. Використай ці функції, щоб оцінити властивості стану GHZ, зокрема намагніченість окремого кубіта, двокубітні кореляції та двокубітну RDM, за зібраними вимірюваннями. Потім порівняй оцінки з точними значеннями.

_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

Візуалізація результатів

Перший графік показує оцінені значення намагніченості окремих кубітів для різних кубітів. Другий графік показує розподіл оцінених значень двоточкових кореляторів для всіх пар кубітів. Останній графік показує величину елементів двокубітної редукованої матриці густини. Для даного розміру системи використано досить велику кількість знімків, тому результати близькі до точних.

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

Оцінка масштабування за кількістю shots

У цьому посібнику для оцінки спостережуваних використовується 5 000 знімків. Масштабування томографії тіней за кількістю shots оцінюється шляхом побудови графіка середньої похибки однієї з оцінюваних функцій залежно від різних бюджетів shots. Похибка двокубітного корелятора спадає як $1/\sqrt{T}$, як і очікується для стандартної статистичної збіжності оцінювача середнього. Важливість томографії тіней полягає в тому, що стала передмножника перед $1/\sqrt{T}$ не зростає зі збільшенням кількості кубітів для локальних спостережуваних, таких як ті, що розглядаються в цьому посібнику.

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

Приклад великомасштабного запуску на обладнанні

Попередня симуляція підтверджує, що локальна томографія тіней Паулі відновлює правильні observable на безшумному векторі стану. Запуск того самого протоколу зі 100 кубітами на реальному обладнанні вносить дві проблеми, які не виникають у симуляції: похибки гейтів і зчитування накопичуються вздовж глибокої схеми GHZ, а пошук оптимального розташування кубітів вимагає перебору графа зв'язності пристрою замість підбору розташування вручну.

Щоб розв'язати першу проблему, використай виявлення помилок з малими накладними витратами за допомогою просторово-часових кодів [6]. Невелика кількість допоміжних кубітів (ancilla), суміжних із деревом GHZ, використовується як перевірки парності. Кожен ancilla заплутується з двома своїми сусідами в GHZ перед вимірюванням, тож одна помилка переворота біта на будь-якому із сусідів змінює результат ancilla з 0 на 1. Знімки, де будь-який перевірочний ancilla показує 1, відкидаються, залишаючи чистішу підмножину, з якої оцінюється $\langle Z_i Z_{i+1}\rangle$.

Другу проблему розв'язує пошук серед кандидатів на розташування дерева GHZ і вибір того, що максимізує покриття перевірками парності, тобто частку ліній схеми, які контролюються принаймні однією перевіркою, за прийнятного збільшення глибини схеми.

Кроки 1-4 (виконуються в одному блоці коду)

Тепер той самий чотирикроковий робочий процес застосовується до 100-кубітного стану GHZ на реальному квантовому обладнанні, з додаванням двох апаратно-специфічних можливостей, які не виникають у симуляції:

• Крок 1 характеризує шум backend (похибки гейтів, похибки зчитування та часи T2), щоб ненадійні кубіти та ребра були виключені з розташування схеми.

• Крок 2 перебирає можливі розташування дерева GHZ, щоб знайти те, яке максимізує покриття перевірками парності за прийнятної глибини схеми, замінюючи тривіальний ланцюжок, використаний у симуляції.

• Крок 3 будує схеми тіней (по одній для кожного випадкового базису Паулі), транспілює їх у нативні для обладнання гейти та надсилає їх до примітива Sampler із shots=8 (кожна схема дає вісім незалежних бітових рядків).

• Крок 4 виконує постселекцію знімків, у яких усі перевірочні кубіти показали 0 (помилок не виявлено), а потім оцінює $\langle Z_i Z_{i+1}\rangle$ з постселекцією та без неї, щоб продемонструвати покращення від виявлення помилок.

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

Інтерпретація результатів

Два графіки показують $\langle Z_i Z_{i+1}\rangle$, оцінене за допомогою томографії тіней для кожної пари сусідніх кубітів у 100-кубітному ланцюжку GHZ, до та після постселекції за перевірками парності.

Для ідеального стану GHZ усі кореляторы найближчих сусідів дорівнюють $1$. На обладнанні шум знижує необроблені оцінки та вносить розкид від кубіта до кубіта. Після виявлення помилок знімки, у яких виникла принаймні одна виявна помилка, видаляються, а решта оцінок наближається до $1$ зі зменшеною дисперсією, що узгоджується з чистішим ефективним станом.

Покращення після виявлення помилок узгоджується з тим, що перевірки виявляють справжні помилки.

Наступні кроки

Ознайомся з наступними ресурсами, щоб дізнатися більше:

Рекомендації

• Загальне формулювання квантової інформації
• Посібник з виявлення помилок з малими накладними витратами за допомогою просторово-часових кодів
• Документація API додатка Qiskit Paulice

Посилання

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