Transverzální Isingův model s řízením výkonu od Q-CTRL

Odhadovaná spotřeba: 2 minuty na procesoru Heron r2. (POZNÁMKA: Jde pouze o odhad. Skutečná doba běhu se může lišit.)

Pozadí

Transverzální Isingův model (TFIM) je důležitý pro studium kvantového magnetismu a fázových přechodů. Popisuje sadu spinů uspořádaných na mřížce, kde každý spin interaguje se svými sousedy a zároveň je ovlivňován externím magnetickým polem, které vyvolává kvantové fluktuace.

Běžným přístupem k simulaci tohoto modelu je použití Trotterovy dekompozice k aproximaci operátoru časového vývoje, přičemž se konstruují obvody, které střídají jednoqubitové rotace a dvoququbitové provázací interakce. Tato simulace na reálném hardwaru je však náročná kvůli šumu a dekoherenci, což vede k odchylkám od skutečné dynamiky. K překonání tohoto problému používáme nástroje Q-CTRL Fire Opal pro potlačení chyb a řízení výkonu, nabízené jako Qiskit Function (viz dokumentace Fire Opal). Fire Opal automaticky optimalizuje provádění obvodů pomocí dynamického oddělování, pokročilého rozmístění, routingu a dalších technik potlačení šumu, přičemž vše směřuje ke snížení šumu. Díky těmto vylepšením výsledky na hardwaru lépe odpovídají bezšumovým simulacím, a tak můžeme studovat magnetizační dynamiku TFIM s vyšší věrností.

V tomto tutoriálu budeme:

• Sestavovat TFIM Hamiltonián na grafu propojených spinových trojúhelníků
• Simulovat časový vývoj pomocí Trotterizovaných obvodů různých hloubek
• Počítat a vizualizovat jednoqubitové magnetizace $\langle Z_i \rangle$ v čase
• Porovnávat základní simulace s výsledky z hardwarových spuštění pomocí řízení výkonu Q-CTRL Fire Opal

Přehled

Transverzální Isingův model (TFIM) je kvantový spinový model, který zachycuje základní rysy kvantových fázových přechodů. Hamiltonián je definován jako:

$$H = -J \sum_{i} Z_i Z_{i+1} - h \sum_{i} X_i$$

kde $Z_i$ a $X_i$ jsou Pauliho operátory působící na Qubit $i$, $J$ je síla vazby mezi sousedními spiny a $h$ je síla transverzálního magnetického pole. První člen reprezentuje klasické feromagnetické interakce, zatímco druhý zavádí kvantové fluktuace prostřednictvím transverzálního pole. K simulaci dynamiky TFIM se používá Trotterova dekompozice unitárního operátoru vývoje $e^{-iHt}$, implementovaná pomocí vrstev hradel RX a RZZ na základě vlastního grafu propojených spinových trojúhelníků. Simulace zkoumá, jak se magnetizace $\langle Z \rangle$ vyvíjí s rostoucím počtem Trotterových kroků.

Výkon navrhované implementace TFIM je hodnocen porovnáním bezšumových simulací s šumovými backendy. Vylepšené funkce provádění a potlačování chyb od Fire Opal se používají ke zmírnění dopadu šumu na reálném hardwaru, čímž se získávají spolehlivější odhady spinových observabilních veličin, jako jsou $\langle Z_i \rangle$ a korelátoři $\langle Z_i Z_j \rangle$.

Požadavky

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

• Qiskit SDK v1.4 nebo novější, s podporou vizualizace
• Qiskit Runtime v0.40 nebo novější (pip install qiskit-ibm-runtime)
• Qiskit Functions Catalog v0.9.0 (pip install qiskit-ibm-catalog)
• Fire Opal SDK v9.0.2 nebo novější (pip install fire-opal)
• Q-CTRL Visualizer v8.0.2 nebo novější (pip install qctrl-visualizer)

Nastavení

Nejprve se ověř pomocí svého IBM Quantum API klíče. Pak vyber Qiskit Function takto. (Tento kód předpokládá, že jsi již uložil/a svůj účet do místního prostředí.)

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

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit import QuantumCircuit
from qiskit_ibm_catalog import QiskitFunctionsCatalog
from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import SamplerV2 as Sampler
from qiskit.quantum_info import SparsePauliOp
from qiskit_aer import AerSimulator

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import qctrlvisualizer as qv

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Access Function
perf_mgmt = catalog.load("q-ctrl/performance-management")

Krok 1: Mapování klasických vstupů na kvantový problém

Generování TFIM grafu

Začínáme definováním mřížky spinů a vazeb mezi nimi. V tomto tutoriálu je mřížka sestavena z propojených trojúhelníků uspořádaných do lineárního řetězce. Každý trojúhelník se skládá ze tří uzlů propojených v uzavřené smyčce a řetězec je tvořen propojením jednoho uzlu každého trojúhelníku s předchozím trojúhelníkem.

Pomocná funkce connected_triangles_adj_matrix sestaví matici sousednosti pro tuto strukturu. Pro řetězec $n$ trojúhelníků má výsledný graf $2n+1$ uzlů.

def connected_triangles_adj_matrix(n):
    """
    Generate the adjacency matrix for 'n' connected triangles in a chain.
    """
    num_nodes = 2 * n + 1
    adj_matrix = np.zeros((num_nodes, num_nodes), dtype=int)

    for i in range(n):
        a, b, c = i * 2, i * 2 + 1, i * 2 + 2  # Nodes of the current triangle

        # Connect the three nodes in a triangle
        adj_matrix[a, b] = adj_matrix[b, a] = 1
        adj_matrix[b, c] = adj_matrix[c, b] = 1
        adj_matrix[a, c] = adj_matrix[c, a] = 1

        # If not the first triangle, connect to the previous triangle
        if i > 0:
            adj_matrix[a, a - 1] = adj_matrix[a - 1, a] = 1

    return adj_matrix

Pro vizualizaci mřížky, kterou jsme právě definovali, můžeme vykreslit řetězec propojených trojúhelníků a označit každý uzel. Níže uvedená funkce sestaví graf pro zvolený počet trojúhelníků a zobrazí jej.

def plot_triangle_chain(n, side=1.0):
    """
    Plot a horizontal chain of n equilateral triangles.
    Baseline: even nodes (0,2,4,...,2n) on y=0
    Apexes: odd nodes (1,3,5,...,2n-1) above the midpoint.
    """
    # Build graph
    A = connected_triangles_adj_matrix(n)
    G = nx.from_numpy_array(A)

    h = np.sqrt(3) / 2 * side
    pos = {}

    # Place baseline nodes
    for k in range(n + 1):
        pos[2 * k] = (k * side, 0.0)

    # Place apex nodes
    for k in range(n):
        x_left = pos[2 * k][0]
        x_right = pos[2 * k + 2][0]
        pos[2 * k + 1] = ((x_left + x_right) / 2, h)

    # Draw
    fig, ax = plt.subplots(figsize=(1.5 * n, 2.5))
    nx.draw(
        G,
        pos,
        ax=ax,
        with_labels=True,
        font_size=10,
        font_color="white",
        node_size=600,
        node_color=qv.QCTRL_STYLE_COLORS[0],
        edge_color="black",
        width=2,
    )
    ax.set_aspect("equal")
    ax.margins(0.2)
    plt.show()

    return G, pos

V tomto tutoriálu použijeme řetězec 20 trojúhelníků.

n_triangles = 20
n_qubits = 2 * n_triangles + 1
plot_triangle_chain(n_triangles, side=1.0)
plt.show()

Barvení hran grafu

Pro implementaci vazby spin–spin je užitečné seskupit hrany, které se nepřekrývají. To nám umožňuje aplikovat dvoququbitová hradla paralelně. Lze to provést pomocí jednoduché procedury barvení hran [1], která přiřadí každé hraně barvu tak, aby hrany sdílející stejný uzel byly v různých skupinách.

def edge_coloring(graph):
    """
    Takes a NetworkX graph and returns a list of lists where each inner list contains
    the edges assigned the same color.
    """
    line_graph = nx.line_graph(graph)
    edge_colors = nx.coloring.greedy_color(line_graph)

    color_groups = {}
    for edge, color in edge_colors.items():
        if color not in color_groups:
            color_groups[color] = []
        color_groups[color].append(edge)

    return list(color_groups.values())

Krok 2: Optimalizace problému pro provádění na kvantovém hardwaru

Generování Trotterizovaných obvodů na spinových grafech

Pro simulaci dynamiky TFIM konstruujeme obvody, které aproximují operátor časového vývoje.

$$U(t) = e^{-i H t}, \quad \text{where} \quad H = -J \sum_{\langle i,j \rangle} Z_i Z_j - h \sum_i X_i .$$

Používáme Trotterovu dekompozici druhého řádu:

$$e^{-i H \Delta t} \approx e^{-i H_X \Delta t / 2}\, e^{-i H_Z \Delta t}\, e^{-i H_X \Delta t / 2},$$

kde $H_X = -h \sum_i X_i$ a $H_Z = -J \sum_{\langle i,j \rangle} Z_i Z_j$.

• Člen $H_X$ je implementován pomocí vrstev rotací RX.
• Člen $H_Z$ je implementován pomocí vrstev hradel RZZ podél hran interakčního grafu.

Úhly těchto hradel jsou určeny transverzálním polem $h$, konstantou vazby $J$ a časovým krokem $\Delta t$. Skládáním více Trotterových kroků generujeme obvody rostoucí hloubky, které aproximují dynamiku systému. Funkce generate_tfim_circ_custom_graph a trotter_circuits konstruují Trotterizovaný kvantový Circuit z libovolného grafu spinových interakcí.

def generate_tfim_circ_custom_graph(
    steps, h, J, dt, psi0, graph: nx.graph.Graph, meas_basis="Z", mirror=False
):
    """
    Generate a second order trotter of the form e^(a+b) ~ e^(b/2) e^a e^(b/2) for simulating a transverse field ising model:
    e^{-i H t} where the Hamiltonian H = -J \\sum_i Z_i Z_{i+1} + h \\sum_i X_i.

    steps: Number of trotter steps
    theta_x: Angle for layer of X rotations
    theta_zz: Angle for layer of ZZ rotations
    theta_x: Angle for second layer of X rotations
    J: Coupling between nearest neighbor spins
    h: The transverse magnetic field strength
    dt: t/total_steps
    psi0: initial state (assumed to be prepared in the computational basis).
    meas_basis: basis to measure all correlators in

    This is a second order trotter of the form e^(a+b) ~ e^(b/2) e^a e^(b/2)
    """
    theta_x = h * dt
    theta_zz = -2 * J * dt
    nq = graph.number_of_nodes()
    color_edges = edge_coloring(graph)
    circ = QuantumCircuit(nq, nq)
    # Initial state, for typical cases in the computational basis
    for i, b in enumerate(psi0):
        if b == "1":
            circ.x(i)
    # Trotter steps
    for step in range(steps):
        for i in range(nq):
            circ.rx(theta_x, i)
        if mirror:
            color_edges = [sublist[::-1] for sublist in color_edges[::-1]]
        for edge_list in color_edges:
            for edge in edge_list:
                circ.rzz(theta_zz, edge[0], edge[1])
        for i in range(nq):
            circ.rx(theta_x, i)

    # some typically used basis rotations
    if meas_basis == "X":
        for b in range(nq):
            circ.h(b)
    elif meas_basis == "Y":
        for b in range(nq):
            circ.sdg(b)
            circ.h(b)

    for i in range(nq):
        circ.measure(i, i)

    return circ

def trotter_circuits(G, d_ind_tot, J, h, dt, meas_basis, mirror=True):
    """
    Generates a sequence of Trotterized circuits, each with increasing depth.
    Given a spin interaction graph and Hamiltonian parameters, it constructs
    a list of circuits with 1 to d_ind_tot Trotter steps

    G: Graph defining spin interactions (edges = ZZ couplings)
    d_ind_tot: Number of Trotter steps (maximum depth)
    J: Coupling between nearest neighboring spins
    h: Transverse magnetic field strength
    dt: (t / total_steps
    meas_basis: Basis to measure all correlators in
    mirror: If True, mirror the Trotter layers
    """
    qubit_count = len(G)
    circuits = []
    psi0 = "0" * qubit_count

    for steps in range(1, d_ind_tot + 1):
        circuits.append(
            generate_tfim_circ_custom_graph(
                steps, h, J, dt, psi0, G, meas_basis, mirror
            )
        )
    return circuits

Odhadni magnetizace jednotlivých Qubitů $\langle Z_i \rangle$

Abychom mohli studovat dynamiku modelu, chceme měřit magnetizaci každého Qubitu, definovanou střední hodnotou $\langle Z_i \rangle = \langle \psi | Z_i | \psi \rangle$.

V simulacích to můžeme vypočítat přímo z výsledků měření. Funkce z_expectation zpracovává počty bitových řetězců a vrací hodnotu $\langle Z_i \rangle$ pro zvolený index Qubitu. Na skutečném hardware vyhodnocujeme tutéž veličinu zadáním Pauliho operátoru pomocí funkce generate_z_observables, a pak Backend vypočítá střední hodnotu.

def z_expectation(counts, index):
    """
    counts: Dict of mitigated bitstrings.
    index: Index i in the single operator expectation value < II...Z_i...I > to be calculated.
    return:  < Z_i >
    """
    z_exp = 0
    tot = 0
    for bitstring, value in counts.items():
        bit = int(bitstring[index])
        sign = 1
        if bit % 2 == 1:
            sign = -1
        z_exp += sign * value
        tot += value

    return z_exp / tot

def generate_z_observables(nq):
    observables = []
    for i in range(nq):
        pauli_string = "".join(["Z" if j == i else "I" for j in range(nq)])
        observables.append(SparsePauliOp(pauli_string))
    return observables

observables = generate_z_observables(n_qubits)

Nyní definujeme parametry pro generování Trotterizovaných Circuitů. V tomto tutoriálu je mřížka řetězec 20 propojených trojúhelníků, což odpovídá systému s 41 Qubity.

all_circs_mirror = []
for num_triangles in [n_triangles]:
    for meas_basis in ["Z"]:
        A = connected_triangles_adj_matrix(num_triangles)
        G = nx.from_numpy_array(A)
        nq = len(G)
        d_ind_tot = 22
        dt = 2 * np.pi * 1 / 30 * 0.25
        J = 1
        h = -7
        all_circs_mirror.extend(
            trotter_circuits(G, d_ind_tot, J, h, dt, meas_basis, True)
        )
circs = all_circs_mirror

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

Spuštění MPS simulace

Seznam Trotterizovaných Circuit je spuštěn pomocí simulátoru matrix_product_state s libovolně zvoleným počtem $4096$ výstřelů. Metoda MPS poskytuje efektivní aproximaci dynamiky Circuit, přičemž přesnost je určena zvolenou rozměrností vazby. Pro velikosti systémů uvažované v tomto tutoriálu je výchozí rozměrnost vazby dostatečná k zachycení dynamiky magnetizace s vysokou věrností. Surové počty jsou normalizovány a z nich vypočítáme střední hodnoty jednoQubitového operátoru $\langle Z_i \rangle$ v každém Trotterově kroku. Nakonec vypočítáme průměr přes všechny Qubity a získáme jedinou křivku ukazující, jak se magnetizace mění v čase.

backend_sim = AerSimulator(method="matrix_product_state")

def normalize_counts(counts_list, shots):
    new_counts_list = []
    for counts in counts_list:
        a = {k: v / shots for k, v in counts.items()}
        new_counts_list.append(a)
    return new_counts_list

def run_sim(circ_list):
    shots = 4096
    res = backend_sim.run(circ_list, shots=shots)
    normed = normalize_counts(res.result().get_counts(), shots)
    return normed

sim_counts = run_sim(circs)

Spuštění na hardwaru

service = QiskitRuntimeService()
backend = service.backend("ibm_marrakesh")

def run_qiskit(circ_list):
    shots = 4096
    pm = generate_preset_pass_manager(backend=backend)
    isa_circuits = [pm.run(qc) for qc in circ_list]
    sampler = Sampler(mode=backend)
    res = sampler.run(isa_circuits, shots=shots)
    res = [r.data.c.get_counts() for r in res.result()]
    normed = normalize_counts(res, shots)
    return normed

qiskit_counts = run_qiskit(circs)

Spuštění na hardwaru s Fire Opal

Vyhodnotíme dynamiku magnetizace na skutečném kvantovém hardwaru. Fire Opal poskytuje Qiskit Function, která rozšiřuje standardní primitivum Qiskit Runtime Estimator o automatické potlačení chyb a správu výkonu. Trotterizované Circuit odesíláme přímo do Backend IBM®, zatímco Fire Opal se stará o spouštění s ohledem na šum.

Připravíme seznam pubs, kde každá položka obsahuje Circuit a odpovídající Pauli-Z observabely. Ty jsou předány funkci Estimator od Fire Opal, která vrátí střední hodnoty $\langle Z_i \rangle$ pro každý Qubit v každém Trotterově kroku. Výsledky lze poté zprůměrovat přes Qubity a získat křivku magnetizace z hardwaru.

backend_name = "ibm_marrakesh"
estimator_pubs = [(qc, observables) for qc in all_circs_mirror[:]]

# Run the circuit using the estimator
qctrl_estimator_job = perf_mgmt.run(
    primitive="estimator",
    pubs=estimator_pubs,
    backend_name=backend_name,
    options={"default_shots": 4096},
)

result_qctrl = qctrl_estimator_job.result()

Krok 4: Post-processing a vrácení výsledku v požadovaném klasickém formátu

Nakonec porovnáme křivku magnetizace ze simulátoru s výsledky získanými na skutečném hardwaru. Vykreslení obou vedle sebe ukazuje, jak přesně spuštění na hardwaru s Fire Opal odpovídá bezšumovému základnímu průběhu napříč Trotterovými kroky.

def make_correlators(test_counts, nq, d_ind_tot):
    mz = np.empty((nq, d_ind_tot))
    for d_ind in range(d_ind_tot):
        counts = test_counts[d_ind]
        for i in range(nq):
            mz[i, d_ind] = z_expectation(counts, i)
    average_z = np.mean(mz, axis=0)
    return np.concatenate((np.array([1]), average_z), axis=0)

sim_exp = make_correlators(sim_counts[0:22], nq=nq, d_ind_tot=22)
qiskit_exp = make_correlators(qiskit_counts[0:22], nq=nq, d_ind_tot=22)

qctrl_exp = [ev.data.evs for ev in result_qctrl[:]]
qctrl_exp_mean = np.concatenate(
    (np.array([1]), np.mean(qctrl_exp, axis=1)), axis=0
)

def make_expectations_plot(
    sim_z,
    depths,
    exp_qctrl=None,
    exp_qctrl_error=None,
    exp_qiskit=None,
    exp_qiskit_error=None,
    plot_from=0,
    plot_upto=23,
):
    import numpy as np
    import matplotlib.pyplot as plt

    depth_ticks = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22]

    d = np.asarray(depths)[plot_from:plot_upto]
    sim = np.asarray(sim_z)[plot_from:plot_upto]

    qk = (
        None
        if exp_qiskit is None
        else np.asarray(exp_qiskit)[plot_from:plot_upto]
    )
    qc = (
        None
        if exp_qctrl is None
        else np.asarray(exp_qctrl)[plot_from:plot_upto]
    )

    qk_err = (
        None
        if exp_qiskit_error is None
        else np.asarray(exp_qiskit_error)[plot_from:plot_upto]
    )
    qc_err = (
        None
        if exp_qctrl_error is None
        else np.asarray(exp_qctrl_error)[plot_from:plot_upto]
    )

    # ---- helper(s) ----
    def rmse(a, b):
        if a is None or b is None:
            return None
        a = np.asarray(a, dtype=float)
        b = np.asarray(b, dtype=float)
        mask = np.isfinite(a) & np.isfinite(b)
        if not np.any(mask):
            return None
        diff = a[mask] - b[mask]
        return float(np.sqrt(np.mean(diff**2)))

    def plot_panel(ax, method_y, method_err, color, label, band_color=None):
        # Noiseless reference
        ax.plot(d, sim, color="grey", label="Noiseless simulation")

        # Method line + band
        if method_y is not None:
            ax.plot(d, method_y, color=color, label=label)
            if method_err is not None:
                lo = np.clip(method_y - method_err, -1.05, 1.05)
                hi = np.clip(method_y + method_err, -1.05, 1.05)
                ax.fill_between(
                    d,
                    lo,
                    hi,
                    alpha=0.18,
                    color=band_color if band_color else color,
                    label=f"{label} ± error",
                )
        else:
            ax.text(
                0.5,
                0.5,
                "No data",
                transform=ax.transAxes,
                ha="center",
                va="center",
                fontsize=10,
                color="0.4",
            )

        # RMSE box (vs sim)
        r = rmse(method_y, sim)
        if r is not None:
            ax.text(
                0.98,
                0.02,
                f"RMSE: {r:.4f}",
                transform=ax.transAxes,
                va="bottom",
                ha="right",
                fontsize=8,
                bbox=dict(
                    boxstyle="round,pad=0.35", fc="white", ec="0.7", alpha=0.9
                ),
            )
        # Axes
        ax.set_xticks(depth_ticks)
        ax.set_ylim(-1.05, 1.05)
        ax.grid(True, which="both", linewidth=0.4, alpha=0.4)
        ax.set_axisbelow(True)
        ax.legend(prop={"size": 8}, loc="best")

    fig, axes = plt.subplots(1, 2, figsize=(10, 4), dpi=300, sharey=True)

    axes[0].set_title("Fire Opal (Q-CTRL)", fontsize=10)
    plot_panel(
        axes[0],
        qc,
        qc_err,
        color="#680CE9",
        label="Fire Opal",
        band_color="#680CE9",
    )
    axes[0].set_xlabel("Trotter step")
    axes[0].set_ylabel(r"$\langle Z \rangle$")
    axes[1].set_title("Qiskit", fontsize=10)
    plot_panel(
        axes[1], qk, qk_err, color="blue", label="Qiskit", band_color="blue"
    )
    axes[1].set_xlabel("Trotter step")

    plt.tight_layout()
    plt.show()

depths = list(range(d_ind_tot + 1))
errors = np.abs(np.array(qctrl_exp_mean) - np.array(sim_exp))

errors_qiskit = np.abs(np.array(qiskit_exp) - np.array(sim_exp))

make_expectations_plot(
    sim_exp,
    depths,
    exp_qctrl=qctrl_exp_mean,
    exp_qctrl_error=errors,
    exp_qiskit=qiskit_exp,
    exp_qiskit_error=errors_qiskit,
)

Reference

[1] Graph coloring. Wikipedia. Retrieved September 15, 2025, from https://en.wikipedia.org/wiki/Graph_coloring

Průzkum tutoriálu

Věnuj chvíli poskytnutí zpětné vazby k tomuto tutoriálu. Tvoje postřehy nám pomohou zlepšit nabídku obsahu a uživatelský zážitek.

Odkaz na průzkum

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.