Simulace kicked Ising Hamiltoniánu s dynamickými obvody

Odhad využití: 7,5 minut na procesoru Heron r3. (POZNÁMKA: Jedná se pouze o odhad. Skutečný čas může být odlišný.) Dynamické obvody jsou obvody s klasickým feedforwardem — jinými slovy jde o mid-circuit měření následovaná klasickými logickými operacemi, které určují kvantové operace podmíněné klasickým výstupem. V tomto tutoriálu simulujeme kicked Ising model na hexagonální mřížce spinů a používáme dynamické obvody k realizaci interakcí přesahujících fyzické propojení hardwaru.

Isingův model byl rozsáhle studován napříč oblastmi fyziky. Modeluje spiny, které procházejí Isingovými interakcemi mezi uzly mřížky, a také „kopy" od lokálního magnetického pole na každém uzlu. Trotterizovaný časový vývoj spinů uvažovaný v tomto tutoriálu, převzatý z [1], je dán následujícím unitárním operátorem:

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

Pro zkoumání spinové dynamiky studujeme průměrnou magnetizaci spinů na každém uzlu jako funkci Trotterových kroků. Konstruujeme tedy následující pozorovatelnou veličinu:

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

Pro realizaci ZZ interakce mezi uzly mřížky představujeme řešení využívající vlastnost dynamických obvodů, která vede k výrazně kratší dvouqubitové hloubce ve srovnání se standardní metodou směrování pomocí SWAP hradel. Na druhou stranu klasické feedforward operace v dynamických obvodech mají typicky delší dobu provádění než kvantová hradla; dynamické obvody tedy mají svá omezení a kompromisy. Také ukážeme způsob, jak přidat sekvenci dynamického oddělování (dynamical decoupling) na nečinné qubity během klasické feedforward operace pomocí trvání stretch.

Požadavky

Než začneš s tímto tutoriálem, ujisti se, že máš nainstalované následující:

• Qiskit SDK v2.0 nebo novější s podporou vizualizace
• Qiskit Runtime v0.37 nebo novější s podporou vizualizace (pip install 'qiskit-ibm-runtime[visualization]')
• Knihovna grafů Rustworkx (pip install rustworkx)
• Qiskit Aer (pip install qiskit-aer)

Nastavení

# 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: Mapování klasických vstupů na kvantový obvod

Začneme definováním mřížky k simulaci. Zvolíme honeycomb (neboli hexagonální) mřížku, což je planární graf s uzly stupně 3. Zde specifikujeme velikost mřížky a relevantní parametry obvodu v Trotterizované dynamice. Simulujeme Trotterizovaný časový vývoj pod Isingovým modelem pro tři různé hodnoty $\theta$ lokálního magnetického pole.

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

Začněme s malým testovacím příkladem:

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)

Pro ilustraci a simulaci budeme používat malý příklad. Níže také konstruujeme velký příklad, abychom ukázali, že pracovní postup lze rozšířit na větší velikosti.

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

Sestavení unitárních obvodů

S definovanou velikostí problému a parametry jsme nyní připraveni sestavit parametrizovaný obvod, který simuluje Trotterizovaný časový vývoj $U(\theta)$ s různými Trotterovými kroky, specifikovanými argumentem depth. Obvod, který konstruujeme, má střídající se vrstvy hradel Rx($\theta$) a Rzz. Hradla Rzz realizují ZZ interakce mezi spojenými spiny, která budou umístěna mezi každý uzel mřížky specifikovaný argumentem 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

Vizualizace malého testovacího obvodu:

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)

Podobně konstruujeme unitární obvody pro velký příklad při různých Trotterových krocích a pozorovatelnou veličinu pro odhad střední hodnoty.

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

Sestavení implementace dynamického obvodu

Tato sekce demonstruje hlavní implementaci dynamického obvodu pro simulaci stejného Trotterizovaného časového vývoje. Všimni si, že honeycomb mřížka, kterou chceme simulovat, neodpovídá heavy lattice qubitů hardwaru. Jedním přímočarým způsobem, jak mapovat obvod na hardware, je zavedení série SWAP operací, aby se interagující qubity dostaly vedle sebe a umožnily realizaci ZZ interakce. Zde zdůrazňujeme alternativní přístup pomocí dynamických obvodů jako řešení, které ilustruje, že v Qiskitu můžeme v rámci obvodu kombinovat kvantové a real-time klasické výpočty k realizaci interakcí přesahujících nejbližší sousedy.

V implementaci dynamického obvodu je ZZ interakce efektivně implementována pomocí ancilla qubitů, mid-circuit měření a feedforwardu. Abys to lépe pochopil/a, vezmi v úvahu, že ZZ rotace aplikují fázový faktor $e^{i\theta}$ na stav na základě jeho parity. Pro dva qubity jsou stavové báze $|00\rangle$, $|01\rangle$, $|10\rangle$ a $|11\rangle$. ZZ rotační hradlo aplikuje fázový faktor na stavy $|01\rangle$ a $|10\rangle$, jejichž parita (počet jedniček ve stavu) je lichá, a ponechává stavy s sudou paritou nezměněné. Níže je popsáno, jak lze efektivně implementovat ZZ interakce na dvou qubitech pomocí dynamických obvodů.

1. Výpočet parity do ancilla qubitu: místo přímého aplikování ZZ na dva qubity zavedeme třetí qubit — ancilla qubit — pro uložení informace o paritě dvou datových qubitů. Ancilla provážeme s každým datovým qubitem pomocí CX hradel z datového qubitu na ancilla qubit.

2. Aplikování jednoqubitové Z rotace na ancilla qubit: je to proto, že ancilla obsahuje informaci o paritě dvou datových qubitů, což efektivně implementuje ZZ rotaci na datových qubitech.

3. Měření ancilla qubitu v X bázi: toto je klíčový krok, který kolabuje stav ancilla qubitu, a výsledek měření nám říká, co se stalo:

   • Měření 0: když je pozorován výsledek 0, správně jsme aplikovali rotaci $ZZ(\theta)$ na datové qubity.

   • Měření 1: když je pozorován výsledek 1, aplikovali jsme místo toho $ZZ(\theta + \pi)$.

4. Aplikování korekčního hradla při měření 1: pokud jsme naměřili 1, aplikujeme Z hradla na datové qubity, abychom „opravili" extra $\pi$ fázi.

Výsledný obvod je následující:

Když přijmeme tento přístup pro simulaci honeycomb mřížky, výsledný obvod se dokonale zabuduje do hardwaru s heavy-hex mřížkou: všechny datové qubity se nacházejí na uzlech stupně 3 mřížky, které tvoří hexagonální mřížku. Každý pár datových qubitů sdílí ancilla qubit umístěný na uzlu stupně 2. Níže konstruujeme qubitovou mřížku pro implementaci dynamického obvodu a zavedeme ancilla qubity (zobrazené jako tmavší fialové kruhy).

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)

Vizualizace heavy-hex mřížky pro datové qubity a ancilla qubity v malém měřítku:

(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

Níže konstruujeme dynamický obvod pro Trotterizovaný časový vývoj. Hradla RZZ jsou nahrazena implementací dynamického obvodu pomocí kroků popsaných výše.

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

Dynamické oddělování (DD) a podpora trvání stretch

Jednou nevýhodou použití implementace dynamického obvodu pro realizaci ZZ interakce je, že mid-circuit měření a klasické feedforward operace typicky trvají déle než kvantová hradla. Aby se potlačila dekoherence qubitů během doby nečinnosti pro klasické operace, přidali jsme sekvenci dynamického oddělování (DD) po operaci měření na ancilla qubitech a před podmíněnou Z operací na datovém qubitu, před příkazem if_test.

Sekvence DD je přidána funkcí add_stretch_dd(), která využívá trvání stretch k určení časových intervalů mezi DD hradly. Trvání stretch je způsob, jak specifikovat roztažitelnou časovou délku pro operaci delay, takže délka zpoždění může narůst, aby zaplnila dobu nečinnosti qubitu. Proměnné trvání specifikované pomocí stretch jsou při kompilaci přeloženy na požadované délky splňující určité omezení. To je velmi užitečné, když je načasování DD sekvencí zásadní pro dobrou výkonnost potlačení chyb. Více podrobností o typu stretch najdeš v dokumentaci OpenQASM. Podpora typu stretch v Qiskit Runtime je v současné době experimentální. Podrobnosti o omezeních jejího použití najdeš v sekci omezení dokumentace stretch.

Pomocí výše definovaných funkcí sestavíme obvody Trotterizovaného časového vývoje s DD i bez DD a odpovídající pozorovatelné veličiny. Začneme vizualizací dynamického obvodu malého příkladu:

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)

Podobně konstruujeme dynamické obvody pro velký příklad:

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: Optimalizace problému pro spuštění na hardwaru

Nyní jsme připraveni transpilovat Circuit na hardware. Transpilujeme jak standardní unitární implementaci, tak implementaci s dynamickými obvody.

Abychom mohli transpilovat na hardware, nejprve vytvoříme instanci Backendu. Pokud je to možné, vybereme Backend, který podporuje instrukci 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,
    )

Transpilace pro dynamické obvody

Nejprve transpilujeme dynamické obvody, a to s DD sekvencí i bez ní. Abychom zajistili použití stejné sady fyzických Qubitů ve všech obvodech pro konzistentnější výsledky, nejprve transpilujeme Circuit jednou a poté použijeme jeho rozvržení pro všechny následující obvody, specifikované pomocí initial_layout v pass manageru. Poté sestavíme primitive unified bloky (PUBs) jako vstup pro primitivní rozhraní 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]

Níže můžeš vizualizovat rozvržení Qubitů transpilovaného obvodu. Černé kruhy znázorňují datové Qubity a ancilla Qubity použité v implementaci s dynamickými obvody.

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

poznámka

Pokud se zobrazí chyby o nenalezeném neato z plot_circuit_layout(), ujisti se, že máš nainstalovaný balíček graphviz a že je dostupný ve tvém PATH. Pokud se nainstaluje do nestandardního umístění (například při použití homebrew na MacOS), možná budeš muset aktualizovat proměnnou prostředí PATH. To lze provést přímo v tomto notebooku takto:

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)

Transpilace pomocí MidCircuitMeasure

MidCircuitMeasure je rozšíření dostupných operací měření, kalibrované speciálně pro provádění měření uprostřed obvodu. Instrukce MidCircuitMeasure se mapuje na instrukci measure_2 podporovanou Backendy. Upozorňujeme, že measure_2 není podporováno na všech Backendech. Pomocí service.backends(filters=lambda b: "measure_2" in b.supported_instructions) lze najít Backendy, které ho podporují. Zde ukážeme, jak transpilovat Circuit tak, aby mid-circuit měření definovaná v obvodu byla prováděna pomocí operace MidCircuitMeasure, pokud to Backend podporuje.

Níže vytiskneme dobu trvání instrukce measure_2 a standardní instrukce 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
]

Transpilace pro unitární obvody

Abychom zajistili spravedlivé srovnání dynamických obvodů s jejich unitárním protějškem, použijeme stejnou sadu fyzických Qubitů, které byly použity v dynamických obvodech pro datové Qubity, jako rozvržení pro transpilaci unitárních obvodů.

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]

Vizualizujeme rozvržení Qubitů transpilovaných unitárních obvodů. Černé kruhy označují fyzické Qubity použité pro transpilaci unitárních obvodů a jejich indexy odpovídají indexům virtuálních Qubitů. Porovnáním s rozvržením vykresleným pro dynamické obvody lze potvrdit, že unitární obvody používají stejnou sadu fyzických Qubitů jako datové Qubity v dynamických obvodech.

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

Nyní přidáme DD sekvenci k transpilovaným obvodům a sestavíme odpovídající PUBs pro odeslání úlohy.

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]

Porovnání hloubky dvoukubitových Gate unitárních a dynamických obvodů

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

Hlavní výhodou obvodu založeného na měření je, že při implementaci více interakcí ZZ lze vrstvy CX paralelizovat a měření mohou probíhat současně. Je tomu tak proto, že všechny interakce ZZ komutují, a výpočet lze tedy provést s hloubkou měření 1. Po transpilaci obvodů pozorujeme, že přístup s dynamickými obvody dosahuje výrazně kratší hloubky dvoukubitových Gate než standardní unitární přístup, s tím upozorněním, že samotné mid-circuit měření a klasická zpětná vazba (feedforward) zabírají čas a přinášejí vlastní zdroje chyb.

Krok 3: Spuštění pomocí primitiv Qiskit

Lokální testovací režim

Před odesláním úloh na hardware můžeme spustit malou testovací simulaci dynamického Circuit pomocí lokálního testovacího režimu.

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

Simulace MPS

Pro velké Circuit lze použít simulátor matrix_product_state (MPS), který poskytuje přibližný výsledek střední hodnoty podle zvolené dimenze vazby. Výsledky MPS simulace později použijeme jako základní linii pro porovnání s výsledky z 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)

Nyní, když jsou Circuit a pozorovatelné připraveny, spustíme je na hardware pomocí primitiva Sampler.

Zde odesíláme tři úlohy pro unitary_pubs, dynamic_pubs a dynamic_pubs_dd. Každá je seznam parametrizovaných Circuit odpovídajících devíti různým Trotterovým krokům se třemi různými parametry $\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: Zpracování výsledků a jejich vrácení v požadovaném klasickém formátu

Po dokončení úloh můžeme z metadat výsledků úloh načíst dobu trvání Circuit a vizualizovat informace o harmonogramu Circuit. Více o vizualizaci informací o plánování Circuit najdeš na této stránce.

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

Zobrazujeme doby trvání Circuit pro unitární Circuit a dynamické Circuit. Z níže uvedeného grafu vidíme, že navzdory času potřebnému pro mid-circuit měření a klasické operace vede implementace dynamického Circuit s measure_2 k podobným dobám trvání Circuit jako unitární implementace.

# 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 dokončení úloh načteme níže uvedená data a vypočítáme průměrnou magnetizaci odhadnutou pozorovatelnými observables_unitary nebo observables_dynamic, které jsme sestavili dříve.

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

Níže zobrazujeme spinovou magnetizaci jako funkci Trotterových kroků při různých hodnotách $\theta$, odpovídajících různým intenzitám lokálního magnetického pole. Zobrazujeme předem vypočtené výsledky MPS simulace pro unitární ideální Circuit spolu s experimentálními výsledky z:

1. spuštění unitárních Circuit s DD
2. spuštění dynamických Circuit s DD a 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()

Při porovnání experimentálních výsledků se simulací vidíme, že implementace dynamického Circuit (tečkovaná čára s hvězdičkami) má celkově lepší výkon než standardní unitární implementace (tečkovaná čára s kroužky). Shrnuto, představujeme dynamické Circuit jako řešení pro simulaci Isingových spinových modelů na mřížce ve tvaru plástve, topologie, která není nativní pro hardware. Řešení s dynamickým Circuit umožňuje ZZ interakce mezi Qubit, které nejsou nejbližšími sousedy, s kratší hloubkou dvouQubitových Gate než při použití SWAP Gate, za cenu zavedení dalších ancilla Qubit a operací klasického zpětného přenosu.

Reference

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