Extrapolação para ruído zero com amplificação probabilística de erros

A extrapolação para ruído zero (ZNE) é uma técnica para mitigar o ruído quântico em cálculos de valores esperados de observáveis. Ela funciona executando um circuito quântico em escalas de ruído crescentes, obtendo valores esperados para essas escalas de ruído e extrapolando para estimar um valor esperado sem ruído. Tanto a etapa de amplificação do ruído quanto a de extrapolação podem ser implementadas de diferentes maneiras. Neste guia, você vai amplificar o ruído com amplificação probabilística de erros (PEA), como sugerido em "Evidence for the utility of quantum computing before fault tolerance". Para amplificar com PEA, você precisa primeiro modelar o ruído que afeta as portas emaranhadoras do circuito. Você também pode implementar a amplificação com gate folding digital, que não precisa de modelo de ruído e é abordado no guia de gate folding. Você vai aprender o ruído com qiskit-noise-learning e extrapolar com a estratégia de ajuste exponencial pronta.

Tenha em mente que, embora a ZNE possa ser eficaz para mitigar o ruído quântico com um custo de amostragem menor do que outros métodos, ela não oferece limites para o viés.

O fluxo de trabalho é o seguinte:

1. Especificar um circuito e um observável
2. Calcular valores esperados exatos usando o simulador de vetor de estado do Qiskit
3. Calcular valores esperados ruidosos, sem mitigação
4. Aprender o ruído das camadas emaranhadoras com qiskit-noise-learning
5. Preparar um programa ZNE com amplificação probabilística de erros e circuitos de amostra com escalas de ruído variadas
6. Extrapolar para o limite de ruído zero
7. Comparar os valores esperados mitigados

Especificar um circuito e um observável

Neste guia, você vai calcular a magnetização média por sítio $\langle Z \rangle_{\text{avg}} = \frac{1}{N} \sum_{i=1}^{N} \langle Z_i \rangle$ de um modelo de Ising de campo transversal (TFIM) trotterizado em oito qubits em uma topologia de anel. Você vai varrer vários valores da intensidade do campo transversal e plotar a magnetização em função da intensidade do campo. Em uma topologia de anel com número par de qubits, as camadas de portas de dois qubits podem ser divididas em duas camadas únicas, como mostram as cores vermelha e verde no mapa de acoplamento abaixo. Projetar circuitos com um conjunto mínimo de camadas emaranhadoras únicas é importante se você quer aprender um modelo de ruído, porque as camadas únicas precisam ser aprendidas de forma independente. Se o seu circuito tiver camadas únicas demais, o aprendizado de ruído pode se tornar inviável devido ao custo de amostragem de aprender muitas camadas.

# 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 exatos usando o simulador de vetor de estado do Qiskit

Use a classe Statevector do Qiskit para calcular o valor esperado exato de $Z_i$ em cada sítio 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 ruidosos, sem mitigação

Um QuantumProgram contém um circuito quântico e um samplex em um formato adequado para a primitiva Executor. MitigationTask.prepare constrói um QuantumProgram a partir do circuito, dos observáveis e dos valores dos parâmetros, e envolve o circuito em caixas com as opções de boxing fornecidas. Defina a opção de boxing measure_annotations como "change_basis", que rotaciona as medições para as bases necessárias, e adicione add_tags: "unique_box" para que o simulador possa mapear um modelo de ruído nas caixas do circuito. O programa passa por Pauli twirling por padrão, o que, em hardware, transforma ruído coerente no ruído de Pauli que os métodos de mitigação assumem. num_randomizations define quantas instâncias twirled de cada circuito são amostradas, e shots_per_randomization define quantos shots cada instância recebe, então o produto dos dois é o orçamento de shots por circuito.

Um usuário com acesso a QPU executaria o programa com um Executor criado a partir de um backend do QiskitRuntimeService. Aqui você vai passar um AerSimulator para o Executor, que executa o programa localmente, e usar a opção de simulador layer_noise_model para injetar um modelo de ruído de Pauli-Lindblad aleatório em cada camada de portas de emaranhamento, gerado com random_pauli_lindblad_noise de qiskit-mitigation com um erro por porta em camada (EPLG) de 0.003, típico dos QPUs Heron atuais.

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

Aprenda o ruído das camadas de emaranhamento com qiskit-noise-learning

O PEA começa aprendendo o modelo de ruído com twirling de cada camada de portas de emaranhamento no circuito. Assumimos um modelo de ruído markoviano e, portanto, só o modelo de ruído de cada camada única de portas de emaranhamento precisa ser aprendido, porque se assume que a profundidade do circuito não afeta o canal de ruído de uma determinada camada de emaranhamento. Quanto mais preciso for o modelo de ruído aprendido, mais precisa será a amplificação.

Usamos qiskit-noise-learning para aprender o ruído. O aprendiz repete cada camada fragment_depths vezes em seus circuitos de aprendizado, aplica twirling a cada circuito num_randomizations vezes e ajusta o decaimento das fidelidades medidas a um modelo de Pauli-Lindblad. O código a seguir desenha o modelo de ruído aprendido, com as duas camadas combinadas em um único mapa. Cada qubit é um gráfico de pizza das suas taxas X, Y e Z de peso 1 somadas nas duas camadas, e cada aresta mostra as taxas de peso 2 aprendidas; uma taxa de 3e-4 significa que esse erro ocorre em aproximadamente três de cada 10.000 usos da camada.

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

Prepare um programa ZNE com amplificação probabilística de erro e circuitos de amostra com diferentes escalas de ruído

Aqui você vai preparar um programa ZNE contendo circuitos com escalas de ruído crescentes. Com o PEA, cada fator de ruído representa o múltiplo do ruído aprendido que o circuito experimenta: um fator de ruído de 2.0 injeta uma cópia adicional do canal de ruído aprendido de cada camada depois dessa camada, além do ruído que a camada já tem. Como o PEA escala diretamente as taxas de ruído aprendidas em vez de repetir portas, qualquer fator de ruído é possível, inclusive não inteiros, e o circuito em si não é modificado: os itens do programa compartilham um único circuito modelo e diferem apenas no ruído injetado pelo samplex. O ruído injetado é amostrado uma vez por randomização, então o erro estatístico do PEA é controlado pelo número de randomizações, e não pelos shots por randomização; por isso usamos muitas randomizações com poucos shots cada, mantendo o número total de shots igual ao do programa sem mitigação. É importante não usar randomizações de menos, pois caso contrário você pode encontrar viés de amostragem na injeção de ruído.

O overhead total de amostragem cresce linearmente com o número de fatores de ruído.

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 fim, execute o programa no simulador com ruído. Os resultados alimentam a extrapolação do próximo passo.

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

Extrapole para o limite de ruído zero

Agora que você tem dados para cada escala de ruído, pode ajustar uma função aos dados e extrapolar de volta até o valor de ruído zero para estimar o valor esperado do observável sem ruído. Sob ruído de Pauli, um observável decai exponencialmente com o fator de ruído, então você escolherá o extrapolador exponencial como primeira opção e recorrerá ao linear se o ajuste exponencial falhar. Para visualizar a extrapolação, plote a magnetização média por sítio calculada em cada fator de ruído, o ajuste exponencial e a média por sítio dos valores esperados extrapolados no ponto de ruído zero. Vemos que o ajuste fica ligeiramente abaixo para este observável em 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()

Compare os valores esperados mitigados

Plote a magnetização média por sítio de cada método ao longo da varredura de ângulos de rotação de Ising. O PEA recupera a magnetização exata dentro de suas barras de erro, enquanto o resultado sem mitigação fica bem abaixo dela.

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