Simulazione dell'Hamiltoniano di Ising con circuiti dinamici

Stima di utilizzo: 7.5 minuti su un processore Heron r3. (NOTA: Si tratta solo di una stima. Il tempo di esecuzione può variare.)

I circuiti dinamici sono circuiti con feedforward classico, ovvero misurazioni effettuate a metà circuito seguite da operazioni logiche classiche che determinano operazioni quantistiche condizionate dall'output classico. In questo tutorial, simuliamo il modello di Ising con kick su un reticolo esagonale di spin e utilizziamo circuiti dinamici per realizzare interazioni che vanno oltre la connettività fisica dell'hardware.

Il modello di Ising è stato studiato approfonditamente in diversi ambiti della fisica. Modella gli spin che subiscono interazioni di Ising tra i siti del reticolo, nonché gli impulsi provenienti dal campo magnetico locale su ciascun sito. L'evoluzione temporale trotterizzata degli spin considerati in questo tutorial, tratta da [1], è data dalla seguente unità:

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)

Per sondare la dinamica degli spin, studiamo la magnetizzazione media degli spin in ciascun sito in funzione dei passi di Trotter. Quindi, costruiamo la seguente osservabile:

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

Per realizzare l'interazione ZZ tra i siti del reticolo, presentiamo una soluzione che utilizza la funzione di circuito dinamico, che porta a una profondità di due qubit significativamente più breve rispetto al metodo di instradamento standard con porte SWAP. D'altra parte, le operazioni feedforward classiche nei circuiti dinamici hanno tipicamente tempi di esecuzione più lunghi rispetto ai gate quantistici; pertanto, i circuiti dinamici presentano limitazioni e compromessi. Presentiamo inoltre un modo per aggiungere una sequenza di disaccoppiamento dinamico sui qubit inattivi durante l'operazione di feedforward classica utilizzando la durata dell'allungamento.

Requisiti

Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:

Qiskit SDK v2.0 o versioni successive con supporto alla visualizzazione
Qiskit Runtime v0.37 o versioni successive con supporto alla visualizzazione (pip install 'qiskit-ibm-runtime[visualization]')
Libreria grafica Rustworkx (pip install rustworkx)
Qiskit Aer (pip install qiskit-aer)

Impostazione di

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

Fase 1: mappare gli input classici su un circuito quantistico

Iniziamo definendo il reticolo da simulare. Abbiamo scelto di lavorare con il reticolo a nido d'ape (chiamato anche esagonale), che è un grafo planare con nodi di grado 3. Qui specifichiamo la dimensione del reticolo, i parametri circuitali rilevanti di interesse nella dinamica trotterizzata. Simuliamo l'evoluzione temporale trotterizzata secondo il modello di Ising con tre diversi valori di campo magnetico locale ( $\theta$ ).

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

Cominciamo con un piccolo esempio di prova:

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:

Useremo il piccolo esempio per illustrazione e simulazione. Di seguito riportiamo anche un esempio di grandi dimensioni per mostrare che il flusso di lavoro può essere esteso a dimensioni maggiori.

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

Costruire circuiti unitari

Una volta specificati la dimensione del problema e i parametri, siamo ora pronti a costruire il circuito parametrizzato che simula l'evoluzione temporale trotterizzata di $U(\theta)$ con diversi passi di Trotter, specificati depth dall'argomento. Il circuito che costruiamo presenta strati alternati di porte Rx ( $\theta$ ) e Rzz porte. I Rzz gate realizzano le interazioni ZZ tra gli spin accoppiati, che saranno posizionati tra ciascun sito del reticolo specificato layer_edges dall'argomento.

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

Visualizza il piccolo circuito di prova:

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:

Allo stesso modo, costruisci i circuiti unitari dell'esempio grande a diversi passi di Trotter e l'osservabile per stimare il valore atteso.

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

Costruire un'implementazione dinamica del circuito

Questa sezione illustra l'implementazione del circuito dinamico principale per simulare la stessa evoluzione temporale trotterizzata. Si noti che il reticolo a nido d'ape che vogliamo simulare non corrisponde al reticolo pesante dei qubit hardware. Un modo semplice per mappare il circuito sull'hardware consiste nell'introdurre una serie di operazioni SWAP per avvicinare i qubit interagenti tra loro, in modo da realizzare l'interazione ZZ. Qui mettiamo in evidenza un approccio alternativo che utilizza circuiti dinamici come soluzione, il quale dimostra che è possibile utilizzare la combinazione di calcolo quantistico e calcolo classico in tempo reale all'interno di un circuito in Qiskit per realizzare interazioni oltre il vicino più prossimo.

Nell'implementazione del circuito dinamico, l'interazione ZZ viene efficacemente implementata utilizzando qubit ancilla, misurazione a metà circuito e feedforward. Per comprendere questo concetto, occorre tenere presente che le rotazioni ZZ applicano un fattore di fase $e^{i\theta}$ allo stato in base alla sua parità. Per due qubit, gli stati di base computazionali sono $|00\rangle$ , $|01\rangle$ , $|10\rangle$ e $|11\rangle$ . Il gate di rotazione ZZ applica un fattore di fase agli stati $|01\rangle$ e $|10\rangle$ la cui parità (il numero di uno nello stato) è dispari e lascia invariati gli stati di parità pari. Di seguito viene descritto come implementare efficacemente le interazioni ZZ su due qubit utilizzando circuiti dinamici.

Calcolare la parità in un qubit ancilla: invece di applicare direttamente ZZ a due qubit, introduciamo un terzo qubit, il qubit ancilla, per memorizzare le informazioni di parità dei due qubit di dati. Intrecciamo l'ancilla con ciascun qubit di dati utilizzando porte CX dal qubit di dati al qubit ancilla.
Applicare una rotazione Z a un singolo qubit all'ancilla qubit: questo perché l'ancilla possiede le informazioni di parità dei due qubit di dati, il che implementa efficacemente la rotazione ZZ sui qubit di dati.
Misura il qubit ancilla nella base X: questo è il passaggio chiave che fa crollare lo stato del qubit ancilla, e il risultato della misurazione ci dice cosa è successo:
- Misura 0: quando si osserva un risultato pari a 0, significa che abbiamo applicato correttamente una rotazione dell' $ZZ(\theta)$ e ai nostri qubit di dati.
- Misura 1: quando si osserva un risultato pari a 1, abbiamo applicato invece l' $ZZ(\theta + \pi)$ .
Applicare il gate di correzione durante la misurazione 1: se abbiamo misurato 1, applichiamo i gate Z ai qubit di dati per "correggere" la fase extra $\pi$ .

Il circuito risultante è il seguente:

Quando adottiamo questo approccio per simulare un reticolo a nido d'ape, il circuito risultante si integra perfettamente nell'hardware con un reticolo esagonale pesante: tutti i qubit di dati risiedono sui siti degree-3 del reticolo, che forma un reticolo esagonale. Ogni coppia di qubit di dati condivide un qubit ancillare che risiede in un sito dell' degree-2. Di seguito, costruiamo il reticolo di qubit per l'implementazione del circuito dinamico, introducendo i qubit ancilla (mostrati nei cerchi viola più scuri).

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)

Visualizza il reticolo esagonale pesante per i qubit di dati e i qubit ancilla su piccola scala:

(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

Di seguito, costruiamo il circuito dinamico per l'evoluzione temporale trotterizzata. I RZZ gate vengono sostituiti con l'implementazione del circuito dinamico utilizzando i passaggi descritti sopra.

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

Disaccoppiamento dinamico (DD) e supporto per `stretch` la durata

Un avvertimento sull'uso dell'implementazione del circuito dinamico per realizzare l'interazione ZZ è che la misurazione a metà circuito e le operazioni classiche di feedforward richiedono in genere più tempo per essere eseguite rispetto ai gate quantistici. Per sopprimere la decoerenza dei qubit durante il tempo di inattività necessario per l'esecuzione delle operazioni classiche, abbiamo aggiunto una sequenza di disaccoppiamento dinamico (DD) dopo l'operazione di misurazione sui qubit ancilla e prima dell'operazione Z condizionale sul qubit dati, prima dell'istruzione if_test .

La sequenza DD viene generata dalla funzione add_stretch_dd(), che utilizza le durate stretch per determinare gli intervalli di tempo tra i gate DD. Una durata stretch è un modo per specificare un intervallo di tempo estensibile per l'operazione delay , in modo tale che la durata del ritardo possa aumentare fino a coprire il tempo di inattività del qubit. Le variabili di durata specificate da stretch vengono risolte in fase di compilazione nelle durate desiderate che soddisfano un determinato vincolo. Ciò risulta molto utile quando la temporizzazione delle sequenze DD è fondamentale per ottenere buone prestazioni in termini di soppressione degli errori. Per ulteriori dettagli su questo tipo stretch , consultare la documentazione di OpenQASM. Al momento, il supporto per questo tipo stretch è in fase sperimentale. Per ulteriori dettagli sui vincoli di utilizzo, si prega di consultare la sezione dedicata alle limitazioni nella documentazione stretch .

Utilizzando le funzioni definite sopra, costruiamo i circuiti di evoluzione temporale trotterizzati, con e senza DD, e le corrispondenti osservabili.

Iniziamo visualizzando il circuito dinamico di un piccolo esempio:

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:

Allo stesso modo, costruisci i circuiti dinamici per l'esempio più 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)

Fase 2: Ottimizzazione del problema per l'esecuzione hardware

Ora siamo pronti per trasporre il circuito sull'hardware. Trasformeremo sia l'implementazione dello standard unitario che l'implementazione del circuito dinamico nell'hardware.

Per eseguire la transpilazione sull'hardware, prima istanziamo il backend. Se disponibile, sceglieremo un backend che supporti l'istruzione 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,
    )

Transpilazione per circuiti dinamici

Per prima cosa, transpiliamo i circuiti dinamici, con e senza l'aggiunta della sequenza DD. Per garantire l'utilizzo dello stesso set di qubit fisici in tutti i circuiti e ottenere risultati più coerenti, prima transpiliamo il circuito una volta, quindi utilizziamo il suo layout per tutti i circuiti successivi, specificati da initial_layout nel pass manager. Quindi costruiamo i blocchi primitivi unificati (PUB) come input primitivo del campionatore.

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]

Possiamo visualizzare la disposizione dei qubit del circuito transpilato qui sotto. I cerchi neri indicano i qubit di dati e i qubit ancilla utilizzati nell'implementazione del circuito dinamico.

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

Se ricevi errori relativi a neato non trovato da plot_circuit_layout(), assicurati di avere il graphviz pacchetto installato e disponibile nel tuo PATH. Se viene installato in una posizione non predefinita (ad esempio, utilizzando homebrew su MacOS ), potrebbe essere necessario aggiornare la variabile PATH di ambiente. Questo può essere fatto all'interno di questo notebook utilizzando quanto segue:

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:

Transpile utilizzando `MidCircuitMeasure`

MidCircuitMeasure è un'aggiunta alle operazioni di misurazione disponibili, calibrata specificamente per eseguire misurazioni a metà circuito. L'istruzione MidCircuitMeasure corrisponde measure_2 all'istruzione supportata dai backend. Si noti che non measure_2 è supportato su tutti i backend. Puoi usare service.backends(filters=lambda b: "measure_2" in b.supported_instructions) per trovare i backend che lo supportano. Qui mostriamo come traspilare il circuito in modo che le misurazioni a metà circuito definite nel circuito vengano eseguite utilizzando MidCircuitMeasure l'operazione, se il backend la supporta.

Di seguito riportiamo la durata measure_2 dell'istruzione e dell'istruzione 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
]

Traspilazione per circuiti unitari

Per stabilire un confronto equo tra i circuiti dinamici e la loro controparte unitaria, utilizziamo lo stesso insieme di qubit fisici utilizzati nei circuiti dinamici per i qubit di dati come layout per la transpilazione dei circuiti unitari.

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]

Visualizziamo la disposizione dei qubit dei circuiti unitari transpilati. I cerchi neri indicano i qubit fisici utilizzati per trasporre i circuiti unitari e i loro indici corrispondono agli indici dei qubit virtuali. Confrontando questo dato con il layout tracciato per i circuiti dinamici, possiamo confermare che i circuiti unitari utilizzano lo stesso insieme di qubit fisici dei qubit di dati nei circuiti dinamici.

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

Output:

Ora aggiungiamo la sequenza DD ai circuiti transpilati e costruiamo i PUB corrispondenti per l'invio dei lavori.

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]

Confronta la profondità del gate a due qubit dei circuiti unitari e dinamici

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

Il vantaggio principale del circuito basato sulla misurazione è che, quando si implementano più interazioni ZZ, gli strati CX possono essere parallelizzati e le misurazioni possono avvenire simultaneamente. Questo perché tutte le interazioni ZZ sono commutative, quindi il calcolo può essere eseguito con profondità di misurazione 1. Dopo aver trasposto i circuiti, osserviamo che l'approccio dinamico produce una profondità di due qubit significativamente più breve rispetto all'approccio unitario standard, con l'avvertenza che la misurazione aggiuntiva a metà circuito e il feedforward classico richiedono tempo e introducono fonti di errore proprie.

Passaggio 3: eseguire utilizzando Qiskit primitives

Modalità di test locale

Prima di inviare i lavori all'hardware, possiamo eseguire una piccola simulazione di prova del circuito dinamico utilizzando la modalità di test locale.

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

Simulazione MPS

Per circuiti di grandi dimensioni, possiamo utilizzare il simulatore matrix_product_state (MPS), che fornisce un risultato approssimativo al valore atteso in base alla dimensione del legame scelta. Successivamente utilizziamo i risultati della simulazione MPS come base di riferimento per confrontare i risultati ottenuti dall'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)

Una volta preparati i circuiti e gli osservabili, li eseguiamo sull'hardware utilizzando la primitiva Sampler.

Qui inviamo tre lavori per unitary_pubs, dynamic_pubs, e dynamic_pubs_dd. Ciascuno è un elenco di circuiti parametrizzati corrispondenti a nove diversi passaggi di Trotter con tre diversi parametri dell' $\theta$ .

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

Fase 4: Post-elaborazione e restituzione dei risultati nel formato classico desiderato

Una volta completati i processi, possiamo ricavare la durata del ciclo dai metadati dei risultati dei processi e visualizzare le informazioni relative alla pianificazione del ciclo. Per ulteriori informazioni sulla visualizzazione delle informazioni di pianificazione di un circuito, consultare questa pagina.

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

Tracciamo la durata dei circuiti unitari e dei circuiti dinamici. Dal grafico sottostante, possiamo vedere che, nonostante il tempo necessario per le misurazioni a metà circuito e le operazioni classiche, l'implementazione dinamica del circuito con measure_2 comporta durate del circuito paragonabili a quelle dell'implementazione 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 volta completati i lavori, recuperiamo i dati riportati di seguito e calcoliamo la magnetizzazione media stimata dagli osservabili observables_unitary o observables_dynamic che abbiamo costruito in precedenza.

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

Di seguito tracciamo la magnetizzazione di spin in funzione dei passi di Trotter a diversi valori di $\theta$ , corrispondenti a diverse intensità del campo magnetico locale. Tracciamo sia i risultati della simulazione MPS precalcolati per i circuiti ideali unitari, sia i risultati sperimentali ottenuti da quanto segue:

esecuzione dei circuiti unitari con DD
eseguendo i circuiti dinamici con DD e 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:

Quando confrontiamo i risultati sperimentali con la simulazione, vediamo che l'implementazione del circuito dinamico (linea tratteggiata con stelle) ha complessivamente prestazioni migliori rispetto all'implementazione unitaria standard (linea tratteggiata con cerchi). In sintesi, presentiamo i circuiti dinamici come soluzione per simulare i modelli di spin di Ising su un reticolo a nido d'ape, una topologia che non è nativa dell'hardware. La soluzione del circuito dinamico consente interazioni ZZ tra qubit che non sono vicini più prossimi, con una profondità di gate a due qubit più breve rispetto all'uso dei gate SWAP, al costo di introdurre qubit ancilla aggiuntivi e operazioni classiche di feedforward.

Riferimenti

[1] Quantum computing con Qiskit, di 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. e Johnson, B.R., 2024. arXiv preprint arXiv:2405.08810 (2024)

Questa pagina è stata utile?

Segnala un bug, un errore di battitura o richiedi contenuti su GitHub.

Simulazione dell'Hamiltoniano di Ising con circuiti dinamici

Requisiti

Impostazione di

Fase 1: mappare gli input classici su un circuito quantistico

Costruire circuiti unitari

Costruire un'implementazione dinamica del circuito

Disaccoppiamento dinamico (DD) e supporto per stretch la durata

Fase 2: Ottimizzazione del problema per l'esecuzione hardware

Transpilazione per circuiti dinamici

Transpile utilizzando MidCircuitMeasure

Traspilazione per circuiti unitari

Confronta la profondità del gate a due qubit dei circuiti unitari e dinamici

Passaggio 3: eseguire utilizzando Qiskit primitives

Modalità di test locale

Simulazione MPS

Fase 4: Post-elaborazione e restituzione dei risultati nel formato classico desiderato

Riferimenti

Disaccoppiamento dinamico (DD) e supporto per `stretch` la durata

Transpile utilizzando `MidCircuitMeasure`