Simulación del hamiltoniano de Ising con patada usando circuitos dinámicos

Estimación de uso: 7,5 minutos en un procesador Heron r3. (NOTA: Esto es solo una estimación. Su tiempo de ejecución puede variar.) Los circuitos dinámicos son circuitos con retroalimentación clásica; en otras palabras, son mediciones a mitad de circuito seguidas de operaciones lógicas clásicas que determinan operaciones cuánticas condicionadas por la salida clásica. En este tutorial, simulamos el modelo de Ising con patada en una red hexagonal de espines y usamos circuitos dinámicos para realizar interacciones más allá de la conectividad física del hardware.

El modelo de Ising ha sido estudiado extensamente en diversas áreas de la física. Modela espines que experimentan interacciones de Ising entre sitios de la red, así como patadas del campo magnético local en cada sitio. La evolución temporal trotterizada de los espines considerada en este tutorial, tomada de [1], está dada por el siguiente unitario:

$$U(\theta)=\left(\prod_{\langle j, k\rangle} \exp \left(i \frac{\pi}{8} Z_j Z_k\right)\right)\left(\prod_j \exp \left(-i \frac{\theta}{2} X_j\right)\right)$$

Para sondear la dinámica de los espines, estudiamos la magnetización promedio de los espines en cada sitio como función de los pasos de Trotter. Por lo tanto, construimos el siguiente observable:

$$\langle O\rangle = \frac{1}{N} \sum_i \langle Z_i \rangle$$

Para realizar la interacción ZZ entre sitios de la red, presentamos una solución usando la funcionalidad de circuitos dinámicos, lo que conduce a una profundidad de puertas de dos qubits significativamente menor en comparación con el método estándar de enrutamiento con puertas SWAP. Por otro lado, las operaciones de retroalimentación clásica en circuitos dinámicos típicamente tienen tiempos de ejecución más largos que las puertas cuánticas; por lo tanto, los circuitos dinámicos tienen limitaciones y compromisos. También presentamos una forma de agregar una secuencia de desacoplamiento dinámico en los qubits inactivos durante la operación de retroalimentación clásica usando la duración stretch.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener instalado lo siguiente:

• Qiskit SDK v2.0 o posterior con soporte de visualización
• Qiskit Runtime v0.37 o posterior con soporte de visualización (pip install 'qiskit-ibm-runtime[visualization]')
• Biblioteca de grafos Rustworkx (pip install rustworkx)
• Qiskit Aer (pip install qiskit-aer)

Configuración

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

import numpy as np
from typing import List
import rustworkx as rx
import matplotlib.pyplot as plt
from rustworkx.visualization import mpl_draw
from qiskit.circuit import (
    Parameter,
    QuantumCircuit,
    QuantumRegister,
    ClassicalRegister,
)
from qiskit.transpiler import CouplingMap
from qiskit.quantum_info import SparsePauliOp
from qiskit.circuit.classical import expr
from qiskit.transpiler.preset_passmanagers import (
    generate_preset_pass_manager,
)
from qiskit.transpiler import PassManager
from qiskit.circuit.library import RZGate, XGate
from qiskit.transpiler.passes import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)

from qiskit.transpiler.basepasses import TransformationPass
from qiskit.circuit.measure import Measure
from qiskit.transpiler.passes.utils.remove_final_measurements import (
    calc_final_ops,
)
from qiskit.circuit import Instruction

from qiskit.visualization import plot_circuit_layout
from qiskit.circuit.tools import pi_check

from qiskit_aer import AerSimulator
from qiskit_aer.primitives import SamplerV2 as Aer_Sampler

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Batch,
    SamplerV2 as Sampler,
)
from qiskit_ibm_runtime.exceptions import QiskitBackendNotFoundError
from qiskit_ibm_runtime.visualization import (
    draw_circuit_schedule_timing,
)

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

Comenzamos definiendo la red a simular. Elegimos trabajar con la red de panal de abeja (también llamada hexagonal), que es un grafo planar con nodos de grado 3. Aquí, especificamos el tamaño de la red, los parámetros de circuito relevantes de interés en la dinámica trotterizada. Simulamos la evolución temporal trotterizada bajo el modelo de Ising con tres valores diferentes de $\theta$ del campo magnético local.

hex_rows = 3  # specify lattice size
hex_cols = 5
depths = range(9)  # specify Trotter steps
zz_angle = np.pi / 8  # parameter for ZZ interaction
max_angle = np.pi / 2  # max theta angle
points = 3  # number of theta parameters

θ = Parameter("θ")
params = np.linspace(0, max_angle, points)

def make_hex_lattice(hex_rows=1, hex_cols=1):
    """Define hexagon lattice."""
    hex_cmap = CouplingMap.from_hexagonal_lattice(
        hex_rows, hex_cols, bidirectional=False
    )
    data = list(hex_cmap.physical_qubits)
    graph = hex_cmap.graph.to_undirected(multigraph=False)
    edge_colors = rx.graph_misra_gries_edge_color(graph)
    layer_edges = {color: [] for color in edge_colors.values()}
    for edge_index, color in edge_colors.items():
        layer_edges[color].append(graph.edge_list()[edge_index])
    return data, layer_edges, hex_cmap, graph

Comencemos con un pequeño ejemplo de prueba:

hex_rows_test = 1
hex_cols_test = 2

data_test, layer_edges_test, hex_cmap_test, graph_test = make_hex_lattice(
    hex_rows=hex_rows_test, hex_cols=hex_cols_test
)

# display a small example for illustration
node_colors_test = ["lightblue"] * len(graph_test.node_indices())
pos = rx.graph_spring_layout(
    graph_test,
    k=5 / np.sqrt(len(graph_test.nodes())),
    repulsive_exponent=1,
    num_iter=150,
)
mpl_draw(graph_test, node_color=node_colors_test, pos=pos)

Usaremos el ejemplo pequeño para ilustración y simulación. A continuación, también construimos un ejemplo grande para mostrar que el flujo de trabajo puede extenderse a tamaños mayores.

data, layer_edges, hex_cmap, graph = make_hex_lattice(
    hex_rows=hex_rows, hex_cols=hex_cols
)
num_qubits = len(data)
print(f"num_qubits = {num_qubits}")

# display the honeycomb lattice to simulate
node_colors = ["lightblue"] * len(graph.node_indices())
pos = rx.graph_spring_layout(
    graph,
    k=5 / np.sqrt(num_qubits),
    repulsive_exponent=1,
    num_iter=150,
)
mpl_draw(graph, node_color=node_colors, pos=pos)
plt.show()

num_qubits = 46

Construir circuitos unitarios

Con el tamaño del problema y los parámetros especificados, estamos listos para construir el circuito parametrizado que simula la evolución temporal trotterizada de $U(\theta)$ con diferentes pasos de Trotter, especificados por el argumento depth. El circuito que construimos tiene capas alternadas de puertas Rx($\theta$) y puertas Rzz. Las puertas Rzz realizan las interacciones ZZ entre espines acoplados, que se colocarán entre cada sitio de la red especificado por el argumento layer_edges.

def gen_hex_unitary(
    num_qubits=6,
    zz_angle=np.pi / 8,
    layer_edges=[
        [(0, 1), (2, 3), (4, 5)],
        [(1, 2), (3, 4), (5, 0)],
    ],
    θ=Parameter("θ"),
    depth=1,
    measure=False,
    final_rot=True,
):
    """Build unitary circuit."""
    circuit = QuantumCircuit(num_qubits)
    # Build trotter layers
    for _ in range(depth):
        for i in range(num_qubits):
            circuit.rx(θ, i)
        circuit.barrier()
        for coloring in layer_edges.keys():
            for e in layer_edges[coloring]:
                circuit.rzz(zz_angle, e[0], e[1])
        circuit.barrier()
    # Optional final rotation, set True to be consistent with Ref. [1]
    if final_rot:
        for i in range(num_qubits):
            circuit.rx(θ, i)
    if measure:
        circuit.measure_all()

    return circuit

Visualiza el circuito de prueba pequeño:

circ_unitary_test = gen_hex_unitary(
    num_qubits=len(data_test),
    layer_edges=layer_edges_test,
    θ=Parameter("θ"),
    depth=1,
    measure=True,
)
circ_unitary_test.draw(output="mpl", fold=-1)

De manera similar, construye los circuitos unitarios del ejemplo grande en diferentes pasos de Trotter y el observable para estimar el valor esperado.

circuits_unitary = []
for depth in depths:
    circ = gen_hex_unitary(
        num_qubits=num_qubits,
        layer_edges=layer_edges,
        θ=Parameter("θ"),
        depth=depth,
        measure=True,
    )
    circuits_unitary.append(circ)

observables_unitary = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
    num_qubits=num_qubits,
)

Construir la implementación con circuitos dinámicos

Esta sección demuestra la implementación principal con circuitos dinámicos para simular la misma evolución temporal trotterizada. Observa que la red de panal de abeja que queremos simular no coincide con la red heavy-hex de los qubits del hardware. Una forma directa de mapear el circuito al hardware es introducir una serie de operaciones SWAP para acercar los qubits que interactúan entre sí, para realizar la interacción ZZ. Aquí destacamos un enfoque alternativo usando circuitos dinámicos como solución, que ilustra que podemos usar la combinación de computación cuántica y clásica en tiempo real dentro de un circuito en Qiskit para realizar interacciones más allá de vecinos más cercanos.

En la implementación con circuitos dinámicos, la interacción ZZ se implementa efectivamente usando qubits ancilla, mediciones a mitad de circuito y retroalimentación. Para entender esto, observa que las rotaciones ZZ aplican un factor de fase $e^{i\theta}$ al estado basándose en su paridad. Para dos qubits, los estados de la base computacional son $|00\rangle$, $|01\rangle$, $|10\rangle$ y $|11\rangle$. La puerta de rotación ZZ aplica un factor de fase a los estados $|01\rangle$ y $|10\rangle$ cuya paridad (el número de unos en el estado) es impar y deja los estados de paridad par sin cambios. A continuación se describe cómo podemos implementar efectivamente interacciones ZZ en dos qubits usando circuitos dinámicos.

1. Calcular la paridad en un qubit ancilla: en lugar de aplicar directamente ZZ a dos qubits, introducimos un tercer qubit, el qubit ancilla, para almacenar la información de paridad de los dos qubits de datos. Entrelazamos el ancilla con cada qubit de datos usando puertas CX desde el qubit de datos hacia el qubit ancilla.

2. Aplicar una rotación Z de un solo qubit al qubit ancilla: esto se debe a que el ancilla tiene la información de paridad de los dos qubits de datos, lo que efectivamente implementa la rotación ZZ en los qubits de datos.

3. Medir el qubit ancilla en la base X: este es el paso clave que colapsa el estado del qubit ancilla, y el resultado de la medición nos indica lo que ha ocurrido:

   • Medir 0: cuando se observa un resultado 0, de hecho hemos aplicado correctamente una rotación $ZZ(\theta)$ a nuestros qubits de datos.

   • Medir 1: cuando se observa un resultado 1, hemos aplicado $ZZ(\theta + \pi)$ en su lugar.

4. Aplicar puerta de corrección al medir 1: Si medimos 1, aplicamos puertas Z a los qubits de datos para "corregir" la fase extra de $\pi$.

El circuito resultante es el siguiente:

Cuando adoptamos este enfoque para simular una red de panal de abeja, el circuito resultante se incrusta perfectamente en el hardware con una red heavy-hex: todos los qubits de datos residen en los sitios de grado 3 de la red, que forma una red hexagonal. Cada par de qubits de datos comparte un qubit ancilla que reside en un sitio de grado 2. A continuación, construimos la red de qubits para la implementación con circuitos dinámicos, introduciendo qubits ancilla (mostrados en los círculos morados más oscuros).

def make_lattice(hex_rows=1, hex_cols=1):
    """Define heavy-hex lattice and corresponding lists of data and ancilla nodes."""
    hex_cmap = CouplingMap.from_hexagonal_lattice(
        hex_rows, hex_cols, bidirectional=False
    )
    data = list(hex_cmap.physical_qubits)

    heavyhex_cmap = CouplingMap()
    for d in data:
        heavyhex_cmap.add_physical_qubit(d)

    # make coupling map
    a = len(data)
    for edge in hex_cmap.get_edges():
        heavyhex_cmap.add_physical_qubit(a)
        heavyhex_cmap.add_edge(edge[0], a)
        heavyhex_cmap.add_edge(edge[1], a)
        a += 1
    ancilla = list(range(len(data), a))
    qubits = data + ancilla

    # color edges
    graph = heavyhex_cmap.graph.to_undirected(multigraph=False)
    edge_colors = rx.graph_misra_gries_edge_color(graph)
    layer_edges = {color: [] for color in edge_colors.values()}
    for edge_index, color in edge_colors.items():
        layer_edges[color].append(graph.edge_list()[edge_index])

    # construct observable
    obs_hex = SparsePauliOp.from_sparse_list(
        [("Z", [i], 1 / len(data)) for i in data],
        num_qubits=len(qubits),
    )

    return (data, qubits, ancilla, layer_edges, heavyhex_cmap, graph, obs_hex)

Visualiza la red heavy-hex para qubits de datos y qubits ancilla a pequeña escala:

(data, qubits, ancilla, layer_edges, heavyhex_cmap, graph, obs_hex) = (
    make_lattice(hex_rows=hex_rows, hex_cols=hex_cols)
)

print(f"number of data qubits = {len(data)}")
print(f"number of ancilla qubits = {len(ancilla)}")

node_colors = []
for node in graph.node_indices():
    if node in ancilla:
        node_colors.append("purple")
    else:
        node_colors.append("lightblue")

pos = rx.graph_spring_layout(
    graph,
    k=1 / np.sqrt(len(qubits)),
    repulsive_exponent=2,
    num_iter=200,
)

# Visualize the graph, blue circles are data qubits and purple circles are ancillas
mpl_draw(graph, node_color=node_colors, pos=pos)
plt.show()

number of data qubits = 46
number of ancilla qubits = 60

A continuación, construimos el circuito dinámico para la evolución temporal trotterizada. Las puertas RZZ se reemplazan con la implementación de circuitos dinámicos usando los pasos descritos anteriormente.

def gen_hex_dynamic(
    depth=1,
    zz_angle=np.pi / 8,
    θ=Parameter("θ"),
    hex_rows=1,
    hex_cols=1,
    measure=False,
    add_dd=True,
):
    """Build dynamic circuits."""
    (data, qubits, ancilla, layer_edges, heavyhex_cmap, graph, obs_hex) = (
        make_lattice(hex_rows=hex_rows, hex_cols=hex_cols)
    )
    # Initialize circuit
    qr = QuantumRegister(len(qubits), "qr")
    cr = ClassicalRegister(len(ancilla), "cr")
    circuit = QuantumCircuit(qr, cr)

    for k in range(depth):
        # Single-qubit Rx layer
        for d in data:
            circuit.rx(θ, d)
        circuit.barrier()

        # CX gates from data qubits to ancilla qubits
        for same_color_edges in layer_edges.values():
            for e in same_color_edges:
                circuit.cx(e[0], e[1])
        circuit.barrier()

        # Apply Rz rotation on ancilla qubits and rotate into X basis
        for a in ancilla:
            circuit.rz(zz_angle, a)
            circuit.h(a)
        # Add barrier to align terminal measurement
        circuit.barrier()

        # Measure ancilla qubits
        for i, a in enumerate(ancilla):
            circuit.measure(a, i)
        d2ros = {}
        a2ro = {}
        # Retrieve ancilla measurement outcomes
        for a in ancilla:
            a2ro[a] = cr[ancilla.index(a)]

        # For each data qubit, retrieve measurement outcomes of neighboring ancilla qubits
        for d in data:
            ros = [a2ro[a] for a in heavyhex_cmap.neighbors(d)]
            d2ros[d] = ros

        # Build classical feedforward operations (optionally add DD on idling data qubits)
        for d in data:
            if add_dd:
                circuit = add_stretch_dd(circuit, d, f"data_{d}_depth_{k}")

            # # XOR the neighboring readouts of the data qubit; if True, apply Z to it
            ros = d2ros[d]
            parity = ros[0]
            for ro in ros[1:]:
                parity = expr.bit_xor(parity, ro)
            with circuit.if_test(expr.equal(parity, True)):
                circuit.z(d)

        # Reset the ancilla if its readout is 1
        for a in ancilla:
            with circuit.if_test(expr.equal(a2ro[a], True)):
                circuit.x(a)
        circuit.barrier()

    # Final single-qubit Rx layer to match the unitary circuits
    for d in data:
        circuit.rx(θ, d)

    if measure:
        circuit.measure_all()
    return circuit, obs_hex

def add_stretch_dd(qc, q, name):
    """Add XpXm DD sequence."""
    s = qc.add_stretch(name)
    qc.delay(s, q)
    qc.x(q)
    qc.delay(s, q)
    qc.delay(s, q)
    qc.rz(np.pi, q)
    qc.x(q)
    qc.rz(-np.pi, q)
    qc.delay(s, q)
    return qc

Desacoplamiento dinámico (DD) y soporte para la duración stretch

Una advertencia al usar la implementación con circuitos dinámicos para realizar la interacción ZZ es que las mediciones a mitad de circuito y las operaciones de retroalimentación clásica típicamente toman más tiempo en ejecutarse que las puertas cuánticas. Para suprimir la decoherencia de los qubits durante el tiempo de inactividad mientras se ejecutan las operaciones clásicas, agregamos una secuencia de desacoplamiento dinámico (DD) después de la operación de medición en los qubits ancilla, y antes de la operación Z condicional en el qubit de datos, antes de la sentencia if_test.

La secuencia DD se agrega mediante la función add_stretch_dd(), que utiliza las duraciones stretch para determinar los intervalos de tiempo entre las puertas DD. Una duración stretch es una forma de especificar una duración de tiempo extensible para la operación delay de modo que la duración del retraso pueda crecer para llenar el tiempo de inactividad del qubit. Las variables de duración especificadas por stretch se resuelven en tiempo de compilación en duraciones deseadas que satisfacen una cierta restricción. Esto es muy útil cuando la temporización de las secuencias DD es esencial para lograr un buen rendimiento de supresión de errores. Para más detalles sobre el tipo stretch, consulta la documentación de OpenQASM. Actualmente, el soporte para el tipo stretch en Qiskit Runtime es experimental. Para detalles sobre sus restricciones de uso, consulta la sección de limitaciones de la documentación de stretch.

Usando las funciones definidas anteriormente, construimos los circuitos de evolución temporal trotterizada, con y sin DD, y los observables correspondientes. Comenzamos visualizando el circuito dinámico de un ejemplo pequeño:

hex_rows_test = 1
hex_cols_test = 1

(
    data_test,
    qubits_test,
    ancilla_test,
    layer_edges_test,
    heavyhex_cmap_test,
    graph_test,
    obs_hex_test,
) = make_lattice(hex_rows=hex_rows_test, hex_cols=hex_cols_test)

node_colors = []
for node in graph_test.node_indices():
    if node in ancilla_test:
        node_colors.append("purple")
    else:
        node_colors.append("lightblue")
pos = rx.graph_spring_layout(
    graph_test,
    k=5 / np.sqrt(len(qubits_test)),
    repulsive_exponent=2,
    num_iter=150,
)

# display a small example for illustration
node_colors_test = ["lightblue"] * len(graph_test.node_indices())
mpl_draw(graph_test, node_color=node_colors, pos=pos)

circuit_dynamic_test, obs_dynamic_test = gen_hex_dynamic(
    depth=1,
    θ=Parameter("θ"),
    hex_rows=hex_rows_test,
    hex_cols=hex_cols_test,
    measure=False,
    add_dd=False,
)
circuit_dynamic_test.draw("mpl", fold=-1)

circuit_dynamic_dd_test, _ = gen_hex_dynamic(
    depth=1,
    θ=Parameter("θ"),
    hex_rows=hex_rows_test,
    hex_cols=hex_cols_test,
    measure=False,
    add_dd=True,
)
circuit_dynamic_dd_test.draw("mpl", fold=-1)

De manera similar, construye los circuitos dinámicos para el ejemplo grande:

circuits_dynamic = []
circuits_dynamic_dd = []
observables_dynamic = []
for depth in depths:
    circuit, obs = gen_hex_dynamic(
        depth=depth,
        θ=Parameter("θ"),
        hex_rows=hex_rows,
        hex_cols=hex_cols,
        measure=True,
        add_dd=False,
    )
    circuits_dynamic.append(circuit)

    circuit_dd, _ = gen_hex_dynamic(
        depth=depth,
        θ=Parameter("θ"),
        hex_rows=hex_rows,
        hex_cols=hex_cols,
        measure=True,
        add_dd=True,
    )
    circuits_dynamic_dd.append(circuit_dd)
    observables_dynamic.append(obs)

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

Ahora estamos listos para transpilar el circuito al hardware. Transpilaremos tanto la implementación unitaria estándar como la implementación con circuitos dinámicos al hardware.

Para transpilar al hardware, primero instanciamos el backend. Si está disponible, elegiremos un backend donde la instrucción MidCircuitMeasure (measure_2) sea compatible.

service = QiskitRuntimeService()
try:
    backend = service.least_busy(
        operational=True,
        simulator=False,
        use_fractional_gates=True,
        filters=lambda b: "measure_2" in b.supported_instructions,
    )
except QiskitBackendNotFoundError:
    backend = service.least_busy(
        operational=True,
        simulator=False,
        use_fractional_gates=True,
    )

Transpilación para circuitos dinámicos

Primero, transpilamos los circuitos dinámicos, con y sin agregar la secuencia DD. Para asegurar que usamos el mismo conjunto de qubits físicos en todos los circuitos y obtener resultados más consistentes, primero transpilamos el circuito una vez, y luego usamos su disposición para todos los circuitos subsiguientes, especificada por initial_layout en el administrador de pases. Luego construimos los bloques unificados de primitivas (PUBs) como entrada de la primitiva Sampler.

pm_temp = generate_preset_pass_manager(
    optimization_level=3,
    backend=backend,
)
isa_temp = pm_temp.run(circuits_dynamic[-1])
dynamic_layout = isa_temp.layout.initial_index_layout(filter_ancillas=True)

pm = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=dynamic_layout
)

dynamic_isa_circuits = [pm.run(circ) for circ in circuits_dynamic]
dynamic_pubs = [(circ, params) for circ in dynamic_isa_circuits]

dynamic_isa_circuits_dd = [pm.run(circ) for circ in circuits_dynamic_dd]
dynamic_pubs_dd = [(circ, params) for circ in dynamic_isa_circuits_dd]

Podemos visualizar la disposición de qubits del circuito transpilado a continuación. Los círculos negros muestran los qubits de datos y los qubits ancilla usados en la implementación con circuitos dinámicos.

def _heron_coords_r2():
    cord_map = np.array(
        [
            [
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                1,
                5,
                9,
                13,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                1,
                5,
                9,
                13,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                1,
                5,
                9,
                13,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
                3,
                7,
                11,
                15,
                0,
                1,
                2,
                3,
                4,
                5,
                6,
                7,
                8,
                9,
                10,
                11,
                12,
                13,
                14,
                15,
            ],
            -1
            * np.array([j for i in range(15) for j in [i] * [16, 4][i % 2]]),
        ],
        dtype=int,
    )

    hcords = []
    ycords = cord_map[0]
    xcords = cord_map[1]
    for i in range(156):
        hcords.append([xcords[i] + 1, np.abs(ycords[i]) + 1])

    return hcords

plot_circuit_layout(
    dynamic_isa_circuits_dd[8],
    backend,
    qubit_coordinates=_heron_coords_r2(),
    view="virtual",
)

nota

Si obtiene errores sobre que neato no se encuentra desde plot_circuit_layout(), asegúrate de tener el paquete graphviz instalado y disponible en su PATH. Si se instala en una ubicación no predeterminada (por ejemplo, usando homebrew en MacOS), es posible que necesite actualizar su variable de entorno PATH. Esto se puede hacer dentro de este notebook usando lo siguiente:

import os
os.environ['PATH'] = f"path/to/neato{os.pathsep}{os.environ['PATH']}"

dynamic_isa_circuits[1].draw(fold=-1, output="mpl", idle_wires=False)

dynamic_isa_circuits_dd[1].draw(fold=-1, output="mpl", idle_wires=False)

Transpilar usando MidCircuitMeasure

MidCircuitMeasure es una adición a las operaciones de medición disponibles, calibrada específicamente para realizar mediciones a mitad de circuito. La instrucción MidCircuitMeasure se mapea a la instrucción measure_2 soportada por los backends. Ten en cuenta que measure_2 no es compatible con todos los backends. Puede usar service.backends(filters=lambda b: "measure_2" in b.supported_instructions) para encontrar backends que la soporten. Aquí, mostramos cómo transpilar el circuito para que las mediciones a mitad de circuito definidas en el circuito se ejecuten usando la operación MidCircuitMeasure, si el backend la soporta.

A continuación, imprimimos la duración de la instrucción measure_2 y la instrucción measure estándar.

print(
    f'Mid-circuit measurement `measure_2` duration: {backend.instruction_durations.get('measure_2',0) * backend.dt * 1e9/1e3} μs'
)
print(
    f'Terminal measurement `measure` duration: {backend.instruction_durations.get('measure',0) * backend.dt *1e9/1e3} μs'
)

Mid-circuit measurement `measure_2` duration:  1.624 μs
Terminal measurement `measure` duration: 2.2 μs

"""Pass that replaces terminal measures in the middle of the circuit with
MidCircuitMeasure instructions."""

class ConvertToMidCircuitMeasure(TransformationPass):
    """This pass replaces terminal measures in the middle of the circuit with
    MidCircuitMeasure instructions.
    """

    def __init__(self, target):
        super().__init__()
        self.target = target

    def run(self, dag):
        """Run the pass on a dag."""
        mid_circ_measure = None
        for inst in self.target.instructions:
            if isinstance(inst[0], Instruction) and inst[0].name.startswith(
                "measure_"
            ):
                mid_circ_measure = inst[0]
                break
        if not mid_circ_measure:
            return dag

        final_measure_nodes = calc_final_ops(dag, {"measure"})
        for node in dag.op_nodes(Measure):
            if node not in final_measure_nodes:
                dag.substitute_node(node, mid_circ_measure, inplace=True)

        return dag

pm = PassManager(ConvertToMidCircuitMeasure(backend.target))

dynamic_isa_circuits_meas2 = [pm.run(circ) for circ in dynamic_isa_circuits]
dynamic_pubs_meas2 = [(circ, params) for circ in dynamic_isa_circuits_meas2]

dynamic_isa_circuits_dd_meas2 = [
    pm.run(circ) for circ in dynamic_isa_circuits_dd
]
dynamic_pubs_dd_meas2 = [
    (circ, params) for circ in dynamic_isa_circuits_dd_meas2
]

Transpilación para circuitos unitarios

Para establecer una comparación justa entre los circuitos dinámicos y su contraparte unitaria, usamos el mismo conjunto de qubits físicos utilizados en los circuitos dinámicos para los qubits de datos como la disposición para transpilar los circuitos unitarios.

init_layout = [
    dynamic_layout[ind] for ind in range(circuits_unitary[0].num_qubits)
]

pm = generate_preset_pass_manager(
    target=backend.target,
    initial_layout=init_layout,
    optimization_level=3,
)

def transpile_minimize(circ: QuantumCircuit, pm: PassManager, iterations=10):
    """Transpile circuits for specified number of iterations and return the one with smallest two-qubit gate depth"""
    circs = [pm.run(circ) for i in range(iterations)]
    circs_sorted = sorted(
        circs,
        key=lambda x: x.depth(lambda x: x.operation.num_qubits == 2),
    )
    return circs_sorted[0]

unitary_isa_circuits = []
for circ in circuits_unitary:
    circ_t = transpile_minimize(circ, pm, iterations=100)
    unitary_isa_circuits.append(circ_t)

unitary_pubs = [(circ, params) for circ in unitary_isa_circuits]

Visualizamos la disposición de qubits de los circuitos unitarios transpilados. Los círculos negros indican los qubits físicos usados para transpilar los circuitos unitarios y sus índices corresponden a los índices de qubits virtuales. Al comparar esto con la disposición graficada para los circuitos dinámicos, podemos confirmar que los circuitos unitarios usan el mismo conjunto de qubits físicos que los qubits de datos en los circuitos dinámicos.

plot_circuit_layout(
    unitary_isa_circuits[-1],
    backend,
    qubit_coordinates=_heron_coords_r2(),
    view="virtual",
)

Ahora agregamos la secuencia DD a los circuitos transpilados y construimos los PUBs correspondientes para el envío de trabajos.

pm_dd = PassManager(
    [
        ALAPScheduleAnalysis(target=backend.target),
        PadDynamicalDecoupling(
            dd_sequence=[
                XGate(),
                RZGate(np.pi),
                XGate(),
                RZGate(-np.pi),
            ],
            spacing=[1 / 4, 1 / 2, 0, 0, 1 / 4],
            target=backend.target,
        ),
    ]
)

unitary_isa_circuits_dd = pm_dd.run(unitary_isa_circuits)
unitary_pubs_dd = [(circ, params) for circ in unitary_isa_circuits_dd]

Comparar la profundidad de puertas de dos qubits de los circuitos unitarios y dinámicos

# compare circuit depth of unitary and dynamic circuit implementations
unitary_depth = [
    unitary_isa_circuits[i].depth(lambda x: x.operation.num_qubits == 2)
    for i in range(len(unitary_isa_circuits))
]

dynamic_depth = [
    dynamic_isa_circuits[i].depth(lambda x: x.operation.num_qubits == 2)
    for i in range(len(dynamic_isa_circuits))
]

plt.plot(
    list(range(len(unitary_depth))),
    unitary_depth,
    label="unitary circuits",
    color="#be95ff",
)
plt.plot(
    list(range(len(dynamic_depth))),
    dynamic_depth,
    label="dynamic circuits",
    color="#ff7eb6",
)
plt.xlabel("Trotter steps")
plt.ylabel("Two-qubit depth")
plt.legend()

<matplotlib.legend.Legend at 0x374225760>

El principal beneficio del circuito basado en mediciones es que, al implementar múltiples interacciones ZZ, las capas CX pueden paralelizarse, y las mediciones pueden ocurrir simultáneamente. Esto se debe a que todas las interacciones ZZ conmutan, por lo que el cálculo puede realizarse con profundidad de medición 1. Después de transpilar los circuitos, observamos que el enfoque de circuitos dinámicos produce una profundidad de puertas de dos qubits significativamente menor que el enfoque unitario estándar, con la advertencia de que las mediciones a mitad de circuito adicionales y la retroalimentación clásica en sí mismas toman tiempo e introducen sus propias fuentes de errores.

Paso 3: Ejecutar usando primitivas de Qiskit

Modo de prueba local

Antes de enviar los trabajos al hardware, podemos ejecutar una pequeña simulación de prueba del circuito dinámico usando el modo de prueba local.

aer_sim = AerSimulator()
pm = generate_preset_pass_manager(backend=aer_sim, optimization_level=1)
circuit_dynamic_test.measure_all()
isa_qc = pm.run(circuit_dynamic_test)
with Batch(backend=aer_sim) as batch:
    sampler = Sampler(mode=batch)
    result = sampler.run([(isa_qc, params)]).result()

print(
    "Simulated average magnetization at trotter step = 1 at three theta values"
)
result[0].data["meas"].expectation_values(obs_dynamic_test[0])

Simulated average magnetization at trotter step = 1 at three theta values

array([ 0.16666667,  0.01855469, -0.13476562])

Simulación MPS

Para circuitos grandes, podemos usar el simulador matrix_product_state (MPS), que proporciona un resultado aproximado del valor esperado de acuerdo con la dimensión de enlaza elegida. Más adelante usamos los resultados de la simulación MPS como línea base para comparar los resultados del hardware.

# The MPS simulation below took approximately 7 minutes to run on a laptop with Apple M1 chip

mps_backend = AerSimulator(
    method="matrix_product_state",
    matrix_product_state_truncation_threshold=1e-5,
    matrix_product_state_max_bond_dimension=100,
)
mps_sampler = Aer_Sampler.from_backend(mps_backend)

shots = 4096

data_sim = []
for j in range(points):
    circ_list = [
        circ.assign_parameters([params[j]]) for circ in circuits_unitary
    ]

    mps_job = mps_sampler.run(circ_list, shots=shots)
    result = mps_job.result()

    point_data = [
        result[d].data["meas"].expectation_values(observables_unitary)
        for d in depths
    ]

    data_sim.append(point_data)  # data at one theta value

data_sim = np.array(data_sim)

Con los circuitos y observables preparados, ahora los ejecutamos en hardware usando la primitiva Sampler.

Aquí enviamos tres trabajos para unitary_pubs, dynamic_pubs y dynamic_pubs_dd. Cada uno es una lista de circuitos parametrizados correspondientes a nueve pasos de Trotter diferentes con tres parámetros $\theta$ diferentes.

shots = 10000

with Batch(backend=backend) as batch:
    sampler = Sampler(mode=batch)

    sampler.options.experimental = {
        "execution": {
            "scheduler_timing": True
        },  # set to True to retrieve circuit timing info
    }

    job_unitary = sampler.run(unitary_pubs, shots=shots)
    print(f"unitary: {job_unitary.job_id()}")

    job_unitary_dd = sampler.run(unitary_pubs_dd, shots=shots)
    print(f"unitary_dd: {job_unitary_dd.job_id()}")

    job_dynamic = sampler.run(dynamic_pubs, shots=shots)
    print(f"dynamic: {job_dynamic.job_id()}")

    job_dynamic_dd = sampler.run(dynamic_pubs_dd, shots=shots)
    print(f"dynamic_dd: {job_dynamic_dd.job_id()}")

    job_dynamic_meas2 = sampler.run(dynamic_pubs_meas2, shots=shots)
    print(f"dynamic_meas2: {job_dynamic_meas2.job_id()}")

    job_dynamic_dd_meas2 = sampler.run(dynamic_pubs_dd_meas2, shots=shots)
    print(f"dynamic_dd_meas2: {job_dynamic_dd_meas2.job_id()}")

unitary: d5dtt0ldq8ts73fvbhj0
unitary: d5dtt11smlfc739onuag
dynamic: d5dtt1hsmlfc739onuc0
dynamic_dd: d5dtt25jngic73avdne0
dynamic_meas2: d5dtt2ldq8ts73fvbhm0
dynamic_dd_meas2: d5dtt2tjngic73avdnf0

Paso 4: Post-procesar y devolver resultados en el formato clásico deseado

Después de que los trabajos se hayan completado, podemos recuperar la duración del circuito de los metadatos de los resultados del trabajo y visualizar la información de programación del circuito. Para leer más sobre la visualización de la información de programación de un circuito, consulta esta página.

# Circuit durations is reported in the unit of `dt` which can be retrieved from `Backend` object
unitary_durations = [
    job_unitary.result()[i].metadata["compilation"]["scheduler_timing"][
        "circuit_duration"
    ]
    for i in depths
]

dynamic_durations = [
    job_dynamic.result()[i].metadata["compilation"]["scheduler_timing"][
        "circuit_duration"
    ]
    for i in depths
]

dynamic_durations_meas2 = [
    job_dynamic_meas2.result()[i].metadata["compilation"]["scheduler_timing"][
        "circuit_duration"
    ]
    for i in depths
]

result_dd = job_dynamic_dd.result()[1]
circuit_schedule_dd = result_dd.metadata["compilation"]["scheduler_timing"][
    "timing"
]

# to visualize the circuit schedule, one can show the figure below
fig_dd = draw_circuit_schedule_timing(
    circuit_schedule=circuit_schedule_dd,
    included_channels=None,
    filter_readout_channels=False,
    filter_barriers=False,
    width=1000,
)

# Save to a file since the figure is large
fig_dd.write_html("scheduler_timing_dd.html")

Graficamos las duraciones de los circuitos unitarios y los circuitos dinámicos. En el gráfico a continuación, podemos ver que, a pesar del tiempo necesario para las mediciones a mitad de circuito y las operaciones clásicas, la implementación con circuitos dinámicos con measure_2 resulta en duraciones de circuito comparables a la implementación unitaria.

# visualize circuit durations

def convert_dt_to_microseconds(circ_duration: List, backend_dt: float):
    dt = backend_dt * 1e6  # dt in microseconds
    return list(map(lambda x: x * dt, circ_duration))

dt = backend.target.dt
plt.plot(
    depths,
    convert_dt_to_microseconds(unitary_durations, dt),
    color="#be95ff",
    linestyle=":",
    label="unitary",
)
plt.plot(
    depths,
    convert_dt_to_microseconds(dynamic_durations, dt),
    color="#ff7eb6",
    linestyle="-.",
    label="dynamic",
)
plt.plot(
    depths,
    convert_dt_to_microseconds(dynamic_durations_meas2, dt),
    color="#ff7eb6",
    linestyle="-.",
    marker="s",
    mfc="none",
    label="dynamic w/ meas2",
)

plt.xlabel("Trotter steps")
plt.ylabel(r"Circuit durations in $\mu$s")
plt.legend()

<matplotlib.legend.Legend at 0x17f73c6e0>

Después de que los trabajos se hayan completado, recuperamos los datos a continuación y calculamos la magnetización promedio estimada por los observables observables_unitary u observables_dynamic que construimos anteriormente.

runs = {
    "unitary": (
        job_unitary,
        [observables_unitary] * len(circuits_unitary),
    ),
    "unitary_dd": (
        job_unitary_dd,
        [observables_unitary] * len(circuits_unitary),
    ),
    # Omitting Dyn w/o DD and Dynamic w/ DD plots for better readability
    # "dynamic": (job_dynamic, observables_dynamic),
    # "dynamic_dd": (job_dynamic_dd, observables_dynamic),
    "dynamic_meas2": (job_dynamic_meas2, observables_dynamic),
    "dynamic_dd_meas2": (
        job_dynamic_dd_meas2,
        observables_dynamic,
    ),
}

data_dict = {}
for key, (job, obs) in runs.items():
    data = []
    for i in range(points):
        data.append(
            [
                job.result()[ind].data["meas"].expectation_values(obs[ind])[i]
                for ind in depths
            ]
        )
    data_dict[key] = data

A continuación graficamos la magnetización de los espines como función de los pasos de Trotter para diferentes valores de $\theta$, correspondientes a diferentes intensidades del campo magnético local. Graficamos tanto los resultados de la simulación MPS precalculados para los circuitos unitarios ideales, junto con los resultados experimentales de lo siguiente:

1. ejecutar los circuitos unitarios con DD
2. ejecutar los circuitos dinámicos con DD y MidCircuitMeasure

plt.figure(figsize=(10, 6))

colors = ["#0f62fe", "#be95ff", "#ff7eb6"]
for i in range(points):
    plt.plot(
        depths,
        data_sim[i],
        color=colors[i],
        linestyle="solid",
        label=f"θ={pi_check(i*max_angle/(points-1))} (MPS)",
    )
    # plt.plot(
    #     depths,
    #     data_dict["unitary"][i],
    #     color=colors[i],
    #     linestyle=":",
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Unitary)",
    # )

    plt.plot(
        depths,
        data_dict["unitary_dd"][i],
        color=colors[i],
        marker="o",
        mfc="none",
        linestyle=":",
        label=f"θ={pi_check(i*max_angle/(points-1))} (Unitary w/DD)",
    )

    # Omitting Dyn w/o DD and Dynamic w/ DD plots for better readability
    # plt.plot(
    #     depths,
    #     data_dict["dynamic"][i],
    #     color=colors[i],
    #     linestyle="-.",
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Dyn w/o DD)",
    # )
    # plt.plot(
    #     depths,
    #     data_dict["dynamic_dd"][i],
    #     marker="D",
    #     mfc="none",
    #     color=colors[i],
    #     linestyle="-.",
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Dynamic w/ DD)",
    # )

    # plt.plot(
    #     depths,
    #     data_dict["dynamic_meas2"][i],
    #     color=colors[i],
    #     marker="s",
    #     mfc="none",
    #     linestyle=':',
    #     label=f"θ={pi_check(i*max_angle/(points-1))} (Dynamic w/ MidCircuitMeas)",
    # )

    plt.plot(
        depths,
        data_dict["dynamic_dd_meas2"][i],
        color=colors[i],
        marker="*",
        markersize=8,
        linestyle=":",
        label=f"θ={pi_check(i*max_angle/(points-1))} (Dynamic w/ DD & MidCircuitMeas)",
    )

plt.xlabel("Trotter steps", fontsize=16)
plt.ylabel("Average magnetization", fontsize=16)
plt.xticks(rotation=45)
handles, labels = plt.gca().get_legend_handles_labels()
plt.legend(
    handles,
    labels,
    loc="upper right",
    bbox_to_anchor=(1.46, 1.0),
    shadow=True,
    ncol=1,
)
plt.title(
    f"{hex_rows}x{hex_cols} hex ring, {num_qubits} data qubits, {len(ancilla)} ancilla qubits \n{backend.name}: Sampler"
)
plt.show()

Cuando comparamos los resultados experimentales con la simulación, vemos que la implementación con circuitos dinámicos (línea punteada con estrellas) en general tiene un mejor rendimiento que la implementación unitaria estándar (línea punteada con círculos). En resumen, presentamos los circuitos dinámicos como una solución para simular modelos de espín de Ising en una red de panal de abeja, una topología que no es nativa del hardware. La solución con circuitos dinámicos permite interacciones ZZ entre qubits que no son vecinos más cercanos, con una profundidad de puertas de dos qubits menor que usando puertas SWAP, a costa de introducir qubits ancilla adicionales y operaciones de retroalimentación clásica.

Referencias

[1] Quantum computing with Qiskit, by Javadi-Abhari, A., Treinish, M., Krsulich, K., Wood, C.J., Lishman, J., Gacon, J., Martiel, S., Nation, P.D., Bishop, L.S., Cross, A.W. and Johnson, B.R., 2024. arXiv preprint arXiv:2405.08810 (2024)