Extrapolace k nulovému šumu s pravděpodobnostním zesilováním chyb

Extrapolace k nulovému šumu (ZNE) je technika pro zmírňování kvantového šumu při výpočtech střední hodnoty pozorovatelných. Funguje tak, že spustí kvantový obvod při rostoucích úrovních šumu, vezme střední hodnoty vzhledem k těmto úrovním šumu a extrapoluje, aby odhadla střední hodnotu při nulovém šumu. Fáze zesilování šumu i extrapolace lze implementovat různými způsoby. V tomto průvodci zesílíš šum pomocí pravděpodobnostního zesilování chyb (PEA), jak je navrženo v "Evidence for the utility of quantum computing before fault tolerance". Pro zesílení pomocí PEA musíš nejprve namodelovat šum ovlivňující provázávací hradla obvodu. Zesilování můžeš implementovat také pomocí digitálního hradlového foldingu, který nepotřebuje model šumu a je popsán v průvodci hradlovým foldingem. Šum se naučíš pomocí qiskit-noise-learning a extrapolovat budeš vestavěnou exponenciální fitovací strategií.

Pamatuj, že ačkoli ZNE může být účinná při zmírňování kvantového šumu za nižší cenu vzorkování než jiné metody, nemá žádné meze na zkreslení (bias).

Pracovní postup je následující:

1. Zadej obvod a pozorovatelnou
2. Vypočítej přesné střední hodnoty pomocí simulátoru stavového vektoru v Qiskitu
3. Vypočítej zašuměné, nezmírněné střední hodnoty
4. Nauč se šum provázávacích vrstev pomocí qiskit-noise-learning
5. Připrav program ZNE s pravděpodobnostním zesilováním chyb a vzorkuj obvody s různými úrovněmi šumu
6. Extrapoluj k limitě nulového šumu
7. Porovnej zmírněné střední hodnoty

Zadání obvodu a pozorovatelné

V tomto průvodci vypočítáš magnetizaci zprůměrovanou přes uzly $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ trotterizovaného Isingova modelu s příčným polem (TFIM) na osmi qubitech v prstencové topologii. Projdeš několik hodnot síly příčného pole a vykreslíš magnetizaci jako funkci síly pole. V prstencové topologii se sudým počtem qubitů lze vrstvy dvouqubitových hradel rozdělit na dvě unikátní vrstvy, jak ukazují červená a zelená barva na mapě propojení níže. Navrhovat obvody tak, aby měly minimální množinu unikátních provázávacích vrstev, je důležité, pokud se chceš učit model šumu, protože unikátní vrstvy se musí učit nezávisle. Pokud má tvůj obvod příliš mnoho unikátních vrstev, může být učení šumu neproveditelné kvůli ceně vzorkování při učení mnoha vrstev.

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

_needed = {"IPython": "IPython", "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", "scipy": "scipy"}
_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 IPython.display
import numpy as np
from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.visualization import plot_coupling_map

warnings.filterwarnings("ignore")

def trotter_circuit(num_qubits, theta, layers, num_steps):
    """Trotterized transverse-field Ising evolution: RX on every site, then a CZ layer per coupling layer."""
    circuit = QuantumCircuit(num_qubits)
    for _ in range(num_steps):
        circuit.rx(theta, range(num_qubits))
        for layer in layers:
            circuit.barrier()
            for edge in layer:
                circuit.cz(*edge)
                circuit.sdg(edge)
            circuit.barrier()
    return circuit

num_qubits = 8
trotter_steps = 15
ring_coupling = [(i, (i + 1) % num_qubits) for i in range(num_qubits)]
layers = [
    ring_coupling[0::2],
    ring_coupling[1::2],
]  # two disjoint layers of CZ gates
theta = Parameter("theta")
circuit = trotter_circuit(num_qubits, theta, layers, trotter_steps)

# Parameter values: theta in [0, 0.25 pi] - values for the transverse-field sweep
data_points = 6
parameter_values = np.linspace(0, 0.25 * np.pi, data_points).reshape(
    (data_points, 1)
)

# Define observables - Z on each site
observables = [
    SparsePauliOp("I" * (num_qubits - q - 1) + "Z" + "I" * q)
    for q in range(num_qubits)
]

print("Edges of the same color may execute their gates simultaneously.")
IPython.display.display(
    plot_coupling_map(
        num_qubits=num_qubits,
        qubit_coordinates=None,
        coupling_map=ring_coupling,
        figsize=(4, 4),
        qubit_size=150,
        line_width=5,
        line_color=[
            "red" if edge in layers[0] else "green" for edge in ring_coupling
        ],
        font_size=90,
    )
)
print("Visualization of 1 Trotter step")
trotter_circuit(num_qubits, theta, layers, 1).draw("mpl", scale=0.7)

Edges of the same color may execute their gates simultaneously.

Visualization of 1 Trotter step

Výpočet přesných středních hodnot pomocí simulátoru stavového vektoru v Qiskitu

Použij třídu Statevector z Qiskitu k výpočtu přesné střední hodnoty $Z_i$ na každém uzlu pro každou hodnotu $\theta$.

from qiskit.quantum_info import Statevector
import matplotlib.pyplot as plt

exact_evs = np.array(
    [
        [
            Statevector(circuit.assign_parameters([value]))
            .expectation_value(obs)
            .real
            for value in parameter_values.ravel()
        ]
        for obs in observables
    ]
)  # shape (observables, parameter values)

x_vals = parameter_values.ravel() / np.pi
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "ko-", label="Exact")
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
plt.tight_layout()
plt.show()

Výpočet zašuměných, nezmírněných středních hodnot

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 ohraničí do boxů podle zadaných voleb boxingu. Nastav volbu boxingu measure_annotations na "change_basis", která otočí měření 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 nastavení Pauliho twirlován, což na hardwaru mění koherentní šum na Pauliho šum, který zmírňovací 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 rozpočet shotů na obvod.

Uživatel s přístupem k QPU by program spustil pomocí Executor vytvořeného z backendu QiskitRuntimeService. Zde předáš Executor simulátor AerSimulator, který program spouští lokálně, a pomocí jeho simulátorové volby layer_noise_model vložíš náhodný Pauli-Lindbladův model šumu do každé vrstvy provázávajících hradel. Model je vygenerován pomocí random_pauli_lindblad_noise z qiskit-mitigation s chybou na vrstvené hradlo (EPLG) 0.003, což je typické pro současné QPU Heron.

from qiskit.circuit import CircuitInstruction
from qiskit.quantum_info import PauliLindbladMap, pauli_basis
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import Executor
from qiskit_mitigation import MitigationTask, find_combined_unique_layers

def generate_random_pauli_lindblad_noise(
    layers,
    num_qubits: int,
    eplg: float = 3e-3,
    seed: int | None = None,
) -> list[tuple[CircuitInstruction, 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 contain_measurement:  # do not add measurement noise
            layer_noise = PauliLindbladMap.identity(num_qubits)
        else:
            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))
    return layer_noise_model

num_randomizations = 50
shots_per_randomization = 64

boxing_options = {
    "measure_annotations": "change_basis",
    "add_tags": "unique_box",
}

# find the unique layers to generate the layered noise for
box_layers = find_combined_unique_layers(
    [circuit], custom_boxing_options=boxing_options
)
layer_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(),
    options={"simulator": {"layer_noise_model": layer_noise_model}},
)

task = MitigationTask()
program = task.prepare(
    circuit=circuit,
    observables=observables,
    parameters=parameter_values,
    num_randomizations=num_randomizations,
    shots_per_randomization=shots_per_randomization,
    custom_boxing_options=boxing_options,
)
program._semantic_role = "simulation"  # flag the program for the simulator
unmitigated_res = task.postprocess(noisy_executor.run(program).result())

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "k-", label="Exact")
ax.errorbar(
    x_vals,
    np.mean(unmitigated_res.data.evs, axis=0),
    np.mean(unmitigated_res.data.stds, axis=0),
    fmt="o-",
    label="Unmitigated",
)
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
plt.tight_layout()
plt.show()

Nauč se šum provázávajících vrstev pomocí qiskit-noise-learning

PEA začíná tím, že se naučí twirlovaný model šumu každé vrstvy provázávajících hradel v obvodu. Předpokládáme markovovský model šumu, a proto je třeba naučit pouze model šumu každé jedinečné vrstvy provázávajících hradel, protože se předpokládá, že hloubka obvodu neovlivňuje šumový kanál dané provázávající vrstvy. Čím přesnější je naučený model šumu, tím přesnější bude zesílení.

K učení šumu používáme 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. 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; sazba 3e-4 znamená, že k této chybě dochází zhruba třikrát z každých 10 000 použití vrstvy.

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 PEA
from qiskit_noise_learning.aer_executor import AerExecutor
from qiskit_noise_learning.noise_learner import LearningOptions, NoiseLearner
from samplomatic import InjectNoise, Tag
from samplomatic.utils import get_annotation

backend = GenericBackendV2(
    num_qubits=num_qubits,
    basis_gates=["id", "rz", "sx", "x", "cz"],
    coupling_map=ring_coupling,
)

task_pea = PEA()
# learn only the layers PEA will inject noise into, i.e. not the measurement layer
pea_layers = [
    layer
    for layer in task_pea.find_unique_layers(circuit, boxing_options)
    if get_annotation(layer.operation, InjectNoise)
]

learning_options = LearningOptions(
    fragment_depths=[2, 4, 8, 16],
    num_randomizations=50,
    shots_per_randomizations=64,
)
# 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(), noise_dict=noise_dict
)
learner = NoiseLearner(
    backend, options=learning_options, executor=noise_learner_simulator
)
learned_noise = (
    learner.run(pea_layers).result().to_dict()
)  # keyed by InjectNoise reference

# Draw the learned model on the ring, laid out as a square with the qubits in ring order.
# Composing the two layers' maps sums their rates, giving one map over all the ring's edges.
embedding = Embedding(
    [(0, 0), (0, 2), (0, 4), (2, 4), (4, 4), (4, 2), (4, 0), (2, 0)],
    ring_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 program ZNE s pravděpodobnostním zesílením chyb a vzorkuj obvody s různými úrovněmi šumu

Zde připravíš program ZNE obsahující obvody s rostoucími úrovněmi šumu. U PEA každý faktor šumu představuje násobek naučeného šumu, který obvod zažívá: faktor šumu 2.0 vloží za každou vrstvu jednu další kopii jejího naučeného šumového kanálu, navíc k šumu, který vrstva už má. Protože PEA škáluje naučené sazby šumu přímo, místo aby opakovala hradla, je možný jakýkoli faktor šumu, včetně neceločíselných, a samotný obvod se nemění: položky programu sdílejí jeden šablonový obvod a liší se pouze šumem vloženým pomocí samplexu. Vložený šum se vzorkuje jednou na randomizaci, takže statistická chyba PEA je řízena počtem randomizací, nikoli počtem shotů na randomizaci; proto používáme mnoho randomizací s několika shoty a celkový počet shotů ponecháváme stejný jako u nemitigovaného programu. Je důležité nepoužít příliš málo randomizací, jinak můžeš narazit na zkreslení vzorkování při vkládání šumu.

Celková režie vzorkování roste lineárně s počtem faktorů šumu.

noise_factors = [1, 2, 3]
program_pea = task_pea.prepare(
    circuit=circuit,
    observables=observables,
    parameters=parameter_values,
    num_randomizations=200,
    shots_per_randomization=16,
    custom_boxing_options=boxing_options,
    noise_maps=learned_noise,
    noise_factors=noise_factors,
)

Nakonec spusť program na zašumělém simulátoru. Jeho výsledky vstupují do extrapolace v dalším kroku.

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

Extrapoluj k limitě nulového šumu

Nyní, když máš data pro každou úroveň šumu, můžeš do dat proložit funkci a extrapolovat zpět k hodnotě při nulovém šumu, abys odhadl střední hodnotu bezšumové pozorovatelné. Při Pauliho šumu pozorovatelná exponenciálně klesá s faktorem šumu, takže jako první volbu zvolíš exponenciální extrapolátor a pokud exponenciální proložení selže, přejdeš na lineární. Pro vizualizaci extrapolace vykresli magnetizaci zprůměrovanou přes místa vypočtenou při každém faktoru šumu, exponenciální proložení a průměr extrapolovaných středních hodnot přes místa v bodě nulového šumu. Vidíme, že proložení pro tuto konkrétní pozorovatelnou mírně podstřeluje.

from scipy.optimize import curve_fit

pea_mitigated_res = task_pea.postprocess(
    job_pea_results, extrapolator=["exponential", "linear"]
)

# The result holds the expectation value at each noise factor and the extrapolated value for
# every observable; we average them over the sites and fit the same model to draw the curve.
angle = 3  # theta = 0.15 pi, the angle with the largest unmitigated error
amplifications = pea_mitigated_res.data.evs_noise_factors[:, angle, :].mean(
    axis=0
)
x_fit = np.linspace(0, max(noise_factors), 51)
params, _ = curve_fit(
    lambda x, a, b: a * np.exp(b * x),
    noise_factors,
    amplifications,
    p0=[1.0, -0.1],
)

plt.plot(
    noise_factors, amplifications, "o", color="red", label="amplifications"
)
plt.plot(
    x_fit,
    params[0] * np.exp(params[1] * x_fit),
    "--",
    color="black",
    label="exponential fit",
)
plt.plot(
    [0],
    pea_mitigated_res.data.evs[:, angle].mean(),
    "*",
    color="black",
    markersize=15,
    label="extrapolation",
)
plt.plot(
    [0],
    exact_evs[:, angle].mean(),
    "_",
    color="red",
    markersize=15,
    markeredgewidth=3,
    label="exact",
)
plt.xlabel("noise factor")
plt.ylabel(r"$\langle Z \rangle_{\text{avg}}$")
plt.legend()
plt.show()

Porovnej mitigované střední hodnoty

Vykresli magnetizaci zprůměrovanou přes místa pro každou metodu napříč sadou rotačních úhlů Isingova modelu. PEA obnoví přesnou magnetizaci v rámci svých chybových úseček, zatímco nemitigovaný výsledek je silně zkreslen směrem dolů.

results = {
    "Unmitigated": unmitigated_res,
    "PEA": pea_mitigated_res,
}

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(x_vals, exact_evs.mean(axis=0), "k-", label="Exact")
for label, res in results.items():
    ax.errorbar(
        x_vals,
        np.mean(res.data.evs, axis=0),
        np.mean(res.data.stds, axis=0),
        fmt="o-",
        label=label,
    )
ax.set_xlabel(r"$\theta/\pi$")
ax.set_ylabel(r"$\langle Z \rangle_{\text{avg}}$")
ax.legend()
ax.set_title(
    f"Error-Mitigated Magnetization [{num_qubits}-qubit, {trotter_steps} Trotter steps]"
)
plt.tight_layout()
plt.show()