Simulation de l'hamiltonien d'Ising avec circuits dynamiques

Estimation d'utilisation : 7.5 minutes sur un processeur Heron r3. (REMARQUE : Il s'agit uniquement d'une estimation. Votre durée d'exécution peut varier.)

Les circuits dynamiques sont des circuits à action directe classique. En d'autres termes, il s'agit de mesures effectuées à mi-parcours, suivies d'opérations logiques classiques qui déterminent les opérations quantiques en fonction du résultat classique. Dans ce tutoriel, nous simulons le modèle d'Ising kické sur un réseau hexagonal de spins et utilisons des circuits dynamiques pour réaliser des interactions au-delà de la connectivité physique du matériel.

Le modèle d'Ising a fait l'objet de nombreuses études dans différents domaines de la physique. Il modélise les spins qui subissent des interactions d'Ising entre les sites du réseau, ainsi que les impulsions provenant du champ magnétique local sur chaque site. L'évolution temporelle trottisée des spins considérés dans ce tutoriel, tirée de [1], est donnée par l'unité suivante :

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)

Pour étudier la dynamique des spins, nous étudions la magnétisation moyenne des spins à chaque site en fonction des étapes de Trotter. Nous construisons donc l'observable suivant :

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

Pour réaliser l'interaction ZZ entre les sites du réseau, nous présentons une solution utilisant la fonctionnalité de circuit dynamique, qui permet d'obtenir une profondeur de deux qubits nettement plus courte par rapport à la méthode de routage standard avec des portes SWAP. D'autre part, les opérations classiques de feedforward dans les circuits dynamiques ont généralement des temps d'exécution plus longs que les portes quantiques; les circuits dynamiques présentent donc des limites et des compromis. Nous présentons également une méthode permettant d'ajouter une séquence de découplage dynamique sur les qubits inactifs pendant l'opération classique de feedforward en utilisant la durée d'étirement.

Exigences

Avant de commencer ce tutoriel, assurez-vous que les éléments suivants sont installés :

Qiskit SDK v2.0 ou ultérieurement avec prise en charge de la visualisation
Qiskit Runtime v0.37 ou ultérieurement avec prise en charge de la visualisation (pip install 'qiskit-ibm-runtime[visualization]')
Bibliothèque graphique Rustworkx (pip install rustworkx)
Qiskit Aer (pip install qiskit-aer)

Configurer

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

Étape 1 : Mapper les entrées classiques à un circuit quantique

Nous commençons par définir le réseau à simuler. Nous avons choisi de travailler avec le réseau en nid d'abeille (également appelé hexagonal), qui est un graphe plan avec des nœuds de degré 3. Ici, nous spécifions la taille du réseau, les paramètres de circuit pertinents qui nous intéressent dans la dynamique trotterisée. Nous simulons l'évolution temporelle selon le modèle de Trotter sous le modèle d'Ising avec trois valeurs d' $\theta$ s différentes du champ magnétique 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

Commençons par un petit exemple de test :

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:

Nous utiliserons ce petit exemple à des fins d'illustration et de simulation. Ci-dessous, nous construisons également un exemple à grande échelle afin de montrer que le flux de travail peut être étendu à des tailles importantes.

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

Construire des circuits unitaires

Une fois la taille du problème et les paramètres spécifiés, nous sommes désormais prêts à construire le circuit paramétré qui simule l'évolution temporelle de l' $U(\theta)$ , selon la méthode de Trotter, avec différentes étapes de Trotter, spécifiées par depth l'argument. Le circuit que nous construisons comporte des couches alternées de portes Rx ( $\theta$ ) et de Rzz portes. Les Rzz portes réalisent les interactions ZZ entre les spins couplés, qui seront placés entre chaque site du réseau spécifié par layer_edges l'argument.

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

Visualisez le petit circuit de test :

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:

De même, construisez les circuits unitaires du grand exemple à différentes étapes de Trotter et l'observable pour estimer la valeur attendue.

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

Construire une implémentation de circuit dynamique

Cette section présente la mise en œuvre du circuit dynamique principal permettant de simuler la même évolution temporelle trotterisée. Notez que le réseau en nid d'abeille que nous voulons simuler ne correspond pas au réseau lourd des qubits matériels. Une manière simple de mapper le circuit sur le matériel consiste à introduire une série d'opérations SWAP afin de rapprocher les qubits en interaction les uns des autres, afin de réaliser l'interaction ZZ. Nous mettons ici en avant une approche alternative utilisant des circuits dynamiques comme solution, qui montre que nous pouvons utiliser la combinaison du calcul quantique et du calcul classique en temps réel dans un circuit dans Qiskit pour réaliser des interactions au-delà du voisinage immédiat.

Dans la mise en œuvre du circuit dynamique, l'interaction ZZ est efficacement mise en œuvre à l'aide de qubits auxiliaires, de mesures à mi-circuit et d'une prédiction. Pour comprendre cela, notez que les rotations ZZ appliquent un facteur de phase $e^{i\theta}$ à l'état en fonction de sa parité. Pour deux qubits, les états de base computationnelle sont $|00\rangle$ , $|01\rangle$ , $|10\rangle$ et $|11\rangle$ . La porte de rotation ZZ applique un facteur de phase aux états $|01\rangle$ et $|10\rangle$ dont la parité (le nombre de uns dans l'état) est impaire et laisse les états de parité paire inchangés. Ce qui suit décrit comment nous pouvons mettre en œuvre efficacement des interactions ZZ sur deux qubits à l'aide de circuits dynamiques.

Calculer la parité dans un qubit auxiliaire : au lieu d'appliquer directement ZZ à deux qubits, nous introduisons un troisième qubit, le qubit auxiliaire, pour stocker les informations de parité des deux qubits de données. Nous entrelacons l'ancilla avec chaque qubit de données à l'aide de portes CX allant du qubit de données au qubit ancilla.
Appliquez une rotation Z à un seul qubit à la qubit auxiliaire : en effet, la qubit auxiliaire contient les informations de parité des deux qubits de données, ce qui permet de mettre en œuvre efficacement la rotation ZZ sur les qubits de données.
Mesurez le qubit auxiliaire dans la base X : c'est l'étape clé qui fait s'effondrer l'état du qubit auxiliaire, et le résultat de la mesure nous indique ce qui s'est passé :
- Mesure 0 : lorsqu'un résultat 0 est observé, cela signifie que nous avons correctement appliqué une rotation d' $ZZ(\theta)$ s à nos qubits de données.
- Mesure 1 : lorsqu'un résultat 1 est observé, nous avons appliqué à la place l' $ZZ(\theta + \pi)$ .
Appliquer une porte de correction lors de la mesure 1 : si nous avons mesuré 1, nous appliquons des portes Z aux qubits de données pour « corriger » la phase d' $\pi$ e supplémentaire.

Le circuit obtenu est le suivant :

Lorsque nous adoptons cette approche pour simuler un réseau en nid d'abeille, le circuit résultant s'intègre parfaitement dans le matériel avec un réseau hexagonal dense : tous les qubits de données résident sur les sites d' degree-3 s du réseau, qui forme un réseau hexagonal. Chaque paire de qubits de données partage un qubit auxiliaire résidant sur un site d' degree-2. Ci-dessous, nous construisons le réseau de qubits pour la mise en œuvre du circuit dynamique, en introduisant des qubits auxiliaires (représentés par les cercles violets plus foncés).

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)

Visualisez le réseau hexagonal dense pour les qubits de données et les qubits auxiliaires à petite échelle :

(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

Ci-dessous, nous construisons le circuit dynamique pour l'évolution temporelle trotterisée. Les RZZ portes sont remplacées par la mise en œuvre du circuit dynamique à l'aide des étapes décrites ci-dessus.

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

Découplage dynamique (DD) et prise en charge de `stretch` la durée

Une mise en garde concernant l'utilisation de la mise en œuvre de circuits dynamiques pour réaliser l'interaction ZZ est que la mesure à mi-circuit et les opérations classiques de feedforward prennent généralement plus de temps à exécuter que les portes quantiques. Afin de supprimer la décohérence des qubits pendant le temps d'inactivité nécessaire à la réalisation des opérations classiques, nous avons ajouté une séquence de découplage dynamique (DD) après l'opération de mesure sur les qubits ancilla, et avant l'opération Z conditionnelle sur le qubit de données, avant if_test l'instruction.

La séquence DD est générée par la fonction add_stretch_dd(), qui utilise les durées stretch pour déterminer les intervalles de temps entre les impulsions DD. Une durée stretch permet de définir une durée flexible pour l'opération delay , de sorte que la durée du délai puisse s'allonger jusqu'à occuper tout le temps d'inactivité du qubit. Les variables de durée spécifiées par stretch sont converties, lors de la compilation, en durées souhaitées qui satisfont à une certaine contrainte. Cela s'avère très utile lorsque la synchronisation des séquences DD est essentielle pour obtenir de bonnes performances en matière de suppression des erreurs. Pour plus d'informations sur ce type stretch , consultez la documentation de OpenQASM. À l'heure actuelle, la prise en charge de ce type stretch est encore au stade expérimental. Pour plus de détails sur ses contraintes d'utilisation, veuillez vous reporter à la section « Limitations » de la documentation stretch .

À l'aide des fonctions définies ci-dessus, nous construisons les circuits d'évolution temporelle trotterisés, avec et sans DD, ainsi que les observables correspondants.

Nous commençons par visualiser le circuit dynamique d'un petit exemple :

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:

De même, construisez les circuits dynamiques pour l'exemple de grande taille :

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)

Étape 2 : Optimiser le problème pour l'exécution matérielle

Nous sommes maintenant prêts à transposer le circuit vers le matériel. Nous transposerons à la fois l'implémentation standard unitaire et l'implémentation dynamique du circuit vers le matériel.

Pour transcompiler vers le matériel, nous instancions d'abord le backend. Si disponible, nous choisirons un backend prenant en charge l'instruction MidCircuitMeasure (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,
    )

Transpilation pour circuits dynamiques

Tout d'abord, nous transpilons les circuits dynamiques, avec et sans ajout de la séquence DD. Afin de garantir l'utilisation du même ensemble de qubits physiques dans tous les circuits pour obtenir des résultats plus cohérents, nous transpilons d'abord le circuit une fois, puis utilisons sa disposition pour tous les circuits suivants, spécifiés par initial_layout dans le gestionnaire de passes. Nous construisons ensuite les blocs unifiés primitifs (PUB) comme entrée primitive de l'échantillonneur.

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]

Nous pouvons visualiser la disposition des qubits du circuit transpilé ci-dessous. Les cercles noirs représentent les qubits de données et les qubits auxiliaires utilisés dans la mise en œuvre du circuit dynamique.

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 vous obtenez des erreurs indiquant que neato n'est pas trouvé dans plot_circuit_layout(), assurez-vous que le graphviz paquet est installé et disponible dans votre PATH. Si l'installation s'effectue dans un emplacement autre que celui par défaut (par exemple, en utilisant homebrew sur MacOS ), vous devrez peut-être mettre à jour votre variable PATH d'environnement. Cela peut être fait dans ce cahier à l'aide de la commande suivante :

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:

Transpiler à l'aide de `MidCircuitMeasure`

MidCircuitMeasure vient compléter les opérations de mesure disponibles; il est spécialement calibré pour effectuer des mesures en cours de circuit. L'instruction MidCircuitMeasure correspond à measure_2 l'instruction prise en charge par les backends. Notez que cette measure_2 fonctionnalité n'est pas prise en charge par tous les backends. Vous pouvez utiliser service.backends(filters=lambda b: "measure_2" in b.supported_instructions) pour trouver les backends qui le prennent en charge. Nous montrons ici comment transcompiler le circuit de manière à ce que les mesures intermédiaires définies dans le circuit soient exécutées à l'aide de MidCircuitMeasure l'opération, si le backend la prend en charge.

Ci-dessous, nous imprimons la durée de measure_2 l'instruction et de l'instruction measure standard.

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
]

Transpilation pour circuits unitaires

Afin d'établir une comparaison équitable entre les circuits dynamiques et leur équivalent unitaire, nous utilisons le même ensemble de qubits physiques que celui utilisé dans les circuits dynamiques pour les qubits de données comme disposition pour la transposition des circuits unitaires.

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]

Nous visualisons la disposition des qubits des circuits unitaires transposés. Les cercles noirs indiquent les qubits physiques utilisés pour transposer les circuits unitaires et leurs indices correspondent aux indices des qubits virtuels. En comparant cela avec la disposition tracée pour les circuits dynamiques, nous pouvons confirmer que les circuits unitaires utilisent le même ensemble de qubits physiques que les qubits de données dans les circuits dynamiques.

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

Output:

Nous ajoutons maintenant la séquence DD aux circuits transpilés et construisons les PUB correspondants pour la soumission des tâches.

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]

Comparer la profondeur de porte à deux qubits des circuits unitaires et dynamiques

# 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>

Le principal avantage du circuit basé sur la mesure est que, lors de la mise en œuvre de multiples interactions ZZ, les couches CX peuvent être parallélisées et les mesures peuvent être effectuées simultanément. En effet, toutes les interactions ZZ commutent, ce qui permet d'effectuer le calcul avec une profondeur de mesure 1. Après avoir transpilé les circuits, nous observons que l'approche par circuit dynamique produit une profondeur à deux qubits nettement plus courte que l'approche unitaire standard, avec toutefois la réserve que la mesure supplémentaire à mi-circuit et la rétroaction classique prennent elles-mêmes du temps et introduisent leurs propres sources d'erreurs.

Étape 3 : Exécutez à l'aide d' Qiskit primitives

Mode de test local

Avant de soumettre les tâches au matériel, nous pouvons effectuer une petite simulation test du circuit dynamique à l'aide du mode de test 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])

Simulation MPS

Pour les circuits de grande taille, nous pouvons utiliser le simulateur matrix_product_state (MPS), qui fournit un résultat approximatif de la valeur attendue en fonction de la dimension de liaison choisie. Nous utilisons ensuite les résultats de la simulation MPS comme référence pour comparer les résultats obtenus à partir du matériel.

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

Une fois les circuits et les observables préparés, nous les exécutons désormais sur le matériel à l'aide de la primitive Sampler.

Nous soumettons ici trois tâches pour unitary_pubs, dynamic_pubs, et dynamic_pubs_dd. Chacune est une liste de circuits paramétrés correspondant à neuf étapes Trotter différentes avec trois paramètres d' $\theta$ s différents.

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

Étape 4 : Post-traitement et restitution des résultats dans le format classique souhaité

Une fois les tâches terminées, nous pouvons extraire la durée du circuit à partir des métadonnées des résultats des tâches et visualiser les informations relatives au calendrier du circuit. Pour en savoir plus sur la visualisation des informations de planification d'un circuit, consultez cette page.

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

Nous représentons graphiquement les durées des circuits unitaires et des circuits dynamiques. Le graphique ci-dessous montre que, malgré le temps nécessaire pour les mesures à mi-circuit et les opérations classiques, la mise en œuvre dynamique du circuit avec measure_2 donne des durées de circuit comparables à celles de la mise en œuvre unitaire.

# 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>

Une fois les tâches terminées, nous récupérons les données ci-dessous et calculons la magnétisation moyenne estimée par les observables observables_unitary ou que observables_dynamic nous avons construits précédemment.

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

Ci-dessous, nous représentons graphiquement la magnétisation de spin en fonction des pas de Trotter à différentes valeurs d' $\theta$ , correspondant à différentes intensités du champ magnétique local. Nous représentons graphiquement les résultats précalculés de la simulation MPS pour les circuits idéaux unitaires, ainsi que les résultats expérimentaux suivants :

Faire fonctionner les circuits unitaires avec DD
faire fonctionner les circuits dynamiques avec DD et 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:

Lorsque nous comparons les résultats expérimentaux avec la simulation, nous constatons que la mise en œuvre du circuit dynamique (ligne pointillée avec des étoiles) offre globalement de meilleures performances que la mise en œuvre unitaire standard (ligne pointillée avec des cercles). En résumé, nous présentons les circuits dynamiques comme une solution pour simuler les modèles de spin d'Ising sur un réseau en nid d'abeille, une topologie qui n'est pas native du matériel. La solution de circuit dynamique permet des interactions ZZ entre des qubits qui ne sont pas des voisins immédiats, avec une profondeur de porte à deux qubits plus courte que celle obtenue avec des portes SWAP, au prix de l'introduction de qubits auxiliaires supplémentaires et d'opérations classiques de feedforward.

Références

[1] L'informatique quantique avec Qiskit, par Javadi-Abhari, A., Treinish, M., Krsulich, K., Wood, C.J Lishman, J., Gacon, J., Martiel, S., Nation, P.D Évêque, L.S Cross, A.W. et Johnson, B.R., 2024. arXiv prépublication arXiv:2405.08810 (2024)

Cette page a-t-elle été utile ?

Signaler un bogue, une coquille ou proposer du contenu sur GitHub.

Simulation de l'hamiltonien d'Ising avec circuits dynamiques

Exigences

Configurer

Étape 1 : Mapper les entrées classiques à un circuit quantique

Construire des circuits unitaires

Construire une implémentation de circuit dynamique

Découplage dynamique (DD) et prise en charge de stretch la durée

Étape 2 : Optimiser le problème pour l'exécution matérielle

Transpilation pour circuits dynamiques

Transpiler à l'aide de MidCircuitMeasure

Transpilation pour circuits unitaires

Comparer la profondeur de porte à deux qubits des circuits unitaires et dynamiques

Étape 3 : Exécutez à l'aide d' Qiskit primitives

Mode de test local

Simulation MPS

Étape 4 : Post-traitement et restitution des résultats dans le format classique souhaité

Références

Découplage dynamique (DD) et prise en charge de `stretch` la durée

Transpiler à l'aide de `MidCircuitMeasure`