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

Estimación de uso: 7.5 minutos en un procesador Heron r3. (NOTA: Esto es solo una estimación. El tiempo de ejecución puede variar.

Los circuitos dinámicos son circuitos con alimentación directa clásica; en otras palabras, son mediciones a mitad del 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 efecto Kicked en una red hexagonal de espines y utilizamos circuitos dinámicos para realizar interacciones más allá de la conectividad física del hardware.

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

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 investigar la dinámica de espín, estudiamos la magnetización media de los espines en cada sitio en función de los pasos de Trotter. Por lo tanto, construimos la siguiente observable:

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

Para realizar la interacción ZZ entre los sitios de la red, presentamos una solución que utiliza la función de circuito dinámico, lo que da lugar a una profundidad de dos qubits significativamente menor en comparación con el método de enrutamiento estándar con puertas SWAP. Por otro lado, las operaciones clásicas de alimentación directa en circuitos dinámicos suelen tener tiempos de ejecución más largos que las puertas cuánticas; por lo tanto, los circuitos dinámicos tienen limitaciones y compensaciones. También presentamos una forma de añadir una secuencia de desacoplamiento dinámico en los qubits inactivos durante la operación clásica de alimentación directa utilizando la duración del estiramiento.

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 gráficos Rustworkx (pip install rustworkx)
Qiskit Aer (pip install qiskit-aer)

Configurar

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.providers.exceptions import QiskitBackendNotFoundError
from qiskit_ibm_runtime.visualization import (
    draw_circuit_schedule_timing,
)

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

Comenzamos definiendo la red que se va a simular. Elegimos trabajar con la red de panal (también llamada hexagonal), que es un grafo plano con nodos de grado 3. Aquí especificamos el tamaño de la red y los parámetros relevantes del circuito que nos interesan 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)

Output:

Utilizaremos el pequeño ejemplo para la ilustración y la simulación. A continuación, también construimos un ejemplo grande para mostrar que el flujo de trabajo se puede ampliar a tamaños grandes.

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

Output:

num_qubits = 46

Construir circuitos unitarios

Una vez especificados el tamaño del problema y los parámetros, ya 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 depth argumento. El circuito que construimos tiene capas alternas de puertas Rx ( $\theta$ ) y Rzz puertas. Las Rzz puertas realizan las interacciones ZZ entre espines acoplados, que se colocarán entre cada sitio de la red especificado por el layer_edges argumento.

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

Visualice el pequeño circuito de prueba:

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)

Output:

Del mismo modo, construya 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,
)

Implementación de circuitos dinámicos

Esta sección muestra la implementación del circuito dinámico principal para simular la misma evolución temporal trotterizada. Tenga en cuenta que la red de panal que queremos simular no coincide con la pesada red de los qubits de hardware. Una forma sencilla de mapear el circuito al hardware es introducir una serie de operaciones SWAP para colocar los qubits que interactúan uno al lado del otro, con el fin de realizar la interacción ZZ. Aquí destacamos un enfoque alternativo que utiliza circuitos dinámicos como solución, lo que ilustra que podemos utilizar la combinación de computación cuántica y computación clásica en tiempo real dentro de un circuito en Qiskit para realizar interacciones más allá del vecino más cercano.

En la implementación del circuito dinámico, la interacción ZZ se implementa de manera eficaz utilizando qubits auxiliares, mediciones a mitad del circuito y alimentación directa. Para entender esto, hay que tener en cuenta que las rotaciones ZZ aplican un factor de fase $e^{i\theta}$ al estado en función de su paridad. Para dos qubits, los estados de 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 sin cambios los estados de paridad par. A continuación se describe cómo podemos implementar eficazmente interacciones ZZ en dos qubits utilizando circuitos dinámicos.

Calcular la paridad en un qubit auxiliar: en lugar de aplicar directamente ZZ a dos qubits, introducimos un tercer qubit, el qubit auxiliar, para almacenar la información de paridad de los dos qubits de datos. Entrelazamos el ancilla con cada qubit de datos utilizando puertas CX desde el qubit de datos al qubit ancilla.
Aplica una rotación Z de un solo qubit al qubit auxiliar: esto se debe a que el qubit auxiliar tiene la información de paridad de los dos qubits de datos, lo que implementa eficazmente la rotación ZZ en los qubits de datos.
Mide el qubit auxiliar en la base X: este es el paso clave que colapsa el estado del qubit auxiliar, y el resultado de la medición nos dice lo que ha sucedido:
- Medida 0: cuando se observa un resultado 0, significa que hemos aplicado correctamente una rotación e $ZZ(\theta)$ a a nuestros qubits de datos.
- Medida 1: cuando se observa un resultado 1, hemos aplicado en su lugar $ZZ(\theta + \pi)$ .
Aplicar puerta de corrección al medir 1: Si medimos 1, aplicamos puertas Z a los qubits de datos para «corregir» la fase e $\pi$ e adicional.

El circuito resultante es el siguiente:

Cuando adoptamos este enfoque para simular una red de panal, el circuito resultante se integra perfectamente en el hardware con una red hexagonal densa: todos los qubits de datos residen en los sitios degree-3 de la red, que forma una red hexagonal. Cada par de qubits de datos comparte un qubit auxiliar que reside en un sitio e degree-2. A continuación, construimos la red de qubits para la implementación del circuito dinámico, introduciendo qubits auxiliares (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 hexagonal pesada para los qubits de datos y los qubits auxiliares 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()

Output:

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 RZZ puertas se sustituyen por la implementación del circuito dinámico utilizando 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 `stretch` duración

Una advertencia sobre el uso de la implementación de circuitos dinámicos para realizar la interacción ZZ es que la medición en medio del circuito y las operaciones clásicas de alimentación directa suelen tardar más tiempo en ejecutarse que las puertas cuánticas. Para suprimir la decoherencia de los qubits durante el tiempo de inactividad necesario para que se produzcan las operaciones clásicas, añadimos una secuencia de desacoplamiento dinámico (DD) después de la operación de medición en los qubits auxiliares y antes de la operación Z condicional en el qubit de datos, antes de la if_test instrucción.

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

Utilizando las funciones definidas anteriormente, construimos los circuitos de evolución temporal trotterizados, con y sin DD, y los observables correspondientes.

Comenzamos visualizando el circuito dinámico de un pequeño ejemplo:

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)

Output:

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)

Output:

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)

Output:

Del mismo modo, construya 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 del hardware

Ahora estamos listos para transpilear el circuito al hardware. Transpilamos tanto la implementación del estándar unitario como la implementación del circuito dinámico al hardware.

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

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

En primer lugar, transpilamos los circuitos dinámicos, con y sin añadir la secuencia DD. Para garantizar que utilizamos el mismo conjunto de qubits físicos en todos los circuitos y obtener resultados más consistentes, primero transpilamos el circuito una vez y, a continuación, utilizamos su diseño para todos los circuitos posteriores, especificados por initial_layout en el gestor de pasadas. A continuación, construimos los bloques primitivos unificados (PUB) como entrada primitiva del muestreador.

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 los qubits del circuito transpuesto a continuación. Los círculos negros muestran los qubits de datos y los qubits auxiliares utilizados en la implementación del circuito dinámico.

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",
)

Output:

Note

Si obtienes errores sobre neato no encontrado desde plot_circuit_layout(), asegúrate de tener el graphviz paquete instalado y disponible en tu PATH. Si se instala en una ubicación no predeterminada (por ejemplo, utilizando homebrew en MacOS ), es posible que tenga que actualizar su variable PATH de entorno. Esto se puede hacer dentro de este cuaderno utilizando 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)

Output:

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

Output:

Transpilar utilizando `MidCircuitMeasure`

MidCircuitMeasure es una ampliación de las operaciones de medición disponibles, calibrada específicamente para realizar mediciones en el punto medio del circuito. La MidCircuitMeasure instrucción se traduce a la measure_2 instrucción compatible con los backends. Ten en cuenta que no measure_2 es compatible con todos los backends. Puedes utilizar service.backends(filters=lambda b: "measure_2" in b.supported_instructions) para encontrar backends que lo admitan. A continuación, mostramos cómo transpilar el circuito para que las mediciones en tiempo real definidas en él se ejecuten mediante la MidCircuitMeasure operación, siempre que el backend la admita.

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

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

Output:

Mid-circuit measurement `measure_2` duration: 1.3800000000000003 μs
Terminal measurement `measure` duration: 2.1800000000000006 μ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, utilizamos el mismo conjunto de qubits físicos utilizados en los circuitos dinámicos para los qubits de datos como diseño para transpilear 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 los qubits de los circuitos unitarios transpilados. Los círculos negros indican los qubits físicos utilizados para transpilar los circuitos unitarios y sus índices corresponden a los índices de los qubits virtuales. Al comparar esto con el diseño trazado para los circuitos dinámicos, podemos confirmar que los circuitos unitarios utilizan 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",
)

Output:

Ahora añadimos la secuencia DD a los circuitos transpilados y construimos los PUB 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]

Compara la profundidad de puerta de dos qubits de 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()

Output:

<matplotlib.legend.Legend at 0x12628b0e0>

La principal ventaja del circuito basado en mediciones es que, al implementar múltiples interacciones ZZ, las capas CX se pueden paralelizar y las mediciones se pueden realizar simultáneamente. Esto se debe a que todas las interacciones ZZ son conmutativas, por lo que el cálculo se puede realizar con una profundidad de medición 1. Tras transpilación de los circuitos, observamos que el enfoque de circuito dinámico produce una profundidad de dos qubits significativamente menor que el enfoque unitario estándar, con la salvedad de que la medición adicional a mitad del circuito y la alimentación directa clásica requieren tiempo e introducen sus propias fuentes de error.

Paso 3: Ejecutar utilizando Qiskit primitives

Modo de prueba local

Antes de enviar los trabajos al hardware, podemos ejecutar una pequeña simulación de prueba del circuito dinámico utilizando 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])

Output:

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

array([ 0.16666667,  0.01529948, -0.14290365])

Simulación MPS

Para circuitos grandes, podemos utilizar el simulador matrix_product_state (MPS), que proporciona un resultado aproximado al valor esperado según la dimensión de enlace elegida. Posteriormente, utilizamos los resultados de la simulación MPS como referencia 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 utilizando la primitiva Sampler.

Aquí presentamos tres trabajos para unitary_pubs, dynamic_pubs, y dynamic_pubs_dd. Cada uno es una lista de circuitos parametrizados correspondientes a nueve pasos Trotter diferentes con tres parámetros e $\theta$ es 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()}")

Output:

unitary: d96s4b52su3c739hakrg
unitary_dd: d96s4bt2su3c739haksg
dynamic: d96s4c0tcv6s73dk55mg
dynamic_dd: d96s4ckqp3as739qvid0
dynamic_meas2: d96s4csqp3as739qvie0
dynamic_dd_meas2: d96s4daf47jc73a5v8s0

Paso 4: Procesar posteriormente y devolver los resultados en el formato clásico deseado

Una vez finalizados los trabajos, podemos obtener la duración del circuito a partir de los metadatos de los resultados de los trabajos y visualizar la información sobre el calendario del circuito. Para obtener más información sobre cómo visualizar 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")

Trazamos las duraciones de los circuitos unitarios y los circuitos dinámicos. En el gráfico siguiente podemos ver que, a pesar del tiempo necesario para las mediciones a mitad del circuito y las operaciones clásicas, la implementación dinámica del circuito con measure_2 da como resultado duraciones del circuito comparables a las de 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()

Output:

<matplotlib.legend.Legend at 0x12bfde270>

Una vez completados los trabajos, recuperamos los datos que se muestran a continuación y calculamos la magnetización media estimada por las observables observables_unitary o observables_dynamic que hemos construido 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, representamos gráficamente la magnetización de espín en función de los pasos de Trotter a diferentes valores de $\theta$ , correspondientes a diferentes intensidades del campo magnético local. Trazamos tanto los resultados de la simulación MPS precalculados para los circuitos ideales unitarios, como los resultados experimentales de lo siguiente:

ejecutar los circuitos unitarios con DD
ejecutando 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))} "
        f"(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, "
    f"{len(ancilla)} ancilla qubits \n{backend.name}: Sampler"
)
plt.show()

Output:

Cuando comparamos los resultados experimentales con la simulación, vemos que la implementación del circuito dinámico (línea punteada con estrellas) tiene, en general, 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 hexagonal, una topología que no es nativa del hardware. La solución de circuito dinámico permite interacciones ZZ entre qubits que no son vecinos más cercanos, con una profundidad de puerta de dos qubits más corta que la que se obtiene con puertas SWAP, a costa de introducir qubits auxiliares adicionales y operaciones clásicas de alimentación directa.

Referencias

[1] Computación cuántica con Qiskit, por 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. y Johnson, B.R 2024. « arXiv » ( Preimpresión) arXiv:2405.08810 (2024)

¿Le ha resultado útil esta página?

Informe de un error, de una errata o solicite contenido en GitHub.

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

Requisitos

Configurar

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

Construir circuitos unitarios

Implementación de circuitos dinámicos

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

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

Transpilación para circuitos dinámicos

Transpilar utilizando MidCircuitMeasure

Transpilación para circuitos unitarios

Compara la profundidad de puerta de dos qubits de circuitos unitarios y dinámicos

Paso 3: Ejecutar utilizando Qiskit primitives

Modo de prueba local

Simulación MPS

Paso 4: Procesar posteriormente y devolver los resultados en el formato clásico deseado

Referencias

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

Transpilar utilizando `MidCircuitMeasure`