Saltar al contenido principal

Simular el modelo de Ising 2D con campo inclinado mediante la función QESEM

nota

Las funciones de Qiskit son una característica experimental disponible únicamente para usuarios del Plan Premium, Plan Flex y Plan On-Prem (a través de la API de IBM Quantum Platform) de IBM Quantum®. Se encuentran en estado de versión preliminar y están sujetas a cambios.

Estimación de uso: 20 minutos en un procesador Heron r2. (NOTA: Esto es solo una estimación. Su tiempo de ejecución puede variar.)

Antecedentes

Este tutorial muestra cómo utilizar QESEM, la función de Qiskit de Qedma, para simular la dinámica de un modelo canónico de espín cuántico, el modelo de Ising 2D con campo inclinado (TFI) con ángulos no Clifford:

H=Ji,jZiZj+gxiXi+gziZi,H = J \sum_{\langle i,j \rangle} Z_i Z_j + g_x \sum_i X_i + g_z \sum_i Z_i ,

donde i,j\langle i,j \rangle denota los vecinos más cercanos en una red. Simular la evolución temporal de sistemas cuánticos de muchos cuerpos es una tarea computacionalmente difícil para las computadoras clásicas. Las computadoras cuánticas, por el contrario, están diseñadas de forma natural para realizar esta tarea de manera eficiente. El modelo TFI, en particular, se ha convertido en un punto de referencia popular en hardware cuántico debido a su rico comportamiento físico y su implementación compatible con el hardware.

En lugar de simular la dinámica de tiempo continuo, adoptamos el modelo de Ising pulsado, estrechamente relacionado. La dinámica puede expresarse exactamente como un circuito cuántico periódico, donde cada paso de evolución consiste en tres capas de compuertas fraccionales de dos qubits RZZ(αZZ)R_{ZZ} (\alpha_{ZZ}), intercaladas con capas de compuertas de un solo qubit RX(αX)R_X (\alpha_X) y RZ(αZ)R_Z (\alpha_Z).

Utilizaremos ángulos genéricos que son desafiantes tanto para la simulación clásica como para la mitigación de errores. Específicamente, elegimos αZZ=1.0\alpha_{ZZ} = 1.0, αX=0.53\alpha_X = 0.53 y αZ=0.1\alpha_Z = 0.1, situando el modelo lejos de cualquier punto integrable.

En este tutorial realizaremos lo siguiente:

  • Estimar el tiempo de ejecución esperado en la QPU para la mitigación completa de errores utilizando las funciones de estimación de tiempo analítica y empírica de QESEM.
  • Construir y simular el circuito del modelo de Ising 2D con campo inclinado utilizando diseños de qubits y capas de compuertas inspirados en el hardware.
  • Visualizar la conectividad de qubits del dispositivo y los subgrafos seleccionados para su experimento.
  • Demostrar el uso de la retropropagación de operadores (OBP) para reducir la profundidad del circuito. Esta técnica elimina operaciones del final del circuito a costa de más mediciones de operadores.
  • Realizar mitigación de errores (EM) sin sesgo para múltiples observables simultáneamente utilizando QESEM, comparando resultados ideales, ruidosos y mitigados.
  • Analizar y graficar el impacto de la mitigación de errores en la magnetización a diferentes profundidades de circuito.

Nota: OBP generalmente devolverá un conjunto de observables posiblemente no conmutantes. QESEM optimiza automáticamente las bases de medición cuando los observables objetivo contienen términos no conmutantes. Genera conjuntos candidatos de bases de medición utilizando varios algoritmos heurísticos y selecciona el conjunto que minimiza el número de bases distintas. Esto significa que QESEM agrupa observables compatibles en bases comunes para reducir el número total de configuraciones de medición requeridas, mejorando la eficiencia.

Acerca de QESEM

QESEM es un software confiable, de alta precisión y basado en caracterización que implementa una mitigación de errores cuasi-probabilística eficiente y sin sesgo. Está diseñado para mitigar errores en circuitos cuánticos genéricos y es independiente de la aplicación. Ha sido validado en diversas plataformas de hardware, incluyendo experimentos a escala de utilidad en dispositivos IBM® Eagle y Heron. Las etapas del flujo de trabajo de QESEM son las siguientes:

  1. Caracterización del dispositivo: mapea las fidelidades de las compuertas e identifica errores coherentes, proporcionando datos de calibración en tiempo real. Esta etapa garantiza que la mitigación aproveche las operaciones de mayor fidelidad disponibles.
  2. Transpilación consciente del ruido: genera y evalúa mapeos alternativos de qubits, conjuntos de operaciones y bases de medición, seleccionando la variante que minimiza el tiempo estimado de QPU, con paralelización opcional para acelerar la recopilación de datos.
  3. Supresión de errores: redefine las compuertas nativas, aplica entrelazado de Pauli y optimiza el control a nivel de pulso (en plataformas compatibles) para mejorar la fidelidad.
  4. Caracterización del circuito: construye un modelo de error local adaptado y lo ajusta a las mediciones de la QPU para cuantificar el ruido residual.
  5. Mitigación de errores: construye descomposiciones cuasi-probabilísticas de múltiples tipos y muestrea de ellas en un proceso adaptativo que minimiza el tiempo de mitigación en la QPU y la sensibilidad a las fluctuaciones del hardware, logrando altas precisiones en grandes volúmenes de circuitos.

Para obtener más información sobre QESEM y un experimento a escala de utilidad de este modelo en un subgrafo de 103 qubits con alta conectividad de la geometría nativa heavy-hex de ibm_marrakesh, consulta Reliable high-accuracy error mitigation for utility-scale quantum circuits. QESEM workflow.

Requisitos

Instala los siguientes paquetes de Python antes de ejecutar el notebook:

  • Qiskit SDK v2.0.0 o posterior (pip install qiskit)
  • Qiskit Runtime v0.40.0 o posterior (pip install qiskit-ibm-runtime)
  • Qiskit Functions Catalog v0.8.0 o posterior ( pip install qiskit-ibm-catalog )
  • Operator Backpropagation Qiskit addon v0.3.0 o posterior ( pip install qiskit-addon-obp )
  • Qiskit Utils addon v0.1.1 o posterior ( pip install qiskit-addon-utils )
  • Qiskit Aer simulator v0.17.1 o posterior ( pip install qiskit-aer )
  • Matplotlib v3.10.3 o posterior ( pip install matplotlib )

Configuración

Primero, importa las bibliotecas relevantes:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-addon-obp qiskit-addon-utils qiskit-aer qiskit-ibm-catalog qiskit-ibm-runtime
%matplotlib inline

from typing import Sequence

import matplotlib.pyplot as plt
import numpy as np

import qiskit
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_catalog import QiskitFunctionsCatalog
from qiskit_aer import AerSimulator
from qiskit_addon_utils.slicing import combine_slices, slice_by_gate_types
from qiskit_addon_obp import backpropagate
from qiskit_addon_obp.utils.simplify import OperatorBudget
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.visualization import (
plot_gate_map,
)

A continuación, autentíquese utilizando su clave de API del panel de control de IBM Quantum Platform. Luego, seleccione la función de Qiskit de la siguiente manera. (Tenga en cuenta que, por seguridad, es mejor guardar las credenciales de su cuenta en su entorno local, si se encuentra en una máquina de confianza, para que no tenga que ingresar su clave de API cada vez que se autentique.)

# Paste here your instance and token strings

instance = "YOUR_INSTANCE"
token = "YOUR_TOKEN"
channel = "ibm_quantum_platform"

catalog = QiskitFunctionsCatalog(
channel=channel, token=token, instance=instance
)
qesem_function = catalog.load("qedma/qesem")

Paso 1: Mapear las entradas clásicas a un problema cuántico

Comenzamos definiendo una función que crea el circuito de Trotter:

def trotter_circuit_from_layers(
steps: int,
theta_x: float,
theta_z: float,
theta_zz: float,
layers: Sequence[Sequence[tuple[int, int]]],
init_state: str | None = None,
) -> qiskit.QuantumCircuit:
"""
Generates an ising trotter circuit
:param steps: trotter steps
:param theta_x: RX angle
:param theta_z: RZ angle
:param theta_zz: RZZ angle
:param layers: list of layers (can be list of layers in device)
:param init_state: Initial state to prepare. If None, will not prepare any state. If "+", will
add Hadamard gates to all qubits.
:return: QuantumCircuit
"""
qubits = sorted({i for layer in layers for edge in layer for i in edge})
circ = qiskit.QuantumCircuit(max(qubits) + 1)

if init_state == "+":
print("init_state = +")
for q in qubits:
circ.h(q)

for _ in range(steps):
for q in qubits:
circ.rx(theta_x, q)
circ.rz(theta_z, q)

for layer in layers:
for edge in layer:
circ.rzz(theta_zz, *edge)
circ.barrier(qubits)

return circ

A continuación, creamos una función para calcular los valores de expectación ideales utilizando AerSimulator.

Ten en cuenta que para circuitos grandes (30 o más qubits) recomendamos utilizar valores precalculados a partir de simulaciones de propagación de creencias (BP) con PEPS. Este código incluye valores precalculados para 35 qubits como ejemplo, basados en el enfoque de BP para evolucionar una red tensorial PEPS introducido en este artículo (al que nos referimos como PEPS-BP), utilizando el paquete de Python para redes tensoriales quimb.

def calculate_ideal_evs(circ, obs, num_qubits, step):
# Predefined results for large circuits - calculated using bppeps for 3, 5, 7, 9 trotter steps
predefined_35 = [
0.79537,
0.78653,
0.79699,
]

if num_qubits == 35:
print(
"Using precalculated ideal values for large circuits calculated with belief propagation PEPS. Currently only for 35 qubits."
)
return predefined_35[step]

else:
simulator = AerSimulator()

# Use Estimator primitive to get expectation value
estimator = Estimator(simulator)
sim_result = estimator.run([(circ, [obs])], precision=0.0001).result()

# Extracting the result
ideal_values = sim_result[0].data.evs[0]
return ideal_values

Utilizamos un mapeo de capas RZZR_{ZZ} basado en el hardware tomado del dispositivo Heron, del cual recortamos las capas según el número de qubits que queremos simular. Definimos subgrafos para 10, 21, 28 y 35 qubits que mantienen una estructura 2D (siéntase libre de cambiar a su subgrafo favorito):

LAYERS_HERON_R2 = [  # the full set of hardware layers for Heron r2
[
(2, 3),
(6, 7),
(10, 11),
(14, 15),
(20, 21),
(16, 23),
(24, 25),
(17, 27),
(28, 29),
(18, 31),
(32, 33),
(19, 35),
(36, 41),
(42, 43),
(37, 45),
(46, 47),
(38, 49),
(50, 51),
(39, 53),
(60, 61),
(56, 63),
(64, 65),
(57, 67),
(68, 69),
(58, 71),
(72, 73),
(59, 75),
(76, 81),
(82, 83),
(77, 85),
(86, 87),
(78, 89),
(90, 91),
(79, 93),
(94, 95),
(100, 101),
(96, 103),
(104, 105),
(97, 107),
(108, 109),
(98, 111),
(112, 113),
(99, 115),
(116, 121),
(122, 123),
(117, 125),
(126, 127),
(118, 129),
(130, 131),
(119, 133),
(134, 135),
(140, 141),
(136, 143),
(144, 145),
(137, 147),
(148, 149),
(138, 151),
(152, 153),
(139, 155),
],
[
(1, 2),
(3, 4),
(5, 6),
(7, 8),
(9, 10),
(11, 12),
(13, 14),
(21, 22),
(23, 24),
(25, 26),
(27, 28),
(29, 30),
(31, 32),
(33, 34),
(40, 41),
(43, 44),
(45, 46),
(47, 48),
(49, 50),
(51, 52),
(53, 54),
(55, 59),
(61, 62),
(63, 64),
(65, 66),
(67, 68),
(69, 70),
(71, 72),
(73, 74),
(80, 81),
(83, 84),
(85, 86),
(87, 88),
(89, 90),
(91, 92),
(93, 94),
(95, 99),
(101, 102),
(103, 104),
(105, 106),
(107, 108),
(109, 110),
(111, 112),
(113, 114),
(120, 121),
(123, 124),
(125, 126),
(127, 128),
(129, 130),
(131, 132),
(133, 134),
(135, 139),
(141, 142),
(143, 144),
(145, 146),
(147, 148),
(149, 150),
(151, 152),
(153, 154),
],
[
(3, 16),
(7, 17),
(11, 18),
(22, 23),
(26, 27),
(30, 31),
(34, 35),
(21, 36),
(25, 37),
(29, 38),
(33, 39),
(41, 42),
(44, 45),
(48, 49),
(52, 53),
(43, 56),
(47, 57),
(51, 58),
(62, 63),
(66, 67),
(70, 71),
(74, 75),
(61, 76),
(65, 77),
(69, 78),
(73, 79),
(81, 82),
(84, 85),
(88, 89),
(92, 93),
(83, 96),
(87, 97),
(91, 98),
(102, 103),
(106, 107),
(110, 111),
(114, 115),
(101, 116),
(105, 117),
(109, 118),
(113, 119),
(121, 122),
(124, 125),
(128, 129),
(132, 133),
(123, 136),
(127, 137),
(131, 138),
(142, 143),
(146, 147),
(150, 151),
(154, 155),
(0, 1),
(4, 5),
(8, 9),
(12, 13),
(54, 55),
(15, 19),
],
]

subgraphs = { # the subgraphs for the different qubit counts such that it's 2D
10: list(range(22, 29)) + [16, 17, 37],
21: list(range(3, 12)) + list(range(23, 32)) + [16, 17, 18],
28: list(range(3, 12))
+ list(range(23, 32))
+ list(range(45, 50))
+ [16, 17, 18, 37, 38],
35: list(range(3, 12))
+ list(range(21, 32))
+ list(range(41, 50))
+ [16, 17, 18, 36, 37, 38],
42: list(range(3, 12))
+ list(range(21, 32))
+ list(range(41, 50))
+ list(range(63, 68))
+ [16, 17, 18, 36, 37, 38, 56, 57],
}

n_qubits = 35 # 21, 28, 35, 42
layers = [
[
edge
for edge in layer
if edge[0] in subgraphs[n_qubits] and edge[1] in subgraphs[n_qubits]
]
for layer in LAYERS_HERON_R2
]

print(layers)
[[(6, 7), (10, 11), (16, 23), (24, 25), (17, 27), (28, 29), (18, 31), (36, 41), (42, 43), (37, 45), (46, 47), (38, 49)], [(3, 4), (5, 6), (7, 8), (9, 10), (21, 22), (23, 24), (25, 26), (27, 28), (29, 30), (43, 44), (45, 46), (47, 48)], [(3, 16), (7, 17), (11, 18), (22, 23), (26, 27), (30, 31), (21, 36), (25, 37), (29, 38), (41, 42), (44, 45), (48, 49), (4, 5), (8, 9)]]

Ahora visualizamos la disposición de qubits en el dispositivo Heron para el subgrafo seleccionado:

service = QiskitRuntimeService(
channel=channel,
token=token,
instance=instance,
)
backend = service.backend("ibm_fez") # or any available device

selected_qubits = subgraphs[n_qubits]
num_qubits = backend.configuration().num_qubits
qubit_color = [
"#ff7f0e" if i in selected_qubits else "#d3d3d3"
for i in range(num_qubits)
]

plot_gate_map(
backend=backend,
figsize=(15, 10),
qubit_color=qubit_color,
)
plt.show()

Output of the previous code cell

Observa que la conectividad de la disposición de qubits elegida no es necesariamente lineal, y puede cubrir grandes regiones del dispositivo Heron dependiendo del número de qubits seleccionado.

Ahora generamos el circuito de Trotter y el observable de magnetización promedio para el número de qubits y los parámetros elegidos:

# Chosen parameters:
theta_x = 0.53
theta_z = 0.1
theta_zz = 1.0
steps = 9

circ = trotter_circuit_from_layers(steps, theta_x, theta_z, theta_zz, layers)
print(
f"Circuit 2q layers: {circ.depth(filter_function=lambda instr: len(instr.qubits) == 2)}"
)
print("\nCircuit structure:")

circ.draw("mpl", scale=0.8, fold=-1, idle_wires=False)
plt.show()

observable = qiskit.quantum_info.SparsePauliOp.from_sparse_list(
[("Z", [q], 1 / n_qubits) for q in subgraphs[n_qubits]],
np.max(subgraphs[n_qubits]) + 1,
) # Average magnetization observable

print(observable)
obs_list = [observable]
Circuit 2q layers: 27

Circuit structure:

Output of the previous code cell

SparsePauliOp(['IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'ZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIZIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII', 'IIIIIIIIIIIZIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII'],
coeffs=[0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j,
0.02857143+0.j, 0.02857143+0.j, 0.02857143+0.j])

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

Estimación del tiempo de QPU con y sin OBP

Los usuarios generalmente desean saber cuánto tiempo de QPU se requiere para su experimento. Sin embargo, esto se considera un problema difícil para las computadoras clásicas.

QESEM ofrece dos modos de estimación de tiempo para informar a los usuarios sobre la viabilidad de sus experimentos:

  1. Estimación de tiempo analítica: proporciona una estimación muy aproximada y no requiere tiempo de QPU. Puede utilizarse para probar si un paso de transpilación reduciría potencialmente el tiempo de QPU.
  2. Estimación de tiempo empírica (demostrada aquí): proporciona una estimación bastante precisa y utiliza unos pocos minutos de tiempo de QPU.

En ambos casos, QESEM genera la estimación de tiempo para alcanzar la precisión requerida para todos los observables.

run_on_real_hardware = True

precision = 0.05
if run_on_real_hardware:
backend_name = "ibm_fez"
else:
backend_name = "fake_fez"

# Start a job for empirical time estimation
estimation_job_wo_obp = qesem_function.run(
pubs=[(circ, obs_list)],
instance=instance,
backend_name=backend_name, # E.g. "ibm_brisbane"
options={
"estimate_time_only": "empirical", # "empirical" - gets actual time estimates without running full mitigation
"max_execution_time": 120, # Limits the QPU time, specified in seconds.
"default_precision": precision,
},
)
print(estimation_job_wo_obp.job_id)
print(estimation_job_wo_obp.status())
17d3828e-9fdb-482e-8e9b-392f3eefe313
DONE
# Get the result object (blocking method). Use job.status() in a loop for non-blocking.
# This takes 1-3 minutes
result = estimation_job_wo_obp.result()
print(
f"Empirical time estimation (sec): {result[0].metadata['time_estimation_sec']}"
)
Empirical time estimation (sec): 1200

Ahora utilizaremos la retropropagación de operadores (OBP). (Consulta la guía de OBP para obtener más detalles sobre el complemento OBP de Qiskit.) Crearemos una función que genera los segmentos del circuito para la retropropagación:

def run_backpropagation(circ_vec, observable, steps_vec, max_qwc_groups=8):
"""
Runs backpropagation for a list of circuits and observables.
Returns lists of backpropagated circuits and observables.
"""
op_budget = OperatorBudget(max_qwc_groups=max_qwc_groups)
bp_circuit_vec = []
bp_observable_vec = []

for i, circ in enumerate(circ_vec):
slices = slice_by_gate_types(circ)
bp_observable, remaining_slices, metadata = backpropagate(
observable,
slices,
operator_budget=op_budget,
)
bp_circuit = combine_slices(remaining_slices, include_barriers=True)
bp_circuit_vec.append(bp_circuit)
bp_observable_vec.append(bp_observable)
print(f"n.o. steps: {steps_vec[i]}")
print(f"Backpropagated {metadata.num_backpropagated_slices} slices.")
print(
f"New observable has {len(bp_observable.paulis)} terms, which can be combined into "
f"{len(bp_observable.group_commuting(qubit_wise=True))} groups.\n"
f"After truncation, the error in our observable is bounded by {metadata.accumulated_error(0):.3e}"
)
print("-----------------")
return bp_circuit_vec, bp_observable_vec

Llamamos a la función:

bp_circ_vec, bp_obs_vec = run_backpropagation([circ], observable, [steps])
n.o. steps: 9
Backpropagated 11 slices.
New observable has 363 terms, which can be combined into 4 groups.
After truncation, the error in our observable is bounded by 0.000e+00
-----------------
print("The remaining circuit after backpropagation looks as follows:")
bp_circ_vec[-1].draw("mpl", scale=0.8, fold=-1, idle_wires=False)
None
The remaining circuit after backpropagation looks as follows:

Output of the previous code cell

Podemos observar que la retropropagación redujo dos capas del circuito. Ahora que tenemos nuestro circuito reducido y los observables expandidos, realicemos la estimación de tiempo para el circuito retropropagado:

# Start a job for empirical time estimation
estimation_job_obp = qesem_function.run(
pubs=[(bp_circ_vec[-1], [bp_obs_vec[-1]])],
instance=instance,
backend_name=backend_name,
options={
"estimate_time_only": "empirical",
"max_execution_time": 120,
"default_precision": precision,
},
)
print(estimation_job_obp.job_id)
print(estimation_job_obp.status())
8bae699d-a16b-4d39-bbd9-d123fbcce55d
DONE
result_obp = estimation_job_obp.result()
print(
f"Empirical time estimation (sec): {result_obp[0].metadata['time_estimation_sec']}"
)
Empirical time estimation (sec): 900

Vemos que OBP reduce el costo de tiempo para la mitigación del circuito.

Paso 3: Ejecutar utilizando primitivas de Qiskit

Ejecutar con backend real

Ahora ejecutamos el experimento completo con varios pasos de Trotter. El número de qubits, la precisión requerida y el tiempo máximo de QPU pueden modificarse según los recursos de QPU disponibles. Ten en cuenta que restringir el tiempo máximo de QPU afectará la precisión final, como verás en el gráfico final a continuación.

Analizamos cuatro circuitos con 5, 7 y 9 pasos de Trotter a una precisión de 0.05, comparando sus valores de expectación ideales, ruidosos y mitigados con corrección de errores:

steps_vec = [5, 7, 9]

circ_vec = []
for steps in steps_vec:
circ = trotter_circuit_from_layers(
steps, theta_x, theta_z, theta_zz, layers
)
circ_vec.append(circ)

Nuevamente, realizamos OBP en cada circuito para reducir el tiempo de ejecución:

bp_circ_vec_35, bp_obs_vec_35 = run_backpropagation(
circ_vec, observable, steps_vec
)
n.o. steps: 5
Backpropagated 11 slices.
New observable has 363 terms, which can be combined into 4 groups.
After truncation, the error in our observable is bounded by 0.000e+00
-----------------
n.o. steps: 7
Backpropagated 11 slices.
New observable has 363 terms, which can be combined into 4 groups.
After truncation, the error in our observable is bounded by 0.000e+00
-----------------
n.o. steps: 9
Backpropagated 11 slices.
New observable has 363 terms, which can be combined into 4 groups.
After truncation, the error in our observable is bounded by 0.000e+00
-----------------

Ahora ejecutamos un lote de trabajos completos de QESEM. Limitamos el tiempo máximo de ejecución en la QPU para cada uno de los puntos para tener un mejor control del presupuesto de QPU.

run_on_real_hardware = True

precision = 0.05
if run_on_real_hardware:
backend_name = "ibm_marrakesh"
else:
backend_name = "fake_fez"
# Running full jobs for:
pubs_list = [
[(bp_circ_vec_35[i], bp_obs_vec_35[i])] for i in range(len(bp_obs_vec_35))
]
# Initiating multiple jobs for different lengths
job_list = []
for pubs in pubs_list:
job_obp = qesem_function.run(
pubs=pubs,
instance=instance,
backend_name=backend_name, # E.g. "ibm_brisbane"
options={
"max_execution_time": 300, # Limits the QPU time, specified in seconds.
"default_precision": 0.05,
},
)
job_list.append(job_obp)

Aquí verificamos el estado de cada trabajo:

for job in job_list:
print(job.status())
DONE
DONE
DONE
DONE

Paso 4: Postprocesar y devolver el resultado en el formato clásico deseado

Cuando todos los trabajos hayan terminado de ejecutarse, podemos comparar sus valores de expectación ruidosos y mitigados.

ideal_values = []
noisy_values = []
error_mitigated_values = []
error_mitigated_stds = []

for i in range(len(job_list)):
job = job_list[i]
result = job.result() # Blocking - takes 3-5 minutes
noisy_results = result[0].metadata["noisy_results"]

ideal_val = calculate_ideal_evs(circ_vec[i], observable, n_qubits, i)
print("---------------------------------")
print(f"Ideal: {ideal_val}")
print(f"Noisy: {noisy_results.evs}")
print(f"QESEM: {result[0].data.evs} \u00b1 {result[0].data.stds}")

ideal_values.append(ideal_val)
noisy_values.append(noisy_results.evs)
error_mitigated_values.append(result[0].data.evs)
error_mitigated_stds.append(result[0].data.stds)
Using precalculated ideal values for large circuits calculated with belief propagation PEPS. Currently only for 35 qubits.
---------------------------------
Ideal: 0.79537
Noisy: 0.7039237951821501
QESEM: 0.7828018244130982 ± 0.013257266977728376
Using precalculated ideal values for large circuits calculated with belief propagation PEPS. Currently only for 35 qubits.
---------------------------------
Ideal: 0.78653
Noisy: 0.6478583812958806
QESEM: 0.7875259197423828 ± 0.02703045139248604
Using precalculated ideal values for large circuits calculated with belief propagation PEPS. Currently only for 35 qubits.
---------------------------------
Ideal: 0.79699
Noisy: 0.6171787879868142
QESEM: 0.6918791909168913 ± 0.0740873782039517

Por último, podemos graficar la magnetización en función del número de pasos. Esto resume el beneficio de utilizar la función de Qiskit QESEM para la mitigación de errores sin sesgo en dispositivos cuánticos ruidosos.

plt.plot(steps_vec, ideal_values, "--", label="ideal")
plt.scatter(steps_vec, noisy_values, label="noisy")
plt.errorbar(
steps_vec,
error_mitigated_values,
yerr=error_mitigated_stds,
fmt="o",
capsize=5,
label="QESEM mitigation",
)
plt.legend()
plt.xlabel("n.o. steps")
plt.ylabel("Magnetization")
Text(0, 0.5, 'Magnetization')

Output of the previous code cell

nota

El noveno paso tiene una barra de error estadístico grande porque limitamos el tiempo de QPU a 5 minutos. Si ejecuta este paso durante 15 minutos (como sugiere la estimación de tiempo empírica), obtendrás una barra de error más pequeña. Por lo tanto, el valor mitigado se acercará más al valor ideal.