Simularea Hamiltonianului Ising cu kick-uri folosind circuite dinamice

Estimare utilizare: 7,5 minute pe un procesor Heron r3. (NOTĂ: Aceasta este doar o estimare. Timpul tău de execuție poate varia.) Circuitele dinamice sunt circuite cu feedforward clasic — cu alte cuvinte, sunt măsurători efectuate în mijlocul circuitului, urmate de operații logice clasice care determină operații cuantice condiționate de rezultatul clasic. În acest tutorial, simulăm modelul Ising cu kick-uri pe o rețea hexagonală de spini și folosim circuite dinamice pentru a realiza interacțiuni dincolo de conectivitatea fizică a hardware-ului.

Modelul Ising a fost studiat extensiv în diverse domenii ale fizicii. Modelează spini care suferă interacțiuni Ising între site-urile rețelei, precum și kick-uri de la câmpul magnetic local de pe fiecare site. Evoluția temporală Trotterizată a spinilor considerată în acest tutorial, preluată din [1], este dată de următorul unitar:

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

Pentru a sonda dinamica spinilor, studiem magnetizarea medie a spinilor pe fiecare site în funcție de pașii Trotter. Prin urmare, construim următoarea observabilă:

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

Pentru a realiza interacțiunea ZZ între site-urile rețelei, prezentăm o soluție folosind funcționalitatea circuitului dinamic, care duce la o adâncime de doi qubiți semnificativ mai mică față de metoda standard de rutare cu porți SWAP. Pe de altă parte, operațiile de feedforward clasic din circuitele dinamice au, de regulă, timpi de execuție mai lungi decât porțile cuantice; prin urmare, circuitele dinamice au limitări și compromisuri. Prezentăm, de asemenea, o modalitate de a adăuga o secvență de decuplare dinamică pe qubiții în stare de repaus în timpul operației de feedforward clasic, folosind durata stretch.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele:

• Qiskit SDK v2.0 sau mai recent cu suport de vizualizare
• Qiskit Runtime v0.37 sau mai recent cu suport de vizualizare (pip install 'qiskit-ibm-runtime[visualization]')
• Biblioteca de grafuri Rustworkx (pip install rustworkx)
• Qiskit Aer (pip install qiskit-aer)

Configurare

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

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

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

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

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

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

Pasul 1: Maparea intrărilor clasice într-un Circuit cuantic

Începem prin definirea rețelei de simulat. Alegem să lucrăm cu rețeaua honeycomb (numită și hexagonală), care este un graf plan cu noduri de grad 3. Specificăm dimensiunea rețelei, parametrii relevanți ai circuitului de interes în dinamica Trotterizată. Simulăm evoluția temporală Trotterizată sub modelul Ising pentru trei valori diferite de $\theta$ ale câmpului magnetic 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

Să începem cu un exemplu de test mic:

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)

Vom folosi exemplul mic pentru ilustrare și simulare. Mai jos construim și un exemplu mare pentru a arăta că fluxul de lucru poate fi extins la dimensiuni mai mari.

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

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

num_qubits = 46

Construirea circuitelor unitare

Cu dimensiunea problemei și parametrii specificați, suntem acum gata să construim circuitul parametrizat care simulează evoluția temporală Trotterizată a $U(\theta)$ cu diferiți pași Trotter, specificați prin argumentul depth. Circuitul pe care îl construim are straturi alternante de porți Rx($\theta$) și porți Rzz. Porțile Rzz realizează interacțiunile ZZ între spinii cuplați, care vor fi plasate între fiecare site al rețelei specificat de argumentul layer_edges.

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

    return circuit

Vizualizează circuitul de test mic:

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)

Similar, construiește circuitele unitare ale exemplului mare la diferiți pași Trotter și observabila pentru a estima valoarea de așteptare.

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

Construirea implementării cu circuit dinamic

Această secțiune demonstrează implementarea principală cu circuit dinamic pentru a simula aceeași evoluție temporală Trotterizată. Rețineți că rețeaua honeycomb pe care dorim să o simulăm nu se potrivește cu rețeaua heavy lattice a qubiților hardware. O metodă directă de a mapa circuitul pe hardware constă în introducerea unei serii de operații SWAP pentru a aduce qubiții care interacționează unul lângă altul, în vederea realizării interacțiunii ZZ. Evidențiem aici o abordare alternativă folosind circuite dinamice ca soluție, care ilustrează că putem folosi combinația de calcul cuantic și calcul clasic în timp real în cadrul unui circuit în Qiskit pentru a realiza interacțiuni dincolo de vecinii cei mai apropiați.

În implementarea cu circuit dinamic, interacțiunea ZZ este implementată efectiv folosind qubiți ancilă, măsurători mid-circuit și feedforward. Pentru a înțelege acest lucru, notăm că rotațiile ZZ aplică un factor de fază $e^{i\theta}$ stării în funcție de paritatea sa. Pentru doi qubiți, stările bazei computaționale sunt $|00\rangle$, $|01\rangle$, $|10\rangle$ și $|11\rangle$. Poarta de rotație ZZ aplică un factor de fază stărilor $|01\rangle$ și $|10\rangle$ a căror paritate (numărul de unu din stare) este impară și lasă stările cu paritate pară neschimbate. Mai jos este descrisă modalitatea de implementare efectivă a interacțiunilor ZZ pe doi qubiți folosind circuite dinamice.

1. Calculează paritatea într-un qubit ancilă: în loc să aplicăm direct ZZ pe doi qubiți, introducem un al treilea qubit, qubit-ul ancilă, pentru a stoca informația de paritate a celor doi qubiți de date. Înlănțuim ancila cu fiecare qubit de date folosind porți CX de la qubit-ul de date la qubit-ul ancilă.

2. Aplică o rotație Z cu un singur qubit pe qubit-ul ancilă: aceasta deoarece ancila conține informația de paritate a celor doi qubiți de date, ceea ce implementează efectiv rotația ZZ pe qubiții de date.

3. Măsoară qubit-ul ancilă în baza X: acesta este pasul cheie care colapsează starea qubit-ului ancilă, iar rezultatul măsurătorii ne spune ce s-a întâmplat:

   • Măsoară 0: când se observă un rezultat 0, am aplicat corect o rotație $ZZ(\theta)$ pe qubiții noștri de date.

   • Măsoară 1: când se observă un rezultat 1, am aplicat $ZZ(\theta + \pi)$ în schimb.

4. Aplică poarta de corecție când se măsoară 1: dacă am măsurat 1, aplicăm porți Z pe qubiții de date pentru a „corecta" faza suplimentară de $\pi$.

Circuitul rezultat este următorul:

Când adoptăm această abordare pentru a simula o rețea honeycomb, circuitul rezultat se integrează perfect în hardware cu o rețea heavy-hex: toți qubiții de date se află pe site-urile de grad 3 ale rețelei, care formează o rețea hexagonală. Fiecare pereche de qubiți de date împarte un qubit ancilă aflat pe un site de grad 2. Mai jos, construim rețeaua de qubiți pentru implementarea cu circuit dinamic, introducând qubiți ancilă (marcați cu cercuri violet mai închis).

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)

Vizualizează rețeaua heavy-hex pentru qubiții de date și qubiții ancilă la o scară mică:

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

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

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

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

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

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

Mai jos, construim circuitul dinamic pentru evoluția temporală Trotterizată. Porțile RZZ sunt înlocuite cu implementarea cu circuit dinamic folosind pașii descriși mai sus.

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

Decuplare dinamică (DD) și suport pentru durata stretch

Un dezavantaj al utilizării implementării cu circuit dinamic pentru a realiza interacțiunea ZZ este că măsurătorile mid-circuit și operațiile de feedforward clasic durează, de obicei, mai mult decât porțile cuantice. Pentru a suprima decoerența qubiților în timpul de repaus aferent operațiilor clasice, am adăugat o secvență de decuplare dinamică (DD) după operația de măsurătoare pe qubiții ancilă și înainte de operația Z condiționată pe qubit-ul de date, înaintea instrucțiunii if_test.

Secvența DD este adăugată de funcția add_stretch_dd(), care folosește duratele stretch pentru a determina intervalele de timp dintre porțile DD. O durată stretch este o modalitate de a specifica o durată de timp extensibilă pentru operația delay, astfel încât durata întârzierii să poată crește pentru a umple timpul de repaus al qubitului. Variabilele de durată specificate prin stretch sunt rezolvate la momentul compilării în duratele dorite care satisfac o anumită constrângere. Aceasta este foarte utilă atunci când sincronizarea secvențelor DD este esențială pentru a obține o bună performanță de suprimare a erorilor. Pentru mai multe detalii despre tipul stretch, consultați documentația OpenQASM. În prezent, suportul pentru tipul stretch în Qiskit Runtime este experimental. Pentru detalii despre constrângerile de utilizare, consultați secțiunea de limitări din documentația stretch.

Folosind funcțiile definite mai sus, construim circuitele de evoluție temporală Trotterizată, cu și fără DD, și observabilele corespunzătoare. Începem prin vizualizarea circuitului dinamic al unui exemplu mic:

hex_rows_test = 1
hex_cols_test = 1

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

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

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

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

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

Similar, construiește circuitele dinamice pentru exemplul mare:

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)

Pasul 2: Optimizarea problemei pentru execuția pe hardware

Suntem acum gata să transpilăm circuitul pe hardware. Vom transpila atât implementarea unitară standard, cât și implementarea cu circuite dinamice pe hardware.

Pentru a transpila pe hardware, instanțiem mai întâi backend-ul. Dacă este disponibil, vom alege un backend pe care instrucțiunea MidCircuitMeasure (measure_2) este suportată.

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

Transpilarea pentru circuite dinamice

Mai întâi, transpilăm circuitele dinamice, cu și fără adăugarea secvenței DD. Pentru a ne asigura că folosim același set de qubiți fizici în toate circuitele, în vederea unor rezultate mai consistente, transpilăm circuitul o singură dată, apoi folosim layout-ul său pentru toate circuitele ulterioare, specificat prin initial_layout în pass manager. Construim apoi primitive unified blocs (PUB-uri) ca input pentru primitiva Sampler.

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

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

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

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

Putem vizualiza layout-ul qubiților din circuitul transpilat mai jos. Cercurile negre reprezintă qubiții de date și qubiții ancilă folosiți în implementarea cu circuite dinamice.

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

notă

Dacă primești erori referitoare la neato nu a fost găsit din plot_circuit_layout(), asigură-te că ai pachetul graphviz instalat și disponibil în PATH-ul tău. Dacă se instalează într-o locație non-implicită (de exemplu, folosind homebrew pe MacOS), poate fi necesar să actualizezi variabila de mediu PATH. Acest lucru poate fi făcut în interiorul acestui notebook folosind următoarele:

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

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

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

Transpilare folosind MidCircuitMeasure

MidCircuitMeasure este un plus față de operațiunile de măsurare disponibile, calibrat special pentru a efectua măsurători în mijlocul circuitului. Instrucțiunea MidCircuitMeasure corespunde instrucțiunii measure_2 suportate de backend-uri. Rețineți că measure_2 nu este suportată pe toate backend-urile. Poți folosi service.backends(filters=lambda b: "measure_2" in b.supported_instructions) pentru a găsi backend-urile care o suportă. Aici, arătăm cum să transpilăm circuitul astfel încât măsurătorile de mid-circuit definite în circuit să fie executate folosind operația MidCircuitMeasure, dacă backend-ul o suportă.

Mai jos, afișăm durata pentru instrucțiunea measure_2 și instrucțiunea standard measure.

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

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

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

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

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

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

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

        return dag

pm = PassManager(ConvertToMidCircuitMeasure(backend.target))

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

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

Transpilarea circuitelor unitare

Pentru a stabili o comparație corectă între circuitele dinamice și corespondentul lor unitar, folosim același set de qubiți fizici utilizați în circuitele dinamice pentru qubiții de date ca layout pentru transpilarea circuitelor unitare.

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]

Vizualizăm layout-ul de qubiți al circuitelor unitare transpilate. Cercurile negre indică qubiții fizici folosiți pentru transpilarea circuitelor unitare, iar indicii lor corespund indicilor qubiților virtuali. Comparând aceasta cu layout-ul reprezentat pentru circuitele dinamice, putem confirma că circuitele unitare folosesc același set de qubiți fizici ca qubiții de date din circuitele dinamice.

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

Acum adăugăm secvența DD la circuitele transpilate și construim PUB-urile corespunzătoare pentru trimiterea job-urilor.

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]

Compararea adâncimii Gate-urilor cu doi qubiți între circuitele unitare și cele dinamice

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

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

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

<matplotlib.legend.Legend at 0x374225760>

Principalul avantaj al circuitului bazat pe măsurători este că, atunci când se implementează mai multe interacțiuni ZZ, straturile CX pot fi paralelizate, iar măsurătorile pot avea loc simultan. Aceasta se datorează faptului că toate interacțiunile ZZ comutează, astfel că calculul poate fi efectuat cu adâncime de măsurare 1. După transpilarea circuitelor, observăm că abordarea cu circuit dinamic produce o adâncime pe doi qubiți semnificativ mai mică decât abordarea unitară standard, cu mențiunea că măsurătoarea intermediară și feedforward-ul clasic necesită ele însele timp și introduc propriile surse de erori.

Pasul 3: Execuție folosind primitivele Qiskit

Modul de testare locală

Înainte de a trimite joburile pe hardware, putem rula o mică simulare de test a circuitului dinamic folosind modul de testare locală.

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

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

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

array([ 0.16666667,  0.01855469, -0.13476562])

Simulare MPS

Pentru circuitele de mari dimensiuni, putem folosi simulatorul matrix_product_state (MPS), care oferă un rezultat aproximativ al valorii de așteptare în funcție de dimensiunea de legătură aleasă. Ulterior, vom folosi rezultatele simulării MPS ca referință pentru a compara cu rezultatele de pe 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)

Cu circuitele și observabilele pregătite, le executăm acum pe hardware folosind primitiva Sampler.

Trimitem trei joburi pentru unitary_pubs, dynamic_pubs și dynamic_pubs_dd. Fiecare este o listă de circuite parametrizate corespunzând la nouă pași Trotter diferiți cu trei parametri $\theta$ diferiți.

shots = 10000

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

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

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

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

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

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

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

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

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

Pasul 4: Post-procesare și returnarea rezultatelor în formatul clasic dorit

După ce joburile s-au finalizat, putem recupera durata circuitului din metadatele rezultatelor jobului și vizualiza informațiile despre planificarea circuitului. Pentru a citi mai multe despre vizualizarea informațiilor de planificare ale unui circuit, consultă această pagină.

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

Reprezentăm grafic duratele circuitelor unitare și ale circuitelor dinamice. Din graficul de mai jos, putem observa că, în ciuda timpului necesar pentru măsurătorile de la mijlocul circuitului și operațiile clasice, implementarea circuitului dinamic cu measure_2 produce durate ale circuitelor comparabile cu implementarea unitară.

# visualize circuit durations

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

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

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

<matplotlib.legend.Legend at 0x17f73c6e0>

După ce joburile s-au finalizat, recuperăm datele de mai jos și calculăm magnetizația medie estimată de observabilele observables_unitary sau observables_dynamic pe care le-am construit anterior.

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

Mai jos reprezentăm grafic magnetizația de spin în funcție de pașii Trotter la diferite valori $\theta$, corespunzând diferitelor intensități ale câmpului magnetic local. Reprezentăm atât rezultatele pre-calculate din simularea MPS pentru circuitele unitare ideale, cât și rezultatele experimentale din:

1. rularea circuitelor unitare cu DD
2. rularea circuitelor dinamice cu DD și MidCircuitMeasure

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

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

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

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

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

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

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

Când comparăm rezultatele experimentale cu simularea, observăm că implementarea cu circuit dinamic (linie punctată cu stele) are în general o performanță mai bună decât implementarea unitară standard (linie punctată cu cercuri). În concluzie, prezentăm circuitele dinamice ca o soluție pentru simularea modelelor de spin Ising pe o rețea în fagure, o topologie care nu este nativă pentru hardware. Soluția cu circuit dinamic permite interacțiuni ZZ între Qubiți care nu sunt vecini imediați, cu o adâncime de porți cu doi Qubiți mai mică decât în cazul folosirii porților SWAP, cu prețul introducerii unor Qubiți ancila suplimentari și a operațiilor clasice de feedforward.

Referințe

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