PEC s TREX

Probabilistic error cancellation (PEC) je technika pro zmírnění kvantového šumu při výpočtech středních hodnot pozorovatelných. Umožňuje vypočítat nezkreslený odhad středních hodnot. Technika funguje tak, že spouští soubor zašuměných obvodů, do nichž je vložena inverze naučeného modelu šumu. Model šumu je poté aproximován Pauliho kanálem pomocí Pauli twirlingu (známého také jako randomized compiling). Takové modely šumu lze naučit efektivně a lze je mitigovat pomocí PEC (více podrobností viz Probabilistic error cancellation with sparse Pauli-Lindblad models on noisy quantum processors).

Twirled Readout Error eXtinction (TREX) je technika, která zmírňuje chyby měření a je užitečná pro odhad středních hodnot Pauliho pozorovatelných. Spoléhá na twirling měření (náhodné vkládání hradla X před měření a převrácení naměřeného výsledku), který diagonalizuje přenosovou matici chyb odečtu a usnadňuje její inverzi. Přeškálovací člen se učí z diagonálního šumového kanálu měřením twirlovaných kalibračních obvodů, které byly inicializovány ve stavu nula.

Tento návod ukazuje, jak použít PEC ke zmírnění šumu hradel a zkombinovat jej s TREX, aby se při odhadu středních hodnot pozorovatelných zmírnil i šum měření.

Postup je následující:

1. Zadej obvod a pozorovatelnou.
2. Vypočítej přesné střední hodnoty pomocí statevector simulátoru Qiskitu.
3. Vypočítej zašuměné, nemitigované střední hodnoty.
4. Nauč se šum provázávajících vrstev pomocí qiskit-noise-learning.
5. Připrav program PEC + TREX a přidej do něj kalibrační obvod TREX.
6. Porovnej mitigované střední hodnoty.

Zadej obvod a pozorovatelnou

Připrav zrcadlový obvod Trotterizovaného časového vývoje lineárního Isingova modelu v příčném poli na osmi qubitech. Obvod zrcadli, aby se jeho unitární operace vyrušila a vznikl identický obvod. Díky tomu snadno vypočítáš ideální střední hodnotu každé pozorovatelné. Zrcadlení navíc zdvojnásobí hloubku obvodu, což zase zdvojnásobí vliv šumu.

Použij úhel $\theta= \frac{\pi}{2}$, aby všechna hradla v obvodu patřila do Cliffordovy grupy, což umožní rychlejší simulaci pomocí simulátoru stabilizer.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "qiskit_mitigation": "qiskit-mitigation", "qiskit_noise_learning": "qiskit-noise-learning", "samplomatic": "samplomatic"}
_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 warnings

import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

warnings.filterwarnings("ignore")

def construct_ising_circuit(
    num_qubits: int,
    num_trotter_steps: int,
    rx_angle: float,
    barrier: bool = True,
) -> QuantumCircuit:
    circuit = QuantumCircuit(num_qubits)

    for _step in range(num_trotter_steps):
        circuit.rx(rx_angle, range(num_qubits))
        if barrier:
            circuit.barrier()
        for first_qubit in (1, 2):
            for idx in range(first_qubit, num_qubits, 2):
                # equivalent to Rzz(-pi/2):
                circuit.sdg([idx - 1, idx])
                circuit.cz(idx - 1, idx)
        if barrier:
            circuit.barrier()

    return circuit

def mirror_circuit(circuit: QuantumCircuit) -> QuantumCircuit:
    mirror_circ = circuit.copy_empty_like()
    mirror_circ.compose(circuit, inplace=True)
    mirror_circ.barrier()
    mirror_circ.compose(circuit.inverse(), inplace=True)
    mirror_circ.measure_active()
    return mirror_circ

# define observables
num_qubits = 8
observables = [
    SparsePauliOp("I" * ((num_qubits - 1) - q) + "Z" + "I" * q)
    for q in range(num_qubits)
]

# Instantiate circuit
num_trotter_steps = 5
rx_angle = np.pi / 2

circuit = construct_ising_circuit(
    num_qubits, num_trotter_steps, rx_angle, barrier=False
)
mirrored_circuit = mirror_circuit(circuit)
print("Visualization of 1 mirrored Trotter step")
mirror_circuit(
    construct_ising_circuit(num_qubits, 1, rx_angle, barrier=False)
).draw("mpl", scale=0.7)

Visualization of 1 mirrored Trotter step

Vypočítej zašuměné, nemitigované střední hodnoty

QuantumProgram obsahuje kvantový obvod a samplex ve formátu vhodném pro primitivu Executor. MitigationTask.prepare sestaví QuantumProgram z obvodu, pozorovatelných a hodnot parametrů a obvod zabalí do boxů podle zadaných voleb boxování. Nastav volbu boxování measure_annotations na "change_basis", aby se měření otočila do požadovaných bází, a přidej add_tags: "unique_box", aby simulátor mohl namapovat model šumu na boxy obvodu. Program je ve výchozím stavu Pauli twirlován, což na hardwaru mění koherentní šum na Pauliho šum, který mitigační metody předpokládají. num_randomizations určuje, kolik twirlovaných instancí každého obvodu se vzorkuje, a shots_per_randomization určuje, kolik shotů každá instance dostane, takže jejich součin je limit shotů na obvod.

Pokud máš přístup k QPU, můžeš program spustit pomocí Executor vytvořeného z backendu QiskitRuntimeService. Tento příklad používá AerSimulator ke spuštění programu lokálně. Jeho simulátorová volba layer_noise_model vkládá náhodný Pauli-Lindbladův model šumu do každé vrstvy provázávajících hradel, vygenerovaný pomocí random_pauli_lindblad_noise z qiskit-mitigation s chybou na vrstvené hradlo (EPLG) 0.003, což je typické pro současné QPU Heron.

import matplotlib.pyplot as plt
from qiskit.circuit import CircuitInstruction
from qiskit.quantum_info import PauliLindbladMap, pauli_basis
from qiskit.transpiler import PassManager
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import Executor
from qiskit_mitigation import MitigationTask, find_combined_unique_layers
from qiskit_noise_learning.aer_executor.insert_noise_pass import (
    InsertNoisePass,
)
from samplomatic import Tag
from samplomatic.utils import get_annotation

def generate_random_pauli_lindblad_noise(
    layers,
    num_qubits: int,
    eplg: float = 3e-3,
    readout_noise_scale: float = 3e-2,
    seed: int | None = None,
) -> tuple[
    list[tuple[CircuitInstruction, PauliLindbladMap]], PauliLindbladMap
]:
    """Generate random Pauli-Lindblad noise over the full Pauli basis."""
    layer_noise_model = []
    for layer in layers:
        layer_instructions = layer.operation.params[0]
        contain_measurement = any(
            inst.name.startswith("measure") for inst in layer_instructions
        )
        if not contain_measurement:
            edges = [
                (inst.qubits[0]._index, inst.qubits[1]._index)
                for inst in layer_instructions
            ]  # extract the edges of the layer

            basis_paulis = [p for p in pauli_basis(2) if np.sum(p.x + p.z)]
            rng = np.random.default_rng(seed=seed)
            num_generators = len(basis_paulis)
            rates = rng.random(len(basis_paulis)) * eplg / num_generators

            layer_noise = PauliLindbladMap.from_sparse_list(
                [
                    (pauli.to_label(), edge, rate)
                    for pauli, rate in zip(basis_paulis, rates, strict=True)
                    for edge in edges
                ],
                num_qubits=num_qubits,
            )
            layer_noise_model.append((layer, layer_noise))
    # generate a readout noise model
    rng = np.random.default_rng(seed=seed)
    rates = rng.random(num_qubits) * readout_noise_scale
    readout_noise_model = PauliLindbladMap.from_sparse_list(
        [("X", [q], rate) for q, rate in enumerate(rates)],
        num_qubits=num_qubits,
    )
    return layer_noise_model, readout_noise_model

def add_readout_noise(program, readout_noise_model):
    """Injects readout noise to the circuits in a quantum program."""
    readout_noise_map = {}
    for item in program.items:
        readout_tag = item.circuit[-1].label[
            -4:
        ]  # extract the readout layer tag from the circuit
        readout_noise_map[readout_tag] = readout_noise_model
        item.circuit = PassManager(
            [InsertNoisePass(noise_dict=readout_noise_map, noise_after=False)]
        ).run(item.circuit)

num_randomizations = 300
shots_per_randomization = 64

boxing_options = {"add_tags": "unique_box"}

# find the unique layers to generate the layered noise for
box_layers = find_combined_unique_layers(
    [mirrored_circuit], custom_boxing_options=boxing_options
)
layer_noise_model, readout_noise_model = generate_random_pauli_lindblad_noise(
    box_layers, num_qubits, seed=42
)

# a local noisy simulator in place of a QPU backend
noisy_executor = Executor(
    AerSimulator(method="stabilizer"),
    options={"simulator": {"layer_noise_model": layer_noise_model}},
)

task = MitigationTask()
program = task.prepare(
    circuit=mirrored_circuit,
    observables=observables,
    parameters=None,
    num_randomizations=num_randomizations,
    shots_per_randomization=shots_per_randomization,
    custom_boxing_options=boxing_options,
)
# For the simulation - we inject readout noise to the circuits of the program
add_readout_noise(program, readout_noise_model)

program._semantic_role = "simulation"  # flag the program for the simulator
unmitigated_res = task.postprocess(noisy_executor.run(program).result())

x_vals = np.array(range(num_qubits))
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, [1] * len(x_vals), "k-", alpha=0.5, label="Exact")
ax.errorbar(
    x_vals,
    unmitigated_res.data.evs,
    unmitigated_res.data.stds,
    fmt="o",
    label="Unmitigated",
)

x_tick_labels = [r"$Z_{ {i} }$" for i in range(num_qubits)]
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel("Expectation value")
ax.legend()
ax.set_xticks(x_vals, x_tick_labels)
plt.tight_layout()
plt.show()

Nauč se šum provázávajících vrstev

PEC začíná tím, že se naučí twirlovaný model šumu každé vrstvy provázávajících hradel v obvodu. Pokud předpokládáš markovovský model šumu, předpokládá se, že hloubka obvodu neovlivňuje šumový kanál dané provázávající vrstvy. Proto je třeba naučit pouze model šumu každé jedinečné vrstvy provázávajících hradel. Obvod použitý v tomto návodu lze rozdělit na dvě jedinečné vrstvy. Navrhovat obvody s minimální sadou jedinečných provázávajících vrstev je důležité, pokud se chceš učit model šumu, protože jedinečné vrstvy se musí učit nezávisle. Pokud má tvůj obvod příliš mnoho jedinečných vrstev, může být učení šumu neproveditelné kvůli nákladům na vzorkování při učení mnoha vrstev.

K učení šumu se používá qiskit-noise-learning. Learner opakuje každou vrstvu fragment_depths krát ve svých učicích obvodech, každý obvod twirluje num_randomizations krát a rozpad naměřených fidelit proloží Pauli-Lindbladovým modelem. Tyto parametry mohou ovlivnit přesnost naučeného modelu. Čím přesnější je naučený model šumu, tím přesnější bude potlačení šumu. Následující kód vykreslí naučený model šumu, přičemž obě vrstvy jsou sloučeny do jedné mapy. Každý qubit je koláč z jeho vah-1 sazeb X, Y a Z sečtených přes obě vrstvy a každá hrana ukazuje naučené sazby váhy 2. Například sazba 3e-4 znamená, že k této chybě dochází zhruba třikrát z každých 10 000 použití vrstvy.

Mitigace šumu odečtu se provádí v post-processingu pomocí TREX, na rozdíl od modelu šumu hradel, který je potřebný pro přípravu programu, aby se vložil anti-šum odpovídající naučenému modelu šumu. Šum odečtu proto lze buď spustit a naučit při učení šumu hradel, nebo při spouštění mitigovaného programu. Učení šumu odečtu spolu se šumem hradel může být přínosné pro naučení vnitřně konzistentního modelu šumu (viz například Disambiguating Pauli noise in quantum computers). V případech, kdy jsou šum hradel a šum odečtu v modelu šumu odděleny, doporučujeme učit šum odečtu při spouštění mitigovaného programu, aby se předešlo driftu šumu mezi okamžikem naučení a okamžikem spuštění.

V tomto návodu se šum odečtu učí během spouštění mitigovaného programu pomocí kalibračního obvodu TREX.

from qiskit.providers.fake_provider import GenericBackendV2
from qiskit.quantum_info import PauliList
from qiskit_ibm_runtime.results import LayerError, PauliLindbladError
from qiskit_ibm_runtime.visualization.embeddings import Embedding
from qiskit_mitigation import PEC
from qiskit_noise_learning.aer_executor import AerExecutor
from qiskit_noise_learning.noise_learner import LearningOptions, NoiseLearner
from samplomatic import InjectNoise

line_coupling = [(i, i + 1) for i in range(num_qubits - 1)]
backend = GenericBackendV2(
    num_qubits=num_qubits,
    basis_gates=["id", "rz", "sx", "x", "cz"],
    coupling_map=line_coupling,
)

task_pec = PEC()
# first learn only the layers PEC will inject noise into, i.e. not the measurement layer
pec_layers = [
    layer
    for layer in task_pec.find_unique_layers(mirrored_circuit, boxing_options)
    if get_annotation(layer.operation, InjectNoise)
]
layers = task_pec.find_unique_layers(mirrored_circuit, boxing_options)

learning_options = LearningOptions(
    fragment_depths=[2, 4, 8, 16],
    num_randomizations=250,
    shots_per_randomizations=16,
)

# qiskit-noise-learning is not compatible with the Executor in local mode,
# and requires the simulator present in the package for running simulation
noise_dict = {
    get_annotation(layer.operation, Tag).ref: noise_layer
    for layer, noise_layer in layer_noise_model
}
noise_learner_simulator = AerExecutor(
    qasm_simulator=AerSimulator(method="stabilizer"), noise_dict=noise_dict
)
learner = NoiseLearner(
    backend, options=learning_options, executor=noise_learner_simulator
)
learned_noise = (
    learner.run(pec_layers).result().to_dict()
)  # keyed by InjectNoise reference

# Draw the learned model on the line, laid out as a line with the qubits in order.
# Composing the two layers' maps sums their rates, giving one map over all the edges.
embedding = Embedding(
    [(0, 0), (0, 2), (0, 4), (2, 4), (4, 4), (4, 2), (4, 0), (2, 0)],
    line_coupling,
)
first, second = learned_noise.values()
learned_model = first.compose(second).simplify()
error = PauliLindbladError(
    PauliList(
        learned_model.get_qubit_sparse_pauli_list_copy().to_pauli_list()
    ),
    learned_model.rates,
)
fig = LayerError(circuit, list(range(num_qubits)), error).draw_map(
    embedding, radius=0.4, height=460, width=620
)
fig.update_traces(
    marker={
        "cmin": 0,
        "cmax": max(learned_model.rates),
        "colorbar": {"tickformat": ".1e", "title": "rate"},
    },
    selector={"mode": "markers"},
)
fig.show()

Připrav kombinovaný program PEC a TREX

Zde připravíme úkol TREX i úkol PEC a úkol TREX přidáme jako proměnnou k úkolu PEC. Tím se zajistí, že všechny vlastnosti potřebné pro mitigaci TREX (například twirling měření) se přidají do položky mitigované pomocí PEC. Vložený šum se vzorkuje jednou na randomizaci, takže statistická chyba PEC je řízena počtem randomizací, nikoli počtem shotů na randomizaci. Proto se používá mnoho randomizací - každá s několika shoty - a celkový počet shotů zůstává stejný jako u nemitigovaného programu.

Kalibrační obvod TREX se do programu přidává proto, aby se šum odečtu učil během spouštění programu. Kalibrační obvod TREX je velmi mělký, takže další randomizace pro něj jsou relativně levné a umožňují naučit šum odečtu přesněji.

from qiskit_mitigation import TREX

trex = TREX()

program_pec = task_pec.prepare(
    circuit=mirrored_circuit,
    observables=observables,
    parameters=None,
    num_randomizations=1200,
    shots_per_randomization=16,
    custom_boxing_options=boxing_options,
    trex=trex,
    noise_maps=learned_noise,  # PEC option
    noise_gain=0,  # PEC option
    max_sampling_overhead=100,  # PEC option
)
program_pec = trex.prepare(
    num_randomizations=1500,
    quantum_program=program_pec,
    custom_boxing_options=boxing_options,
)

# For the simulation - we inject readout noise to the circuits of the program
add_readout_noise(program_pec, readout_noise_model)

Prozkoumej režii shotů, spusť program a proveď post-processing

PEC vrací nezkreslený odhad střední hodnoty za cenu velké režie shotů. Požadovanou režii shotů lze aproximovat výpočtem gamma faktoru obvodu, který odráží vliv naučeného modelu šumu na obvod. Ve výchozím stavu se počet randomizací úkolu PEC škáluje přibližnou režií shotů potřebnou pro nezkreslený odhad (~$\gamma^2$). Chceš-li se tomuto automatickému škálování vyhnout, můžeš nastavit parametr scale_randomizations_by_gamma na False. Režie shotů je omezena parametrem max_sampling_overhead.

Před spuštěním bys měl zkontrolovat počet shotů, které PEC potřebuje k dosažení nezkresleného odhadu, a rozhodnout, zda je toto množství přijatelné, protože roste exponenciálně se sazbami šumu a hloubkou obvodu.

print(f"The circuit calculated gamma is: {task_pec.gamma}")
print(f"The circuit shot overhead is ~{task_pec.gamma**2}")

The circuit calculated gamma is: 1.2774846896089953
The circuit shot overhead is ~1.631967132185391

program_pec._semantic_role = (
    "simulation"  # flag the program for the simulator
)
job_pec_results = noisy_executor.run(program_pec).result()

# post-processing
pec_res = task_pec.postprocess(job_pec_results)
# In order to post process the results without TREX, we call `postprocess` again, but instead of the default options we
# pass only the PEC item result and the `measure_noise_data` as `None`
pec_without_trex_res = task_pec.postprocess(
    job_pec_results[task_pec._program_item_index], measure_noise_data=None
)

Porovnej mitigované střední hodnoty

Nakonec vykresli střední hodnotu každé pozorovatelné. PEC + TREX obnoví přesné střední hodnoty v rámci svých chybových úseček, zatímco nemitigovaný výsledek je silně zkreslen směrem dolů.

x = np.array(range(num_qubits))
x_tick_labels = [r"$Z_{ {i} }$" for i in range(num_qubits)]

fig, ax = plt.subplots(figsize=(8, 5))
ax.errorbar(
    x,
    unmitigated_res.data.evs,
    unmitigated_res.data.stds,
    fmt="o",
    label="Twirl only",
)
ax.errorbar(
    x,
    pec_without_trex_res.data.evs,
    pec_without_trex_res.data.stds,
    fmt="o",
    label="PEC only",
)
ax.errorbar(
    x,
    pec_res.data.evs,
    pec_res.data.stds,
    fmt="o",
    label="PEC + TREX",
)

ax.plot(x, [1] * len(x), "--", color="black", alpha=0.5, label="Exact")

tick_labels = [f"$Z_{ {i} }$" for i in range(num_qubits)]
ax.set_xticks(x, tick_labels)

ax.set_ylabel("Expectation value")
ax.legend()
ax.set_title("PEC+TREX Mitigated Expectation Values")
plt.tight_layout()
plt.show()