التصوير الظلي (Shadow tomography) لحالة Greenberger-Horne-Zeilinger

تقدير الاستخدام: أقل من 10 ثوانٍ على معالج Heron r2 (ملاحظة: هذا مجرد تقدير. قد يختلف وقت التشغيل لديك.)

مخرجات التعلّم

• كيف يبني التصوير الظلي ظلًا كلاسيكيًا لحالة كمية ولماذا يتطلب قياسات أقل أسيًا من التصوير الكامل للحالة

• كيفية تنفيذ التصوير الظلي المحلي بمؤثرات باولي باستخدام Qiskit لتقدير المغنطة أحادية الكيوبت، والمترابطات ثنائية النقاط، ومصفوفات الكثافة المختزلة ثنائية الكيوبت

• كيفية استخدام qiskit-paulice لكشف الأخطاء بحمل إضافي منخفض باستخدام الشفرات الزمكانية في تجربة عتاد واسعة النطاق

المتطلبات الأساسية

• القياسات

• التصوير المقطعي للحالة الكمية (Quantum state tomography)

• كشف الأخطاء باستخدام الشفرات الزمكانية

الخلفية

يتطلب التصوير الكامل للحالة الكمية لنظام من $N$ كيوبت عددًا أسيًا من القياسات لتقدير جميع المعاملات المستقلة البالغ عددها $4^N -1$ لمصفوفة الكثافة. فعلى سبيل المثال، لنظام من 50 كيوبت، تتطلب هذه العملية ما يقارب $10^{30}$ لقطة، وهو أكثر من مليار ضعف عدد النجوم في الكون المرصود. التصوير الظلي [1] هو بروتوكول في معلومات الكم يبني "ظلًا كلاسيكيًا" للحالة الكمية، يمكن منه تقدير المرصودات المستهدفة باستخدام عدد من اللقطات يتناسب لوغاريتميًا مع عدد المرصودات ولا يعتمد على حجم النظام $N$، مما يجعله ذا صلة بمجالات كثيرة في الحوسبة الكمية، منها تقدير الطاقة، وحسابات إنتروبيا التشابك ودوال الارتباط، وتقدير النواة في تعلّم الآلة الكمي، وغيرها.

يوضح هذا الدرس كيفية استخدام التصوير الظلي لتقدير المرصودات المحلية، مثل المغنطة أحادية الكيوبت $\langle Z_i \rangle$ للكيوبت $i$، ودالة الارتباط ثنائية النقاط $\langle Z_i Z_j\rangle$، ومصفوفة الكثافة المختزلة ثنائية الكيوبت (2-RDM)، لحالة Greenberger-Horne-Zeilinger (GHZ) عالية التشابك. يوضح الدرس أولًا سير عمل تصوير ظلي صغير النطاق على حالة GHZ من 10 كيوبت باستخدام محاكاة متجه الحالة. ثم يعرض تجربة بحجم 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. تقدير المرصود $o$ للمؤثر $O$ على أنه $\frac{1}{T}\sum_{t=1}^{T}\text{Tr}(O\hat{\rho}_t)$.

متغيرات التصوير الظلي

يؤثر اختيار المجموعة الوحدوية، أي الأساس، في الخطوة 1 من البروتوكول على تعقيد العينات وتباين التقديرات واكتمال المعلومات في تصوير ظلي معيّن. وقد دُرست مجموعات كثيرة في الأدبيات، لكل منها مقايضات بين التعقيد وكفاءة تقدير المرصودات. تُلخَّص خصائص المتغيرات الأكثر استخدامًا في الجدول التالي، رغم أن هذه القائمة ليست شاملة.

المتغير	المجموعة	تعقيد العينات	عمق الدائرة	المرصودات المناسبة
باولي المحلي [1]	مؤثرات باولي أحادية الكيوبت $\{X, Y, Z\}$	$\mathcal{O}(3^k\text{log}M/\varepsilon^2)$	1	محلية/قليلة الأجسام
Clifford العام [1]	Clifford على $N$ كيوبت	$\mathcal{O}(\text{Tr}(O^2)\text{log}M/\varepsilon^2)$	$\mathcal{O}(N\text{log}N)$	منخفضة الرتبة، الدقة (fidelity)
غاوسي فرميوني [2], [3]	دوائر matchgate	$\mathcal{O}(\binom{\mu}{k}\frac{1}{\varepsilon^2})$	$\mathcal{O}(N^2)$	فرميونية
الظلال الضحلة [4], [5]	دوائر منخفضة العمق	$\mathcal{O}(k2^k\text{log}M/\varepsilon^2)$	$\mathcal{O}(\text{log}N)$	غير محلية (ملائمة للعتاد)

تُعرَّف المصطلحات المستخدمة في الجدول على النحو التالي:

• "تعقيد العينات" يشير إلى عدد القياسات اللازمة لمتغير تصوير ظلي معيّن.

• "عمق الدائرة" يصف مدى كلفة تنفيذ القياسات على العتاد الكمي.

• $N$ هو عدد الكيوبتات.

• $k$ هو محلية المرصود (على سبيل المثال، $\langle Z_i \rangle$ له $k=1$ و$\langle Z_i Z_j\rangle$ له $k=2$).

• $M$ هو عدد المرصودات المراد تقديرها في وقت واحد.

• $\varepsilon$ هو تسامح الخطأ الجمعي.

• $\mu$ هو عدد الفرميونات.

مع أن مجموعة باولي المحلية هي الأكثر ملاءمة للعتاد، فإنها تدفع عقوبة تعقيد عينات قدرها $3^k$، مما يجعلها غير مناسبة لتقدير الخصائص غير المحلية. وبشكل أعم، يتناسب حد تعقيد العينات للظلال المحلية مع $4^k$، ويتحسن إلى $3^k$ في الحالة الخاصة للمرصودات ذات الجداء التنسوري.

في المقابل، بينما يمكن استخدام مجموعة Clifford العامة لتقدير المرصودات غير المحلية، فإن عمق دائرتها غير مناسب. تقدم الظلال الغاوسية الفرميونية فائدة في محاكاة فيزياء الأجسام المتعددة والكيمياء والمواد، لكنها تتطلب دوائر matchgate، قد تحمل حملًا إضافيًا كبيرًا على الأجهزة العامة القائمة على البوابات.

تتميز الظلال الضحلة بعمق دائرة بين باولي المحلي وClifford العام، مما يجعلها متغيرًا واعدًا جدًا للتنفيذ العملي، رغم وجود بعض الأسئلة المفتوحة بشأن العمق الأمثل للدائرة لمرصود معيّن.

مع استمرار توسع حدود أبحاث التصوير الظلي، تشمل بعض المشكلات الأكثر صلة التي ينبغي معالجتها اختيار مجموعة القياس الأمثل لمجموعة مرصودات مستهدفة، وتحسين حدود تعقيد العينات، وتوسيع البروتوكول ليشمل المرصودات المتطورة زمنيًا. ومع ذلك، لا يزال التصوير الظلي أحد أبرز سير العمل لاستخراج أقصى قدر من المعلومات الكلاسيكية من حالة كمية باستخدام أقل عدد من القياسات، مما يجعله مثاليًا لتوصيف الحالات الكمية على نطاق يستحيل فيه كل من التصوير الكامل والمحاكاة الكلاسيكية.

مجموعة باولي المحلية

يستخدم هذا الدرس أساس باولي العشوائي، أي بوابات Clifford أحادية الكيوبت، لجمع الظلال وتقدير خصائص حالة GHZ. لاحظ أن قياسات باولي خيار ممتاز لتقدير المرصودات المحلية (مثل كثافة اللف، ودوال الارتباط المحلية، وكثافة احتمال الموضع)، لكنها ليست كذلك للمرصودات العامة (مثل الطاقة الكلية، والزخم، والشحنة).

هنا، يُختار أساس باولي عشوائي $\{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$، لكل منها محلية $k$، بخطأ جمعي $\varepsilon$ باحتمال عالٍ هو $\mathcal{O}(\frac{3^k\text{log}(M)}{\varepsilon^2})$. لا يؤثر حجم النظام في تكلفة القياس ما دامت المرصودات محلية، لذا يمكن تقدير عدد أسي من المرصودات في وقت واحد بتكلفة لوغاريتمية فقط.

الوسيط للمتوسطات (Median of means)

يقلل المتوسط البسيط على $T$ لقطة تباين التقدير أحادي اللقطة لمرصود بمقدار $\frac{1}{T}$، لكنه حساس للقيم الشاذة. يمكن استخدام الوسيط للمتوسطات لعلاج ذلك. وتتضمن العملية ما يلي:

• تقسيم اللقطات البالغ عددها $T$ إلى $K$ مجموعة متساوية الحجم
• حساب المتوسط داخل كل مجموعة
• إرجاع الوسيط عبر متوسطات المجموعات

يجري هذا الدرس تصويرًا ظلليًا بمؤثرات باولي المحلية لتقدير خصائص حالة GHZ، أولًا باستخدام مثال محاكى من 10 كيوبت ثم بعرض تجربة عتاد من 100 كيوبت.

المتطلبات

قبل البدء بهذا الدرس، تأكد من تثبيت ما يلي:

• Qiskit SDK الإصدار 2.0 أو أحدث، مع دعم التصوير المرئي

• IBM Quantum Compute الإصدار 0.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 لتحديد البوابات والكيوبتات والحواف التي تتجاوز عتبات الخطأ أو الترابط، بحيث يمكن استبعادها من تخطيط الدائرة.

• parallel_ghz — تبني دائرة GHZ بحجم معيّن باستخدام بحث بالعرض أولًا (breadth-first) على خريطة الاقتران، موزعةً طبقات CNOT بالتوازي لتقليل عمق الدائرة.

• weighted_coverage، active_wires، z_trace_backward — تحسب مقدار حماية سطح أخطاء دائرة GHZ بمجموعة معيّنة من فحوص التكافؤ الزمكانية، باستخدام تقنية الانتشار العكسي من [6].

• add_measurements — تُلحق قياسات كيوبتات GHZ وكيوبتات الفحص في سجلات كلاسيكية منفصلة بحيث يمكن انتقاؤها لاحقًا بشكل مستقل.

• build_shadow_circuits — تنتج دائرة واحدة لكل أساس باولي عشوائي بإلحاق الدورانات أحادية الكيوبت المناسبة قبل القياس.

• 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

مثال المحاكي صغير النطاق

أولًا، حضّر حالة GHZ من 10 كيوبت باستخدام محاكاة متجه الحالة وقدّر خصائصها باستخدام ظلال باولي المحلية. حالة GHZ هي حالة كمية ذات تشابك أقصى مع ثلاثة كيوبتات أو أكثر، وتعمل كمعيار للتشابك متعدد الأطراف. وهي لبنة بناء شائعة لبروتوكولات الاتصال الكمي والتشفير. وإضافةً إلى ذلك، فإن لدائرة GHZ خصائص تحليلية معروفة بالكامل، مما يجعلها مثالية للتحقق في هذا الدرس.

الخطوة 1: ربط المدخلات الكلاسيكية بمسألة كمية

حضّر حالة GHZ مثالية من 10 كيوبت.

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 primitive على العتاد في المثال واسع النطاق.

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 primitives

في هذه المحاكاة، استخدم متجه الحالة مباشرةً.

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 لقطة لتقدير القيم المرصودة. ونقيّم تدرّج عدد اللقطات في التصوير الظلي برسم متوسط الخطأ لإحدى الدوال المقدّرة مقابل ميزانيات لقطات مختلفة. ينخفض الخطأ في مرتبط الكيوبتين بنسبة $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

مثال على عتاد واسع النطاق

تؤكد المحاكاة السابقة أن التصوير الظلي لباولي المحلي يستعيد المرصودات الصحيحة على متجه حالة خالٍ من الضوضاء. أما تشغيل البروتوكول نفسه باستخدام 100 كيوبت على عتاد حقيقي فيُدخل تحديين لا يظهران في المحاكاة: تتراكم أخطاء البوابات والقراءة عبر دائرة GHZ العميقة، ويتطلب إيجاد أفضل تخطيط للكيوبتات البحث في الرسم البياني للربط في الجهاز بدلًا من اختيار التخطيط يدويًا.

لمعالجة التحدي الأول، استخدم كشف الأخطاء بحمل إضافي منخفض باستخدام الشفرات الزمكانية [6]. نستخدم عددًا صغيرًا من كيوبتات ancilla المجاورة لشجرة GHZ كفحوصات تكافؤ. يتشابك كل ancilla مع اثنين من جيرانه في GHZ قبل القياس، بحيث يؤدي خطأ انقلاب بت واحد في أيٍّ من الجارين إلى قلب نتيجة الـ ancilla من 0 إلى 1. تُهمل اللقطات التي يقرأ فيها أي ancilla فحص القيمة 1، فتبقى مجموعة أنظف نقدّر منها $\langle Z_i Z_{i+1}\rangle$.

أما التحدي الثاني فنعالجه بالبحث بين تخطيطات شجرة GHZ المرشحة واختيار التخطيط الذي يعظّم تغطية فحوصات التكافؤ، أي نسبة أسلاك الدارة التي يراقبها فحص واحد على الأقل، ضمن زيادة مقبولة في عمق الدارة.

الخطوات 1-4 (تُنفَّذ في كتلة شيفرة واحدة)

نطبّق الآن سير العمل نفسه ذي الخطوات الأربع على حالة GHZ بـ 100 كيوبت على عتاد كمي حقيقي، مع إضافة ميزتين خاصتين بالعتاد لا تظهران في المحاكاة:

• الخطوة 1 تحدد خصائص ضوضاء الـ Backend (أخطاء البوابات، وأخطاء القراءة، وأزمنة T2) لاستبعاد الكيوبتات والحواف غير الموثوقة من تخطيط الدارة.

• الخطوة 2 تبحث بين تخطيطات شجرة GHZ الممكنة لإيجاد التخطيط الذي يعظّم تغطية فحوصات التكافؤ ضمن عمق دائرة مقبول، بدلًا من السلسلة البسيطة المستخدمة في المحاكاة.

• الخطوة 3 تبني دوائر الظل (واحدة لكل أساس باولي عشوائي)، وتحوّلها بالـ transpile إلى بوابات أصلية للعتاد، وترسلها إلى Sampler مع shots=8 (تنتج كل دائرة ثماني سلاسل بتات مستقلة).

• الخطوة 4 تُجري postselection للقطات التي تقرأ فيها كل كيوبتات الفحص القيمة 0 (أي لم يُكتشف خطأ)، ثم تقدّر $\langle Z_i Z_{i+1}\rangle$ مع هذا الـ postselection وبدونه لإظهار التحسن الناتج عن كشف الأخطاء.

# ── 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$ المقدّرة من توموغرافيا الظل لكل زوج كيوبتات متجاور في سلسلة GHZ ذات 100 كيوبت، قبل إجراء postselection على فحوصات التكافؤ وبعده.

في حالة 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