Symulacja hamiltoniana Isinga z kopnięciami przy użyciu obwodów dynamicznych

Szacowane użycie: 7,5 minuty na procesorze Heron r3. (UWAGA: To jedynie szacunek. Rzeczywisty czas wykonania może się różnić.) Obwody dynamiczne to obwody z klasycznym sprzężeniem zwrotnym — innymi słowy, są to pomiary w trakcie obwodu, po których następują klasyczne operacje logiczne warunkujące operacje kwantowe na podstawie klasycznego wyniku. W tym samouczku symulujemy model Isinga z kopnięciami na hexagonalnej sieci spinów i używamy obwodów dynamicznych do realizacji oddziaływań wykraczających poza fizyczną łączność sprzętu.

Model Isinga był szeroko badany w różnych dziedzinach fizyki. Opisuje spiny podlegające oddziaływaniom Isinga między węzłami sieci, a także kopnięcia od lokalnego pola magnetycznego w każdym węźle. Trotteryzowana ewolucja czasowa spinów rozważana w tym samouczku, zaczerpnięta z [1], dana jest przez następującą unitarność:

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

Aby zbadać dynamikę spinów, badamy średnie namagnesowanie spinów w każdym węźle w funkcji kroków Trottera. W tym celu konstruujemy następującą obserwowalną:

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

Aby zrealizować oddziaływanie ZZ między węzłami sieci, przedstawiamy rozwiązanie wykorzystujące funkcję obwodów dynamicznych, prowadzące do znacznie krótszej głębokości dwuqubitowej w porównaniu ze standardową metodą trasowania z bramkami SWAP. Z drugiej strony klasyczne operacje sprzężenia zwrotnego w obwodach dynamicznych mają zazwyczaj dłuższy czas wykonania niż bramki kwantowe; dlatego obwody dynamiczne mają swoje ograniczenia i kompromisy. Pokazujemy również sposób dodawania sekwencji dynamicznego rozsprzęgania na bezczynnych qubitach podczas operacji klasycznego sprzężenia zwrotnego przy użyciu czasu trwania stretch.

Wymagania

Przed rozpoczęciem tego samouczka upewnij się, że masz zainstalowane następujące elementy:

• Qiskit SDK v2.0 lub nowszy z obsługą wizualizacji
• Qiskit Runtime v0.37 lub nowszy z obsługą wizualizacji (pip install 'qiskit-ibm-runtime[visualization]')
• Biblioteka grafów Rustworkx (pip install rustworkx)
• Qiskit Aer (pip install qiskit-aer)

Konfiguracja

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

Krok 1: Odwzorowanie klasycznych danych wejściowych na obwód kwantowy

Zaczynamy od zdefiniowania sieci do symulacji. Wybieramy sieć plastra miodu (zwaną też hexagonalną), która jest planarnym grafem o węzłach stopnia 3. Tutaj określamy rozmiar sieci oraz istotne parametry obwodu dla trotteryzowanej dynamiki. Symulujemy trotteryzowaną ewolucję czasową w modelu Isinga dla trzech różnych wartości $\theta$ lokalnego pola magnetycznego.

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

Zacznijmy od małego przykładu testowego:

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)

Użyjemy małego przykładu do ilustracji i symulacji. Poniżej konstruujemy też duży przykład, aby pokazać, że przepływ pracy można rozszerzyć do większych rozmiarów.

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

Budowanie obwodów unitarnych

Po określeniu rozmiaru problemu i parametrów jesteśmy gotowi do zbudowania sparametryzowanego obwodu symulującego trotteryzowaną ewolucję czasową $U(\theta)$ przy różnych krokach Trottera, określonych przez argument depth. Konstruowany obwód ma naprzemienne warstwy bramek Rx($\theta$) i bramek Rzz. Bramki Rzz realizują oddziaływania ZZ między sprzężonymi spinami, które będą umieszczone między każdym węzłem sieci określonym przez argument 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

Zwizualizuj mały obwód testowy:

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)

Podobnie, skonstruuj obwody unitarne dużego przykładu przy różnych krokach Trottera oraz obserwowalną do szacowania wartości oczekiwanej.

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

Budowanie implementacji z obwodem dynamicznym

Ta sekcja demonstruje główną implementację z obwodem dynamicznym służącą do symulacji tej samej trotteryzowanej ewolucji czasowej. Zauważ, że sieć plastra miodu, którą chcemy symulować, nie odpowiada ciężkiej sieci qubitów sprzętowych. Jednym prostym sposobem odwzorowania obwodu na sprzęt jest wprowadzenie serii operacji SWAP w celu zbliżenia oddziałujących qubitów, aby zrealizować oddziaływanie ZZ. Tutaj wyróżniamy alternatywne podejście z użyciem obwodów dynamicznych, które ilustruje, że możemy połączyć obliczenia kwantowe z klasycznymi obliczeniami w czasie rzeczywistym wewnątrz obwodu w Qiskit, aby realizować oddziaływania wykraczające poza sąsiedztwo.

W implementacji z obwodem dynamicznym oddziaływanie ZZ jest efektywnie realizowane przy użyciu qubitów pomocniczych, pomiaru w trakcie obwodu i sprzężenia zwrotnego. Aby to zrozumieć, zauważmy, że rotacje ZZ nadają czynnik fazowy $e^{i\theta}$ stanowi w zależności od jego parzystości. Dla dwóch qubitów stany bazowe obliczeniowe to $|00\rangle$, $|01\rangle$, $|10\rangle$ i $|11\rangle$. Bramka rotacji ZZ nadaje czynnik fazowy stanom $|01\rangle$ i $|10\rangle$, których parzystość (liczba jedynek w stanie) jest nieparzysta, i pozostawia stany parzyste bez zmian. Poniżej opisujemy, jak można efektywnie zrealizować oddziaływania ZZ na dwóch qubitach przy użyciu obwodów dynamicznych.

1. Oblicz parzystość do qubitu pomocniczego: zamiast bezpośrednio stosować ZZ do dwóch qubitów, wprowadzamy trzeci qubit — qubit pomocniczy — do przechowywania informacji o parzystości dwóch qubitów danych. Splątujemy qubit pomocniczy z każdym qubitem danych za pomocą bramek CX od qubitu danych do qubitu pomocniczego.

2. Zastosuj jednoQubitową rotację Z do qubitu pomocniczego: ponieważ qubit pomocniczy zawiera informację o parzystości dwóch qubitów danych, efektywnie realizuje to rotację ZZ na qubitach danych.

3. Zmierz qubit pomocniczy w bazie X: to kluczowy krok, który powoduje kolaps stanu qubitu pomocniczego, a wynik pomiaru mówi nam, co się stało:

   • Wynik 0: gdy obserwowany jest wynik 0, poprawnie zastosowaliśmy rotację $ZZ(\theta)$ do naszych qubitów danych.

   • Wynik 1: gdy obserwowany jest wynik 1, zastosowaliśmy $ZZ(\theta + \pi)$.

4. Zastosuj bramkę korekcyjną przy wyniku 1: jeśli zmierzono 1, stosujemy bramki Z do qubitów danych, aby „naprawić" dodatkową fazę $\pi$.

Wynikowy obwód jest następujący:

Gdy przyjmiemy to podejście do symulacji sieci plastra miodu, wynikowy obwód doskonale wbudowuje się w sprzęt z siecią heavy-hex: wszystkie qubity danych leżą na węzłach stopnia 3 sieci, tworząc sieć hexagonalną. Każda para qubitów danych dzieli qubit pomocniczy leżący na węźle stopnia 2. Poniżej konstruujemy sieć qubitów dla implementacji z obwodem dynamicznym, wprowadzając qubity pomocnicze (pokazane ciemnofioletowymi okręgami).

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)

Zwizualizuj sieć heavy-hex z qubitami danych i qubitami pomocniczymi w małej skali:

(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

Poniżej konstruujemy obwód dynamiczny dla trotteryzowanej ewolucji czasowej. Bramki RZZ są zastępowane implementacją z obwodem dynamicznym zgodnie z opisanymi powyżej krokami.

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

Dynamiczne rozsprzęganie (DD) i obsługa czasu trwania stretch

Jednym zastrzeżeniem dotyczącym użycia implementacji z obwodem dynamicznym do realizacji oddziaływania ZZ jest to, że pomiar w trakcie obwodu i klasyczne operacje sprzężenia zwrotnego zazwyczaj wymagają dłuższego czasu wykonania niż bramki kwantowe. Aby tłumić dekoherencję qubitów podczas czasu bezczynności potrzebnego na przeprowadzenie operacji klasycznych, dodaliśmy sekwencję dynamicznego rozsprzęgania (DD) po operacji pomiaru qubitów pomocniczych, a przed warunkową operacją Z na qubicie danych, przed instrukcją if_test.

Sekwencja DD jest dodawana przez funkcję add_stretch_dd(), która używa czasu trwania stretch do określenia przedziałów czasowych między bramkami DD. Czas trwania stretch to sposób na określenie rozciągliwego czasu trwania dla operacji delay, tak aby czas opóźnienia mógł rosnąć i wypełniać czas bezczynności qubitu. Zmienne czasu trwania określone przez stretch są rozwiązywane w czasie kompilacji do pożądanych czasów trwania spełniających określone ograniczenie. Jest to bardzo przydatne, gdy synchronizacja sekwencji DD jest istotna dla dobrej wydajności tłumienia błędów. Więcej szczegółów na temat typu stretch znajdziesz w dokumentacji OpenQASM. Obecnie obsługa typu stretch w Qiskit Runtime jest eksperymentalna. Szczegóły dotyczące ograniczeń użycia znajdziesz w sekcji ograniczeń dokumentacji stretch.

Używając zdefiniowanych powyżej funkcji, budujemy obwody trotteryzowanej ewolucji czasowej z DD i bez DD oraz odpowiadające im obserwowalne. Zaczniemy od wizualizacji obwodu dynamicznego na małym przykładzie:

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)

Podobnie, skonstruuj obwody dynamiczne dla dużego przykładu:

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)

Krok 2: Optymalizacja problemu pod kątem wykonania na sprzęcie

Jesteśmy teraz gotowi do transpilacji układu na sprzęt. Przetranspilujemy zarówno standardową implementację unitarną, jak i implementację opartą na układach dynamicznych.

Aby dokonać transpilacji na sprzęt, najpierw tworzymy instancję backendu. Jeśli jest dostępny, wybierzemy backend, który obsługuje instrukcję 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,
    )

Transpilacja dla układów dynamicznych

Najpierw transpilujemy układy dynamiczne — z sekwencją DD i bez niej. Aby zapewnić spójność wyników, używamy tego samego zestawu fizycznych Qubitów we wszystkich układach. W tym celu transpilujemy układ jeden raz, a następnie używamy jego rozkładu dla wszystkich kolejnych układów, podając go przez parametr initial_layout w menedżerze przebiegów. Następnie konstruujemy primitive unified blocs (PUBs) jako dane wejściowe dla prymitywu 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]

Poniżej możemy zwizualizować rozkład Qubitów w przetranspilowanym układzie. Czarne kółka oznaczają Qubity danych oraz Qubity pomocnicze (ancilla) użyte w implementacji z układami dynamicznymi.

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

uwaga

Jeśli pojawią się błędy związane z nieznalezieniem neato w plot_circuit_layout(), upewnij się, że pakiet graphviz jest zainstalowany i dostępny w zmiennej PATH. Jeśli zainstaluje się w niestandardowej lokalizacji (np. przy użyciu homebrew na MacOS), może być konieczna ręczna aktualizacja zmiennej środowiskowej PATH. Można to zrobić wewnątrz tego notebooka w następujący sposób:

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)

Transpilacja z użyciem MidCircuitMeasure

MidCircuitMeasure to rozszerzenie dostępnych operacji pomiarowych, skalibrowane specjalnie do wykonywania pomiarów w trakcie układu. Instrukcja MidCircuitMeasure mapuje się na instrukcję measure_2 obsługiwaną przez backendy. Należy zauważyć, że measure_2 nie jest obsługiwane przez wszystkie backendy. Możesz użyć service.backends(filters=lambda b: "measure_2" in b.supported_instructions), aby znaleźć backendy, które go obsługują. Poniżej pokazujemy, jak przetranspilować układ tak, aby pomiary w trakcie jego trwania były wykonywane przy użyciu operacji MidCircuitMeasure, jeśli backend to obsługuje.

Poniżej wyświetlamy czas trwania instrukcji measure_2 oraz standardowej instrukcji 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
]

Transpilacja dla układów unitarnych

Aby zapewnić uczciwe porównanie między układami dynamicznymi a ich unitarnym odpowiednikiem, używamy tego samego zestawu fizycznych Qubitów co w układach dynamicznych (dla Qubitów danych) jako rozkładu przy transpilacji układów unitarnych.

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]

Wizualizujemy rozkład Qubitów przetranspilowanych układów unitarnych. Czarne kółka wskazują fizyczne Qubity użyte do transpilacji układów unitarnych, a ich indeksy odpowiadają indeksom wirtualnych Qubitów. Porównując to z rozkładem przedstawionym dla układów dynamicznych, możemy potwierdzić, że układy unitarne używają tego samego zestawu fizycznych Qubitów co Qubity danych w układach dynamicznych.

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

Teraz dodajemy sekwencję DD do przetranspilowanych układów i konstruujemy odpowiadające im PUBs do przesłania zadania.

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]

Porównanie głębokości dwuqubitowej układów unitarnych i dynamicznych

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

Główną zaletą układu opartego na pomiarach jest to, że przy implementowaniu wielu interakcji ZZ warstwy CX mogą być zrównoleglone, a pomiary mogą odbywać się jednocześnie. Wynika to z faktu, że wszystkie interakcje ZZ komutują ze sobą, dzięki czemu obliczenia można wykonać przy głębokości pomiarowej równej 1. Po przetranspilowaniu układów obserwujemy, że podejście z dynamicznymi układami daje znacznie mniejszą głębokość dwuqubitową niż standardowe podejście unitarne, choć należy pamiętać, że dodatkowe pomiary w trakcie układu oraz klasyczne sprzężenie zwrotne same zajmują czas i wprowadzają własne źródła błędów.

Krok 3: Wykonanie przy użyciu prymitywów Qiskit

Lokalny tryb testowy

Przed wysłaniem zadań na sprzęt możemy przeprowadzić małą symulację testową obwodu dynamicznego przy użyciu lokalnego trybu testowego.

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

Symulacja MPS

W przypadku dużych obwodów można użyć symulatora matrix_product_state (MPS), który dostarcza przybliżony wynik wartości oczekiwanej zgodnie z wybranym wymiarem więzów. Wyniki symulacji MPS wykorzystujemy później jako punkt odniesienia do porównania z wynikami ze sprzętu.

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

Mając przygotowane obwody i obserwable, wykonujemy je teraz na sprzęcie przy użyciu prymitywu Sampler.

Przesyłamy trzy zadania dla unitary_pubs, dynamic_pubs oraz dynamic_pubs_dd. Każde z nich to lista sparametryzowanych obwodów odpowiadających dziewięciu różnym krokom Trottera z trzema różnymi parametrami $\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()}")

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

Krok 4: Przetwarzanie końcowe i zwracanie wyników w żądanym formacie klasycznym

Po zakończeniu zadań możemy pobrać czas trwania obwodu z metadanych wyników zadania i zwizualizować informacje o harmonogramie obwodu. Aby dowiedzieć się więcej o wizualizacji informacji o harmonogramowaniu obwodu, zajrzyj na tę stronę.

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

Poniżej wykreślamy czasy trwania obwodów unitarnych i dynamicznych. Z wykresu wynika, że pomimo czasu potrzebnego na pomiary w trakcie obwodu i operacje klasyczne, implementacja obwodu dynamicznego z measure_2 daje porównywalne czasy trwania obwodu jak implementacja unitarna.

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

Po zakończeniu zadań pobieramy poniższe dane i obliczamy średnie namagnesowanie oszacowane przez obserwable observables_unitary lub observables_dynamic, które zbudowaliśmy wcześniej.

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

Poniżej wykreślamy namagnesowanie spinowe w funkcji kroków Trottera dla różnych wartości $\theta$, odpowiadających różnym natężeniom lokalnego pola magnetycznego. Przedstawiamy zarówno wstępnie obliczone wyniki symulacji MPS dla unitarnych, idealnych obwodów, jak i wyniki eksperymentalne z:

1. uruchamiania obwodów unitarnych z DD
2. uruchamiania obwodów dynamicznych z 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()

Porównując wyniki eksperymentalne z symulacją, widzimy, że implementacja obwodu dynamicznego (linia przerywana z gwiazdkami) ogólnie osiąga lepsze wyniki niż standardowa implementacja unitarna (linia przerywana z kółkami). Podsumowując, przedstawiamy obwody dynamiczne jako rozwiązanie do symulacji modeli spinowych Isinga na siatce plastra miodu — topologii, która nie jest natywna dla sprzętu. Rozwiązanie z obwodami dynamicznymi umożliwia interakcje ZZ między Qubitami, które nie są najbliższymi sąsiadami, przy krótszej głębokości bramek dwu-Qubitowych niż przy użyciu bramek SWAP, kosztem wprowadzenia dodatkowych Qubitów pomocniczych i operacji klasycznego sprzężenia zwrotnego.

Odniesienia

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