Saltar al contenido principal

Experimento a escala de utilidad I

nota

Tamiya Onodera (5 de julio de 2024)

Descargar PDF de la clase original. Nota: Algunos fragmentos de código podrían estar desactualizados, ya que son imágenes estáticas.

El tiempo aproximado de QPU para este experimento es de 45 segundos.

1. Introducción al Utility Paper​

En esta lección ejecutamos un circuito a escala de utilidad que aparece en el informalmente denominado "Utility Paper", publicado en Nature, volumen 618, del 15 de junio de 2023. El artículo trata sobre la evolución temporal del modelo de Ising transversal bidimensional. Concretamente, los autores consideran la dinámica de espines del hamiltoniano:

H=HZZ+HX=−J∑(i,j)ZiZj+h∑iXiH = H_{ZZ} + H_X = - J \sum_{(i,j)} Z_i Z_j + h \sum_{i} X_i

donde J>0J > 0 es el acoplamiento de espines de vecinos más cercanos con i<ji < j y hh es el campo transversal global. La dinámica de espines se simula a partir de un estado inicial utilizando la descomposición de Trotter de primer orden del operador de evolución temporal:

exp⁡(−iHZZδt)=∏(i,j)exp⁡(iJδtZiZj)=∏(i,j)RZiZj(−2Jδt)exp⁡(−iHXδt)=∏iexp⁡(−ihδtXi)=∏iRXi(2hδt)\begin{aligned} \exp(-i H_{ZZ} \delta t) &= \prod_{(i,j)} \exp (i J \delta t Z_i Z_j) = \prod_{(i,j)} \mathrm{R}_{Z_i Z_j} ( - 2 J \delta t) \\ \exp(-i H_X \delta t) &= \prod_{i} \exp (-i h \delta t X_i ) = \prod_{i} \mathrm{R}_{X_i} ( 2 h \delta t) \end{aligned}

donde el tiempo de evolución TT se divide en T/δtT / \delta t pasos de Trotter y RZiZj(θJ)\mathrm{R}_{Z_i Z_j}(\theta_J) así como RXi(θh)\mathrm{R}_{X_i}(\theta_h) son puertas de rotación ZZZZ y XX respectivamente.

Los experimentos se realizaron en un procesador IBM Quantum® Eagle, un dispositivo de 127 qubits con conectividad heavy-hex, donde las interacciones XX se aplicaron a todos los qubits y las interacciones ZZZZ a todas las aristas del mapa de acoplamiento. Dado que no todas las interacciones ZZZZ pueden aplicarse simultáneamente debido a "dependencias de datos", el mapa de acoplamiento se colorea para agruparlas en capas. Las interacciones de una capa reciben el mismo color y pueden aplicarse en paralelo.

Además, los autores se centraron por simplicidad en el caso θJ=−π/2\theta_J=-\pi /2.

La contribución esencial del artículo consiste en crear circuitos cuánticos a una escala que va más allá de la simulación de vectores de estado, ejecutarlos en computadoras cuánticas ruidosas y aun así extraer resultados fiables. De este modo demostraron la utilidad de las computadoras cuánticas ruidosas. Para ello utilizaron la extrapolación a ruido cero (ZNE) con amplificación probabilística de errores (PEA) para la mitigación de errores.

Desde entonces, llamamos a tales experimentos y circuitos "a escala de utilidad".

1.1 Tu objetivo​

Tu objetivo en esta lección es crear un circuito a escala de utilidad y ejecutarlo en un procesador Eagle. Este notebook no trata de extraer resultados fiables, en parte porque PEA era aún una característica experimental de Qiskit en el momento de escribir este texto, y en parte porque la aplicación de ZNE con PEA requeriría un tiempo considerable.

Concretamente, debes crear y ejecutar el circuito correspondiente a la Figura 4b del artículo y representar tú mismo los puntos "sin mitigar". Como puedes ver, se trata de un circuito de 127 qubits ×\times 60 capas (20 pasos de Trotter) con ⟨Z62⟩\langle Z_{62} \rangle como observable. image.png ¿Parece una tarea grande?   No te preocupes. Las últimas tres lecciones de este curso proporcionan los pasos necesarios. Primero demostraremos un experimento más pequeño: crearemos un circuito de 27 qubits ×\times 6 capas (2 pasos de Trotter) con ⟨Z13⟩\langle Z_{13} \rangle como observable y lo ejecutaremos en un dispositivo simulado (fake device).

Esa fue la introducción. ¡Comencemos la aventura a escala de utilidad!

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx
import qiskit

qiskit.__version__
'2.0.2'
#!pip install qiskit_ibm_runtime
#!pip install qiskit_aer
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx

from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit.circuit.library import YGate
from qiskit.quantum_info import SparsePauliOp
from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    fake_provider,
    EstimatorV2 as Estimator,
)
from qiskit_aer import AerSimulator
service = QiskitRuntimeService()

2. Preparación​

2.1 Construir RZZ(-π\pi / 2)​

En primer lugar, observamos que la puerta RZZ en general requiere dos puertas CXCX.

from qiskit.circuit.library import RZZGate

θ_h = Parameter("$\\theta_h$")
qc1 = QuantumCircuit(2)
qc1.append(RZZGate(θ_h), [0, 1])
qc1.decompose(reps=1).draw("mpl")

Salida de la celda de código anterior

Como se menciónó anteriormente, para este experimento nos centramos en la puerta RZZ con el ángulo específico -π\pi / 2. Como se muestra en el artículo, se puede realizar con una sola puerta CXCX.

qc2 = QuantumCircuit(2)

qc2.sdg([0, 1])
qc2.append(YGate().power(1 / 2), [1])
qc2.cx(0, 1)
qc2.append(YGate().power(1 / 2).adjoint(), [1])

qc2.draw("mpl")

Salida de la celda de código anterior

Definimos una puerta a partir de este circuito para uso posterior.

rzz = qc2.to_gate(label="RZZ")

Hagamos un primer uso del rzz recién definido.

qc3 = QuantumCircuit(3)
qc3.append(rzz, [0, 1])
qc3.append(rzz, [0, 2])
display(qc3.draw("mpl"))
# display(qc.decompose(reps=1).draw("mpl"))

Salida de la celda de código anterior

Antes de seguir usando esta puerta, verifiquemos la equivalencia lógica de qc1 (la puerta RZZ) para -pi/2 y nuestra puerta recién definida rzz o qc2:

from qiskit.quantum_info import Operator

op1 = Operator(qc1.assign_parameters([-np.pi / 2]))
op2 = Operator(qc2)

op1.equiv(op2)
True

2.2 Colorear el mapa de acoplamiento​

Veamos cómo colorear el mapa de acoplamiento de un backend. Esto es necesario para agrupar las interacciones ZZZZ en capas.

Primero visualizamos el mapa de acoplamiento de un backend. Observa que los mapas de acoplamiento de todos los dispositivos IBM Quantum actuales son heavy-hexagonales.

backend = service.least_busy(operational=True, simulator=False)

backend.coupling_map.draw()

Salida de la celda de código anterior

Para colorear el mapa de acoplamiento usamos rustworkx, un paquete de Python para trabajar con grafos y redes complejas. Proporciona varios algoritmos de coloración, todos heurísticos y por lo tanto sin garantía de coloración mínima.

Dado que los grafos heavy-hex son bipartitos, usamos graph_bipartite_edge_color, que debería proporcionar una coloración mínima para estos grafos.

def color_coupling_map(backend):
    graph = backend.coupling_map.graph
    undirected_graph = graph.to_undirected(multigraph=False)
    edge_color_map = rx.graph_bipartite_edge_color(undirected_graph)
    if edge_color_map is None:
        edge_color_map = rx.graph_greedy_edge_color(undirected_graph)
    # build a map from color to a list of edges
    edge_index_map = undirected_graph.edge_index_map()
    color_edges_map = {color: [] for color in edge_color_map.values()}
    for edge_index, color in edge_color_map.items():
        color_edges_map[color].append(
            (edge_index_map[edge_index][0], edge_index_map[edge_index][1])
        )
    return edge_color_map, color_edges_map

Los grafos heavy-hexagonales deberían colorearse con tres colores. Verifiquémoslo para el mapa de acoplamiento anterior.

edge_color_map, color_edges_map = color_coupling_map(backend)
print(
    f"{backend.name}, {backend.num_qubits}-qubit device, {len(color_edges_map.keys())} colors assigned."
)
ibm_strasbourg, 127-qubit device, 3 colors assigned.

¡Exactamente correcto!

Por diversión, coloreemos el mapa de acoplamiento con la coloración obtenida, utilizando las funciones de visualización de rustworkx.

color_str_map = {0: "green", 1: "red", 2: "blue"}

undirected_graph = backend.coupling_map.graph.to_undirected(multigraph=False)
for i in undirected_graph.edge_indices():
    undirected_graph.get_edge_data_by_index(i)["color"] = color_str_map[
        edge_color_map[i]
    ]

rx.visualization.graphviz_draw(
    undirected_graph, method="neato", edge_attr_fn=lambda edge: {"color": edge["color"]}
)

Salida de la celda de código anterior

3. Resolver la evolución temporal trotterizada de un modelo de Ising 2D.​

Definimos una rutina para crear el circuito del Utility Paper para la evolución temporal de un modelo de Ising 2D. La rutina recibe tres parámetros: un backend, un entero para el número de pasos de Trotter y un booleano para controlar la inserción de barreras.

def get_utility_circuit(backend, num_steps: int, barrier: bool = False):
    num_qubits = backend.num_qubits
    _, color_edges_map = color_coupling_map(backend)
    θ_h = Parameter("$\\theta_h$")
    qc = QuantumCircuit(num_qubits)

    for i in range(num_steps):
        qc.rx(θ_h, range(num_qubits))

        for _, edge_list in color_edges_map.items():
            for edge in edge_list:
                qc.append(rzz, edge)

        if barrier:
            qc.barrier()
    return qc

Observa que ya hemos realizado manualmente el mapeo de qubits y el enrutamiento para el circuito creado. Cuando transpilemos el circuito más adelante, el transpilador no debe (no) volver a realizar el mapeo de qubits ni el enrutamiento. Como verás a continuación, lo llamamos con nivel de optimización 1 y el método de layout "trivial".

A continuación definimos una rutina sencilla para obtener información sobre el circuito creado, para una verificación rápida.

def get_circuit_info(qc: QuantumCircuit, reps: int = 0):
    qc0 = qc.decompose(reps=reps)
    return (
        f"{qc0.num_qubits} qubits × {qc0.depth(lambda x: x.operation.num_qubits == 2)} layers ({qc0.depth()}-depth)"
        + ",  "
        + f"""Gate breakdown: {", ".join([f"{k.upper()} {v}" for k, v in qc0.count_ops().items()])}"""
    )

Probemos estas rutinas. Deberías ver un circuito de 27 qubits ×\times 15 capas (5 pasos de Trotter). Dado que el dispositivo simulado tiene 28 aristas, debería haber 28*5 puertas de entrelazamiento.

backend = fake_provider.FakeTorontoV2()
num_steps = 5
qc = get_utility_circuit(backend, num_steps, True)

display(qc.draw(output="mpl", fold=-1))
print(get_circuit_info(qc, reps=0))
print(get_circuit_info(qc, reps=1))

Salida de la celda de código anterior

27 qubits × 15 layers (20-depth),  Gate breakdown: CIRCUIT-165 140, RX 135, BARRIER 5
27 qubits × 15 layers (60-depth),  Gate breakdown: SDG 280, UNITARY 280, CX 140, R 135, BARRIER 5

4. Resolver la versión de 27 qubits del problema.​

Ahora demostramos una versión más pequeña del experimento a escala de utilidad. Creamos un circuito de 27 qubits ×\times 6 capas (2 pasos de Trotter) con ⟨Z13⟩\langle Z_{13} \rangle como observable y lo ejecutamos tanto en el AerSimulator como en un dispositivo simulado (fake device).

Para ello seguimos, por supuesto, nuestro flujo de trabajo de cuatro pasos, el "patrón Qiskit", que consiste en Map, Optimize, Execute y Post-Process. Concretamente:

  • Mapear entradas clásicas a una computación cuántica.
  • Optimizar circuitos para la computación cuántica.
  • Ejecutar circuitos con Primitives.
  • Post-procesar resultados y devolverlos en formato clásico.

A continuación tenemos el paso Map para crear un circuito para el experimento más pequeño. Luego siguen un conjunto de Optimize y Execute para el AerSimulator y otro para el dispositivo simulado. Finalmente, tenemos el paso Post-Process para representar los resultados.

4.1 Paso 1: Map​

backend = fake_provider.FakeTorontoV2()  # a 27 qubit fake device.
num_steps = 2
qc = get_utility_circuit(backend, num_steps)
obs = SparsePauliOp.from_sparse_list(
    [("Z", [13], 1)], num_qubits=backend.num_qubits
)  # Falcon
angles = [
    0,
    0.1,
    0.2,
    0.3,
    0.4,
    0.5,
    0.6,
    0.7,
    0.8,
    1.0,
    np.pi / 2,
]  # We try 11 angles for theta_h.

4.2 Pasos 2 y 3: Optimizar y ejecutar (Simulador)​

backend_sim = AerSimulator()
transpiled_qc_sim = transpile(
    qc, backend_sim, optimization_level=1, layout_method="trivial"
)
transpiled_obs_sim = obs.apply_layout(layout=transpiled_qc_sim.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc_sim, reps=1))
27 qubits × 6 layers (23-depth),  Gate breakdown: SDG 112, UNITARY 112, CX 56, R 54
27 qubits × 6 layers (16-depth),  Gate breakdown: U3 80, CX 56, R 54, U1 32, U 28

Un usuario ejecutó la siguiente celda en un MacBook Pro con procesador Intel Core i7 Quad-Core de 2,3 GHz y 32 GB de RAM 3LPDDR4X bajo macOS 14.5. El tiempo de ejecución fue de 161 ms. Puede variar ligeramente según el portátil.

%%time
params = [[p] for p in angles]
estimator = Estimator(mode=backend_sim)
pub = (transpiled_qc_sim, transpiled_obs_sim, params)
result_sim = estimator.run([pub]).result()
CPU times: user 231 ms, sys: 186 ms, total: 417 ms
Wall time: 111 ms

4.3 Pasos 2 y 3: Optimizar y ejecutar (Fake Device)​

backend_fake = fake_provider.FakeTorontoV2()
transpiled_qc_fake = transpile(
    qc, backend_fake, optimization_level=1, layout_method="trivial"
)
transpiled_obs_fake = obs.apply_layout(layout=transpiled_qc_fake.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc_fake, reps=1))
27 qubits × 6 layers (23-depth),  Gate breakdown: SDG 112, UNITARY 112, CX 56, R 54
27 qubits × 6 layers (49-depth),  Gate breakdown: SDG 324, U1 274, H 162, CX 56, U3 14

Cuando el mismo usuario ejecutó la siguiente celda en el mismo entorno, el tiempo de ejecución fue de 2 min 19 s. La ejecución de un circuito en un fake device invoca una simulación ruidosa, que tarda considerablemente más que una simulación exacta. Recomendamos no ejecutar un circuito más grande (por ejemplo, 27 qubits ×\times 9 capas con 3 pasos de Trotter) en un fake device.

%%time
params = [[p] for p in angles]
estimator = Estimator(mode=backend_fake)
pub = (transpiled_qc_fake, transpiled_obs_fake, params)
result_fake = estimator.run([pub]).result()
CPU times: user 4min 42s, sys: 9.35 s, total: 4min 51s
Wall time: 38.3 s

4.4 Paso 4: Post-procesamiento​

Representamos los resultados de la simulación exacta y la ruidosa. Puedes ver los graves efectos del ruido en FakeToronto.

plt.plot(angles, result_fake[0].data.evs, "o", label="Fake Device")
plt.plot(angles, result_sim[0].data.evs, "o", label="AerSimulator")
plt.xlabel("$\\mathrm{R_x}$ angle $\\theta_h$")
plt.title("$\\langle Z_{13} \\rangle$")
plt.legend()
plt.show()

Salida de la celda de código anterior

5. Resolver la versión de 127 qubits del problema​

Tu objetivo es realizar el experimento a escala de utilidad como se describió al principio. Crearás y ejecutarás un circuito de 127 qubits ×\times 60 capas (20 pasos de Trotter) con ⟨Z62⟩\langle Z_{62} \rangle como observable. Te recomendamos que lo intentes por tu cuenta usando el código de la versión de 27 qubits como plantilla. Sin embargo, la solución se proporciona aquí.

Solución:

5.1 Paso 1: Map​

# backend_map = service.backend("ibm_brisbane")
backend_map = service.least_busy(operational=True, simulator=False)

num_steps = 20
qc = get_utility_circuit(backend_map, num_steps)
obs = SparsePauliOp.from_sparse_list(
    [("Z", [62], 1)], num_qubits=backend_map.num_qubits
)  # Eagle
angles = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, np.pi / 2]

5.2 Pasos 2 y 3: Optimizar y ejecutar​

Nota: El mapa de acoplamiento del procesador Eagle tiene 144 aristas.

# backend = service.backend("ibm_brisbane")
backend = backend_map

transpiled_qc = transpile(qc, backend, optimization_level=1, layout_method="trivial")
transpiled_obs = obs.apply_layout(layout=transpiled_qc.layout)

print(get_circuit_info(qc, reps=1))
print(get_circuit_info(transpiled_qc))
156 qubits × 60 layers (221-depth),  Gate breakdown: SDG 7040, UNITARY 7040, CX 3520, R 3120
156 qubits × 60 layers (201-depth),  Gate breakdown: RZ 11933, SX 6240, CZ 3520
params = [[p] for p in angles]
estimator = Estimator(mode=backend)
pub = (transpiled_qc, transpiled_obs, params)
job = estimator.run([pub])

job_id = job.job_id()
print(f"job id={job_id}")
job id=d1479n6qf56g0081sxa0

5.3 Post-procesamiento​

Proporcionamos los valores para los puntos "mitigados" de la Figura 4b del Utility Paper. Represéntalos junto con tus propios resultados.

result_paper = [
    1.0171,
    1.0044,
    0.9563,
    0.9602,
    0.8394,
    0.8120,
    0.5466,
    0.4556,
    0.1953,
    0.0141,
    0.0117,
]

# REPLACE WITH YOUR OWN JOB ID
job = service.job(job_id)

plt.plot(angles, job.result()[0].data.evs, "o", label=f"{job.backend().name}")
plt.plot(angles, result_paper, "o", label="Utility Paper")
plt.xlabel("$\\mathrm{R_x}$ angle $\\theta_h$")
plt.title("$\\langle Z_{62} \\rangle$")
plt.legend()
plt.show()

Salida de la celda de código anterior

¿Se parecen tus resultados a los puntos "sin mitigar" de la Figura 4b?   Podrían ser muy diferentes, dependiendo del dispositivo y su estado en el momento del experimento. No te preocupes por los resultados en sí. Lo que verificamos es si has realizado la programación correctamente. Si es así, ¡enhorabuena! Has llegado a la línea de salida de la era de utilidad.

Como se describe en el Utility Paper, científicos y científicas de todo el mundo han invertido gran ingenio para extraer resultados significativos incluso en presencia de ruido. El objetivo general de estos esfuerzos colectivos es la ventaja cuántica: un estado en el que las computadoras cuánticas pueden resolver ciertos problemas relevantes para la industria más rápido, con mayor precisión o de forma más económica que las computadoras clásicas. Probablemente no se tratará de un evento único, sino de una era en la que la reproducción clásica de resultados cuánticos tardará cada vez más, hasta que esa ventaja se vuelva decisivamente importante. Una cosa está clara: solo alcanzaremos la ventaja cuántica a través de experimentos a escala de utilidad. Si este curso te lleva a unirte al camino hacia ese objetivo, lleno de desafíos y diversión, nos alegraríamos mucho.

Referencia​