Estrapolazione a rumore zero con gate folding ed estrapolazione personalizzata

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. Questa guida amplifica il rumore con il "digital gate folding", il che significa che i Gate a due qubit vengono sostituiti da sequenze equivalenti del Gate e del suo inverso. Per esempio, sostituire un unitario $U$ con $U U^† U$ produrrebbe un fattore di amplificazione del rumore pari a 3. Puoi anche implementare l'amplificazione con l'amplificazione probabilistica degli errori (PEA), che richiede un passaggio aggiuntivo di apprendimento di un modello di rumore ed è trattata nella guida all'amplificazione probabilistica degli errori. Questa guida estrapola usando un paio di strategie di fitting integrate (lineare ed esponenziale), e mostra anche come implementare una propria funzione di fitting personalizzata.

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. Preparare un programma ZNE con digital gate folding e campionare circuiti con scale di rumore diverse
5. Estrapolare al limite di rumore zero con funzioni di fitting integrate e personalizzate
6. Confrontare i valori di aspettazione mitigati

Specificare un circuito e un'osservabile

Calcola 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.

# 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", "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 numpy as np
from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

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("Visualization of 1 Trotter step")
trotter_circuit(num_qubits, theta, layers, 1).draw("mpl", scale=0.7)

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$.

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

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. Imposti 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 da un Backend del 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 livello di Gate entangling, generato con random_pauli_lindblad_noise di qiskit-mitigation con un errore per Gate a livello (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()

Preparare un programma ZNE con digital gate folding e campionare circuiti con scale di rumore diverse

Prepara un programma ZNE contenente circuiti con scale di rumore crescenti. Ogni fattore di rumore rappresenta il multiplo dei Gate entangling che il circuito contiene (un fattore di rumore di 3.0 contiene tre volte più Gate entangling del circuito originale). Il circuito di ciascun fattore di rumore è sottoposto a twirling con lo stesso numero di randomizzazioni e shot del programma non mitigato.

from qiskit_mitigation import GateFolding

task_gf = GateFolding()
program_zne = task_gf.prepare(
    circuit=circuit,
    observables=observables,
    parameters=parameter_values,
    num_randomizations=num_randomizations,
    shots_per_randomization=shots_per_randomization,
    custom_boxing_options=boxing_options,
    noise_factors=[1, 3, 5],  # ZNE option
)

Disegna il primo passo di Trotter del template twirled con fattore di rumore 1 e con fattore di rumore 3 per visualizzare l'amplificazione. Per il fattore di rumore 3.0, ogni livello del circuito viene sostituito dal livello, dal suo inverso e di nuovo dal livello.

import IPython.display

def first_trotter_step(template, noise_factor):
    """The template up to the closing barrier of its first Trotter step (two layers, each folded ``noise_factor`` times)."""
    toy = template.copy_empty_like()
    for instruction in template.data:
        toy.append(instruction)
        if (
            instruction.operation.name == "barrier"
            and instruction.operation.label.startswith(
                f"R{2 * noise_factor - 1}@"
            )
        ):
            break
    return toy

IPython.display.display(
    first_trotter_step(program_zne.items[0].circuit, 1).draw(
        "mpl", scale=0.7, plot_barriers=False
    )
)
first_trotter_step(program_zne.items[1].circuit, 3).draw(
    "mpl", fold=-1, plot_barriers=False
)

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

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

Estrapolare al limite di rumore zero con funzioni di fitting integrate e personalizzate

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. Esistono molti modi per adattare i dati, e la scelta giusta della funzione dipende dalla natura dei dati campionati. GateFolding.postprocess fornisce diverse funzioni di estrapolazione integrate, tra cui linear, exponential e double_exponential. Il metodo postprocess accetta un elenco di metodi di estrapolazione, li valuta in ordine e restituisce il primo valore con valori fisicamente validi e varianza finita. Questo significa che l'ordine in cui specifichi gli estrapolatori ha un grande impatto su quale estrapolatore viene scelto. Passare fallback come estrapolatore significa che, se non viene trovata alcuna risposta valida, verrà scelto e restituito il valore di aspettazione del fattore di rumore più piccolo. Se fallback non viene passato e nessun estrapolatore produce un valore fisicamente valido, viene restituito il risultato con il std più piccolo. Sotto rumore di Pauli un'osservabile decade esponenzialmente con il fattore di rumore, quindi con 30 livelli (15 passi di Trotter, 2 livelli per passo) di Gate entangling l'estrapolatore esponenziale è la scelta naturale; puoi anche estrapolare linearmente, una scelta comune quando il decadimento è piccolo.

Per usare una funzione di fitting personalizzata, passa l'estrapolatore "custom" insieme a una tupla custom_fit contenente la funzione di fitting, un elenco con una stima iniziale per ciascuno dei suoi parametri e i limiti. Questi vengono inoltrati a scipy.optimize.curve_fit. Il post-processing può essere ripetuto con estrapolatori diversi senza rieseguire il programma. Qui usi una funzione di fitting personalizzata a radice quadrata. Per visualizzare l'estrapolazione, traccia la magnetizzazione mediata sui siti calcolata a ciascun fattore di rumore per $\theta = 0.15\pi$, i tre fit e la media sui siti dei valori di aspettazione estrapolati nel punto a rumore zero. I tre punti dati sembrano quasi lineari, eppure l'estrapolazione lineare finisce visibilmente al di sotto del valore esatto. Come previsto, il fit esponenziale stima nel modo più efficace il valore di aspettazione privo di rumore.

from scipy.optimize import curve_fit

zne_mitigated_exp = task_gf.postprocess(
    job_zne_results, extrapolator=["exponential", "linear"]
)
zne_mitigated_lin = task_gf.postprocess(
    job_zne_results, extrapolator=["linear"]
)

def sqrt_fit(x, a, b):
    return a * np.sqrt(x) + b

custom_fit = (sqrt_fit, [1.0, 1.0], (-np.inf, np.inf))
zne_mitigated_custom = task_gf.postprocess(
    job_zne_results,
    extrapolator=["custom", "exponential", "linear"],
    custom_fit=custom_fit,
)

# The result holds the expectation value at each noise factor and the extrapolated value, but
# not the fit parameters, so we refit the same models here to draw the curves.
noise_factors = np.array([1, 3, 5])
angle = 3  # theta = 0.15 pi, the angle with the largest unmitigated error
amplifications = zne_mitigated_exp.data.evs_noise_factors[:, angle, :].mean(
    axis=0
)
x_fit = np.linspace(0, 5, 51)
fits = {
    "exponential": (
        lambda x, a, b: a * np.exp(b * x),
        zne_mitigated_exp,
        "black",
    ),
    "linear": (lambda x, a, b: a * x + b, zne_mitigated_lin, "blue"),
    "square-root": (sqrt_fit, zne_mitigated_custom, "green"),
}

plt.plot(
    noise_factors, amplifications, "o", color="red", label="amplifications"
)
for name, (fit, res, color) in fits.items():
    params, _ = curve_fit(fit, noise_factors, amplifications, p0=[1.0, -0.1])
    plt.plot(
        x_fit, fit(x_fit, *params), "--", color=color, label=f"{name} fit"
    )
    plt.plot(
        [0],
        res.data.evs[:, angle].mean(),
        "*",
        color=color,
        markersize=15,
        label=f"{name} 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()

Confrontare i valori di aspettazione mitigati

Traccia la magnetizzazione mediata sui siti di ciascun metodo al variare degli angoli di rotazione di Ising. L'estrapolazione esponenziale recupera la magnetizzazione esatta entro le sue barre di errore, mentre l'estrapolazione lineare elimina gran parte del bias ma sottostima costantemente, perché il vero decadimento è esponenziale. Il fit a radice quadrata sovrastima in modo marcato.

results = {
    "Unmitigated": unmitigated_res,
    "ZNE exponential": zne_mitigated_exp,
    "ZNE linear": zne_mitigated_lin,
    "ZNE square-root": zne_mitigated_custom,
}

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