Extrapolación a ruido cero con amplificación probabilística de errores

La extrapolación a ruido cero (ZNE) es una técnica para mitigar el ruido cuántico en el cálculo de valores esperados de observables. Funciona ejecutando un circuito cuántico con escalas de ruido crecientes, tomando valores esperados respecto a esas escalas de ruido y extrapolando para estimar un valor esperado a ruido cero. Tanto la amplificación del ruido como las etapas de extrapolación se pueden implementar de diferentes maneras. En esta guía amplificarás el ruido con amplificación probabilística de errores (PEA), como se sugiere en "Evidence for the utility of quantum computing before fault tolerance". Para amplificar con PEA, primero debes modelar el ruido que afecta a las puertas entrelazantes del circuito. También puedes implementar la amplificación con plegado digital de puertas, que no necesita un modelo de ruido y se explica en la guía de plegado de puertas. Aprenderás el ruido con qiskit-noise-learning y extrapolarás con la estrategia de ajuste exponencial integrada.

Ten en cuenta que, aunque la ZNE puede ser eficaz para mitigar el ruido cuántico con un coste de muestreo menor que otros métodos, no ofrece ninguna cota sobre el sesgo.

El flujo de trabajo es el siguiente:

1. Especificar un circuito y un observable
2. Calcular valores esperados exactos con el simulador de vector de estado de Qiskit
3. Calcular valores esperados con ruido y sin mitigar
4. Aprender el ruido de las capas entrelazantes con qiskit-noise-learning
5. Preparar un programa de ZNE con amplificación probabilística de errores y circuitos de muestra con distintas escalas de ruido
6. Extrapolar al límite de ruido cero
7. Comparar los valores esperados mitigados

Especificar un circuito y un observable

En esta guía calcularás la magnetización promediada por sitios $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ de un modelo de Ising de campo transversal (TFIM) trotterizado en ocho qubits con topología de anillo. Barrerás varios valores de la intensidad del campo transversal y representarás la magnetización en función de la intensidad del campo. En una topología de anillo con un número par de qubits, las capas de puertas de dos qubits se pueden dividir en dos capas únicas, como muestran los colores rojo y verde en el mapa de acoplamiento de abajo. Diseñar circuitos con un conjunto mínimo de capas entrelazantes únicas es importante si quieres aprender un modelo de ruido, porque las capas únicas deben aprenderse de forma independiente. Si tu circuito tiene demasiadas capas únicas, el aprendizaje de ruido puede volverse inviable por el coste de muestreo de aprender muchas capas.

# 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

Calcular valores esperados exactos con el simulador de vector de estado de Qiskit

Usa la clase Statevector de Qiskit para calcular el valor esperado exacto de $Z_i$ en cada sitio para cada valor de $\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()

Calcular valores esperados con ruido y sin mitigar

Un QuantumProgram contiene un circuito cuántico y un samplex en un formato adecuado para la primitiva Executor. MitigationTask.prepare construye un QuantumProgram a partir del circuito, los observables y los valores de los parámetros, y encierra el circuito en cajas con las opciones de boxing indicadas. Establece la opción de boxing measure_annotations en "change_basis", que rota las mediciones a las bases necesarias, y añade add_tags: "unique_box" para que el simulador pueda asignar un modelo de ruido a las cajas del circuito. El programa se somete a Pauli twirling de forma predeterminada, lo que en hardware convierte el ruido coherente en el ruido de Pauli que asumen los métodos de mitigación. num_randomizations establece cuántas instancias con twirling de cada circuito se muestrean, y shots_per_randomization establece cuántos shots recibe cada instancia, de modo que su producto es el presupuesto de shots por circuito.

Un usuario con acceso a QPU ejecutaría el programa con un Executor creado a partir de un backend de QiskitRuntimeService. Aquí pasarás un AerSimulator al Executor, que ejecuta el programa localmente, y usarás su opción de simulador layer_noise_model para inyectar un modelo de ruido de Pauli-Lindblad aleatorio en cada capa de puertas de entrelazamiento, generado con random_pauli_lindblad_noise de qiskit-mitigation con un error por puerta en capa (EPLG) de 0.003, típico de los QPU Heron actuales.

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

Aprende el ruido de las capas de entrelazamiento con qiskit-noise-learning

PEA comienza aprendiendo el modelo de ruido con twirling de cada capa de puertas de entrelazamiento del circuito. Asumimos un modelo de ruido markoviano y, por lo tanto, solo debe aprenderse el modelo de ruido de cada capa única de puertas de entrelazamiento, porque se supone que la profundidad del circuito no afecta al canal de ruido de una capa de entrelazamiento dada. Cuanto más preciso sea el modelo de ruido aprendido, más precisa será la amplificación.

Usamos qiskit-noise-learning para aprender el ruido. El aprendiz repite cada capa fragment_depths veces en sus circuitos de aprendizaje, aplica twirling a cada circuito num_randomizations veces y ajusta el decaimiento de las fidelidades medidas a un modelo de Pauli-Lindblad. El siguiente código dibuja el modelo de ruido aprendido, con las dos capas fusionadas en un solo mapa. Cada qubit es un gráfico circular con sus tasas X, Y y Z de peso 1 sumadas sobre las dos capas, y cada arista muestra las tasas de peso 2 aprendidas; una tasa de 3e-4 significa que ese error ocurre aproximadamente en tres de cada 10 000 usos de la capa.

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 programa ZNE con amplificación probabilística de errores y circuitos de muestra con distintas escalas de ruido

Aquí prepararás un programa ZNE que contiene circuitos con escalas de ruido crecientes. Con PEA, cada factor de ruido representa el múltiplo del ruido aprendido que experimenta el circuito: un factor de ruido de 2.0 inyecta una copia adicional del canal de ruido aprendido de cada capa después de esa capa, además del ruido que la capa ya tiene. Como PEA escala directamente las tasas de ruido aprendidas en lugar de repetir puertas, es posible cualquier factor de ruido, incluidos los no enteros, y el circuito en sí no se modifica: los elementos del programa comparten un circuito plantilla y solo se diferencian en el ruido inyectado por el samplex. El ruido inyectado se muestrea una vez por aleatorización, de modo que el error estadístico de PEA se controla mediante el número de aleatorizaciones y no mediante los shots por aleatorización; por ello usamos muchas aleatorizaciones con pocos shots cada una, manteniendo el número total de shots igual que en el programa sin mitigar. Es importante no usar muy pocas aleatorizaciones, ya que de lo contrario podrías encontrar sesgo de muestreo en la inyección de ruido.

La sobrecarga total de muestreo escala linealmente con el número de factores de ruido.

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

Por último, ejecuta el programa en el simulador con ruido. Sus resultados alimentan la extrapolación del siguiente paso.

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

Extrapola al límite de ruido cero

Ahora que tienes datos para cada escala de ruido, puedes ajustar una función a los datos y extrapolar hacia atrás hasta el valor de ruido cero para estimar el valor esperado del observable sin ruido. Con ruido de Pauli, un observable decae exponencialmente con el factor de ruido, por lo que elegirás el extrapolador exponencial como primera opción y recurrirás al lineal si el ajuste exponencial falla. Para visualizar la extrapolación, representa la magnetización promediada sobre los sitios calculada en cada factor de ruido, el ajuste exponencial y el promedio sobre los sitios de los valores esperados extrapolados en el punto de ruido cero. Vemos que el ajuste se queda ligeramente corto para este observable en particular.

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

Compara los valores esperados mitigados

Representa la magnetización promediada sobre los sitios de cada método a lo largo del barrido de ángulos de rotación de Ising. PEA recupera la magnetización exacta dentro de sus barras de error, mientras que el resultado sin mitigar queda sesgado muy por debajo.

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