Estrapolazione a rumore zero con amplificazione probabilistica degli errori

L'estrapolazione a rumore zero (ZNE) è una tecnica per mitigare il rumore quantistico nel calcolo dei valori di aspettazione di osservabili. Funziona eseguendo un circuito quantistico a scale di rumore crescenti, calcolando i valori di aspettazione rispetto a tali scale di rumore ed estrapolando per stimare un valore di aspettazione a rumore zero. Sia la fase di amplificazione del rumore sia quella di estrapolazione possono essere implementate in modi diversi. In questa guida amplificherai il rumore con l'amplificazione probabilistica degli errori (PEA), come suggerito in "Evidence for the utility of quantum computing before fault tolerance". Per amplificare con la PEA, devi prima modellare il rumore che affligge i Gate entangling del circuito. Puoi anche implementare l'amplificazione con il digital gate folding, che non richiede alcun modello di rumore ed è trattato nella guida al gate folding. Apprenderai il rumore con qiskit-noise-learning ed estrapolerai con la strategia di fitting esponenziale integrata.

Tieni presente che, sebbene la ZNE possa essere efficace nel mitigare il rumore quantistico con un costo di campionamento inferiore rispetto ad altri metodi, non offre alcun limite sul bias.

Il flusso di lavoro è il seguente:

1. Specificare un circuito e un'osservabile
2. Calcolare i valori di aspettazione esatti usando il simulatore di vettori di stato di Qiskit
3. Calcolare i valori di aspettazione rumorosi e non mitigati
4. Apprendere il rumore dei livelli entangling con qiskit-noise-learning
5. Preparare un programma ZNE con amplificazione probabilistica degli errori e campionare circuiti con scale di rumore diverse
6. Estrapolare al limite di rumore zero
7. Confrontare i valori di aspettazione mitigati

Specificare un circuito e un'osservabile

In questa guida calcolerai la magnetizzazione mediata sui siti $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ di un modello di Ising a campo trasverso (TFIM) trotterizzato su otto qubit in una topologia ad anello. Farai variare diversi valori dell'intensità del campo trasverso e traccerai la magnetizzazione in funzione dell'intensità del campo. In una topologia ad anello su un numero pari di qubit, i livelli di Gate a due qubit possono essere suddivisi in due livelli unici, come mostrano i colori rosso e verde nella mappa di accoppiamento qui sotto. Progettare circuiti con un insieme minimo di livelli entangling unici è importante se vuoi apprendere un modello di rumore, perché i livelli unici devono essere appresi in modo indipendente. Se hai troppi livelli unici nel tuo circuito, l'apprendimento del rumore può diventare impraticabile a causa del costo di campionamento necessario per apprendere molti livelli.

# 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

Calcolare i valori di aspettazione esatti usando il simulatore di vettori di stato di Qiskit

Usa la classe Statevector di Qiskit per calcolare il valore di aspettazione esatto di $Z_i$ su ciascun sito per ogni valore di $\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()

Calcolare i valori di aspettazione rumorosi e non mitigati

Un QuantumProgram contiene un circuito quantistico e un samplex in un formato adatto alla primitiva Executor. MitigationTask.prepare costruisce un QuantumProgram a partire dal circuito, dalle osservabili e dai valori dei parametri, e racchiude il circuito in box con le opzioni di boxing indicate. Imposta l'opzione di boxing measure_annotations su "change_basis", che ruota le misure nelle basi richieste, e aggiungi add_tags: "unique_box" in modo che il simulatore possa mappare un modello di rumore sui box del circuito. Il programma è sottoposto a Pauli twirling per impostazione predefinita, il che su hardware trasforma il rumore coerente nel rumore di Pauli assunto dai metodi di mitigazione. num_randomizations imposta quante istanze twirled di ciascun circuito vengono campionate, e shots_per_randomization imposta quanti shot riceve ciascuna istanza, quindi il loro prodotto è il budget di shot per circuito.

Un utente con accesso a una QPU eseguirebbe il programma con un Executor costruito a partire da un backend di QiskitRuntimeService. Qui passerai un AerSimulator all'Executor, che esegue il programma in locale, e userai la sua opzione di simulatore layer_noise_model per iniettare un modello di rumore di Pauli-Lindblad casuale su ogni layer di gate di entanglement, generato con random_pauli_lindblad_noise di qiskit-mitigation con un errore per gate a livello di layer (EPLG) di 0.003, tipico delle attuali 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()

Apprendi il rumore dei layer di entanglement con qiskit-noise-learning

PEA inizia apprendendo il modello di rumore twirled di ciascun layer di gate di entanglement nel circuito. Assumiamo un modello di rumore markoviano e quindi occorre apprendere solo il modello di rumore di ciascun layer univoco di gate di entanglement, perché si assume che la profondità del circuito non influenzi il canale di rumore di un dato layer di entanglement. Più accurato è il modello di rumore appreso, più accurata sarà l'amplificazione.

Usiamo qiskit-noise-learning per apprendere il rumore. Il learner ripete ogni layer fragment_depths volte nei suoi circuiti di apprendimento, applica il twirling a ciascun circuito num_randomizations volte e adatta il decadimento delle fedeltà misurate a un modello di Pauli-Lindblad. Il codice seguente disegna il modello di rumore appreso, con i due layer uniti in un'unica mappa. Ogni qubit è un grafico a torta dei suoi tassi X, Y e Z di peso 1 sommati sui due layer, e ogni arco mostra i tassi di peso 2 appresi; un tasso di 3e-4 significa che l'errore si verifica circa tre volte ogni 10.000 utilizzi del layer.

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

Prepara un programma ZNE con amplificazione probabilistica degli errori e circuiti campione con scale di rumore variabili

Qui preparerai un programma ZNE contenente circuiti con scale di rumore crescenti. Con PEA, ogni fattore di rumore rappresenta il multiplo del rumore appreso che il circuito subisce: un fattore di rumore di 2.0 inietta una copia aggiuntiva del canale di rumore appreso di ciascun layer dopo quel layer, in aggiunta al rumore che il layer ha già. Poiché PEA scala direttamente i tassi di rumore appresi anziché ripetere i gate, è possibile qualsiasi fattore di rumore, compresi quelli non interi, e il circuito stesso non viene modificato: gli elementi del programma condividono un unico circuito modello e differiscono solo per il rumore iniettato dal samplex. Il rumore iniettato viene campionato una volta per randomizzazione, quindi l'errore statistico di PEA è controllato dal numero di randomizzazioni anziché dagli shot per randomizzazione; usiamo perciò molte randomizzazioni con pochi shot ciascuna, mantenendo il numero totale di shot uguale a quello del programma non mitigato. È importante non usare troppo poche randomizzazioni, altrimenti potresti incorrere in un bias di campionamento nell'iniezione del rumore.

L'overhead di campionamento totale scala linearmente con il numero di fattori di rumore.

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

Infine, esegui il programma sul simulatore rumoroso. I suoi risultati alimentano l'estrapolazione del passaggio successivo.

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

Estrapola al limite di rumore zero

Ora che hai i dati per ogni scala di rumore, puoi adattare una funzione ai dati ed estrapolare fino al valore a rumore zero per stimare un valore di aspettazione dell'osservabile privo di rumore. Con il rumore di Pauli un'osservabile decade esponenzialmente con il fattore di rumore, quindi sceglierai l'estrapolatore esponenziale come prima scelta e ripiegherai su quello lineare se l'adattamento esponenziale fallisce. Per visualizzare l'estrapolazione, rappresenta graficamente la magnetizzazione mediata sui siti calcolata a ciascun fattore di rumore, l'adattamento esponenziale e la media sui siti dei valori di aspettazione estrapolati nel punto a rumore zero. Vediamo che per questa particolare osservabile l'adattamento sottostima leggermente.

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

Confronta i valori di aspettazione mitigati

Rappresenta graficamente la magnetizzazione mediata sui siti di ciascun metodo al variare degli angoli di rotazione di Ising. PEA recupera la magnetizzazione esatta entro le sue barre di errore, mentre il risultato non mitigato è affetto da un bias che lo porta molto al di sotto.

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