Simulace zašuměných kvantových systémů pomocí Pauliho propagace

V tomto průvodci používáme balíček pauli-prop ke klasické simulaci časové dynamiky zašuměného devítiqubitového Isingova modelu s příčným polem (TFIM) na čtvercové mřížce 3x3. Pomocí instrukcí PauliLindbladError definujeme šumový kanál $\Lambda$ působící na množinu provázávacích vrstev $\mathcal{U}$. Poté propagujeme observabli $O$ zpětně zašuměným obvodem a odhadujeme střední hodnoty pro různé modely šumu i pro bezšumový případ.

Při zpětné propagaci observabli obvodem každý šumový kanál $\Lambda_k$ příslušející provázávací vrstvě $\mathcal{U}_k$ tlumí Pauliho členy v $O$, které antikomutují s jeho Pauli-Lindbladovými generátory. Přesněji, pokud je $G_{k,i}$ Pauliho generátor kanálu $\Lambda_k$ s mírou $\gamma_{k,i}$, pak se Pauliho člen $P$ v $O$ transformuje takto: $c_P \mapsto c_P e^{-2\gamma_{k,i}} \quad \text{if } \{P, G_{k,i}\}=0$, kde $c_P$ je koeficient členu $P$. Jakmile je $O$ propagována na začátek obvodu, lze střední hodnotu vůči nulovému stavu, $|0\rangle^{\otimes N}$, triviálně vypočítat sečtením koeficientů každého diagonálního členu v $O$ (členy obsahující na všech qubitech $Z$ nebo $I$).

Postup:

• Zadej mřížku TFIM a pomocí barvení hran urči minimální množinu provázávacích vrstev
• Vygeneruj syntetické modely šumu, $\Lambda_k$, pro každou jedinečnou provázávací vrstvu, $U_k$
  • Vytvoř modely šumu různých úrovní, abys prozkoumal vliv šumu hradel na systém
• Vytvoř bezšumové a zašuměné kvantové obvody pro různé hloubky a úrovně šumu, které tě zajímají
  • V zašuměných obvodech jsou instrukce PauliLindbladError vkládány před každou provázávací vrstvu
• Pomocí Pauliho propagace simuluj přesné střední hodnoty systému v různých hloubkách
  • Pro devět qubitů se to provádí tak, že necháme $O$ narůst na $4^9$ členů, což pokrývá celý Pauliho prostor
• Pomocí Pauliho propagace simuluj zašuměné střední hodnoty
• Pozoruj, jak rostoucí šum hradel zhoršuje přesnost kvantového modelu

Vygenerování čtvercové mřížky 3x3 a nalezení 4-obarvení hran

Vrcholy grafu představují qubity a hrany představují spojení mezi dvěma qubity. Obarvení hran odpovídá jedinečným provázávacím vrstvám v kvantovém obvodu tak, že hradla na spojeních s různými barvami nelze aplikovat současně.

Určení minimální množiny jedinečných provázávacích vrstev je často důležité pro efektivní protokoly učení šumu, protože šum každé vrstvy se musí učit nezávisle. Čím více vrstev se musíme učit, tím více shotů musíme z QPU odebrat. V tomto dema používáme informace o vrstvách k sestavení zašuměných obvodů a před každou provázávací vrstvu vkládáme instrukce PauliLindbladError z qiskit-aer, abychom modelovali šum hradel QPU.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "pauli_prop": "pauli-prop", "qiskit": "qiskit", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "rustworkx": "rustworkx"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

from collections import defaultdict

import numpy as np
from qiskit.transpiler import CouplingMap
from qiskit_addon_utils.coloring import auto_color_edges

# Define rectangular square-lattice on 20 qubits
num_rows = 3
num_cols = 3
num_qubits = num_rows * num_cols

coupling_map = CouplingMap.from_grid(
    num_rows=num_rows, num_columns=num_cols, bidirectional=False
)

# Create mapping from color to edge list
coloring = auto_color_edges(coupling_map.get_edges())
color_to_edge = defaultdict(list)
for edge, color in coloring.items():
    color_to_edge[color].append(edge)

from rustworkx import PyDiGraph
from rustworkx.visualization import graphviz_draw

# Inspect graph coupling and unique entangling layers
print(
    f"The circuit will have {num_qubits} qubits and {len(color_to_edge)} unique entangling layers."
)
sq_lattice = PyDiGraph()
sq_lattice.extend_from_weighted_edge_list(
    [
        (source, target, color)
        for ((source, target), color) in coloring.items()
    ]
)

def color_edge_4color(edge):
    color_dict = {0: "red", 1: "green", 2: "blue", 3: "orange"}
    return {"color": color_dict[edge]}

graphviz_draw(sq_lattice, edge_attr_fn=color_edge_4color, method="neato")

The circuit will have 9 qubits and 4 unique entangling layers.

Vygenerování syntetických modelů šumu

Před vytvořením kvantových obvodů vygenerujeme pro každou z provázávacích vrstev model šumu (instanci PauliLindbladError). Později je vložíme jako instrukce do našich kvantových obvodů. Pro každou vrstvu generujeme šumové kanály různých úrovní. Konkrétně generujeme modely šumu s Error Per Layered Gate (EPLG) přibližně .0004, .0008, .0012, .0016, a .002.

from qiskit.quantum_info import SparsePauliOp, pauli_basis
from qiskit_aer.noise import PauliLindbladError

# Pauli-Lindblad noise parameters
seed = 1764
target_EPLGs = [0.0004, 0.0008, 0.0012, 0.0016, 0.002]

def generate_random_pauli_lindblad_noise(
    edges,
    num_qubits: int | None = None,
    noise_scale: float = 1e-3,
    seed: int | None = None,
) -> PauliLindbladError:
    """Generate random Pauli-Lindblad noise over the full Pauli basis."""
    if num_qubits is None:
        num_qubits = np.max(edges)

    basis_paulis = [p for p in pauli_basis(2) if np.sum(p.x + p.z)]
    basis_paulis = SparsePauliOp.from_sparse_list(
        [
            (pauli.to_label(), edge, 1)
            for pauli in basis_paulis
            for edge in edges
        ],
        num_qubits=num_qubits,
    )
    basis_paulis = basis_paulis.simplify()
    basis_paulis = basis_paulis.paulis

    rng = np.random.default_rng(seed=seed)
    rates = rng.random(len(basis_paulis)) * noise_scale

    return PauliLindbladError(generators=basis_paulis, rates=rates)

num_generators = (
    (num_rows * num_cols)
    + (num_rows - 1) * num_cols
    + num_rows * (num_cols - 1)
)
noise_scales = [
    EPLG * (num_rows * num_cols) / num_generators for EPLG in target_EPLGs
]
noise_models_per_EPLG = [
    [
        generate_random_pauli_lindblad_noise(
            color_to_edge[color],
            num_qubits=num_qubits,
            noise_scale=noise_scale,
            seed=seed,
        )
        for color in range(len(color_to_edge))
    ]
    for noise_scale in noise_scales
]

Vytvoření kvantových obvodů

V tomto dema simulujeme časovou dynamiku Isingova modelu s příčným polem (TFIM) pro rostoucí počet Trotterových kroků (1–10 kroků). Pro každou z 10 hloubek obvodu simulujeme vliv šumu hradel při modelech šumu různých úrovní (EPLGs = .0004, .0008, .0012, .0016, .002). Šum je do QuantumCircuit vložen jako instrukce PauliLindbladError z Qiskit Aer. Uvažovaný hamiltonián je:

$H = -J\sum\limits_{\langle i,j \rangle} Z_iZ_j + h\sum\limits_iX_i$

kde $J>0$ popisuje vazbu mezi nejbližšími sousedními spiny, $i<j$, a $h$ je globální příčné pole.

Zde implementujeme časově vyvíjený hamiltonián napříč různými časovými a šumovými škálami. Vytvoříme celkem 60 obvodů: 10 bezšumových obvodů lišících se Trotterovou hloubkou a 50 zašuměných obvodů pro 10 Trotterových hloubek napříč pěti úrovněmi šumu. Pro daný graf konektivity je model parametrizován několika proměnnými:

• num_steps: Počet Trotterových kroků
• J: Síla vazby propojených míst
• h: Síla vnějšího magnetického pole
• dt: Změna času během jednoho Trotterova kroku
• initial_state_angle: Počáteční excitace, $R_y(\theta)$, aplikovaná rovnoměrně na všechny qubity

from typing import Any

from qiskit import QuantumCircuit

# Ising model parameters
num_steps = 10
J = -1.0
dt = 0.25 / abs(J)
h = 2.0 * abs(J)
initial_state_angle = np.pi / 18.0
rx_angle = 2.0 * h * dt
rzz_angle = 2.0 * J * dt

def generate_ising_circuit(
    num_qubits: int,
    num_steps: int,
    rx_angle: float,
    rzz_angle: float,
    coloring: dict[Any, list[tuple[int, int]]],
    layer_noise_models: list[PauliLindbladError] | None = None,
    initial_state_angle: float | None = None,
) -> QuantumCircuit:
    """Generate a quantum circuit implementing a transverse-field Ising model"""
    qc = QuantumCircuit(num_qubits)
    if initial_state_angle:
        qc.ry(initial_state_angle, range(num_qubits))
    qc.rx(rx_angle / 2, range(num_qubits))
    for i in range(num_steps):
        for j, layer in enumerate(coloring):
            edges = coloring[layer]
            if layer_noise_models:
                qc.append(layer_noise_models[j], qargs=range(num_qubits))
            for edge in edges:
                qc.rzz(rzz_angle, *edge)
        if i == num_steps - 1:
            qc.rx(rx_angle / 2, range(num_qubits))
        else:
            qc.rx(rx_angle, range(num_qubits))
    return qc

# Create the noiseless and noisy circuits
noiseless_circs = []
noisy_circs = []
for steps in range(1, num_steps + 1):
    noiseless_circs.append(
        generate_ising_circuit(
            num_qubits,
            steps,
            rx_angle,
            rzz_angle,
            color_to_edge,
            initial_state_angle=initial_state_angle,
        )
    )
    noisy_circs_per_step = []
    for noise_models in noise_models_per_EPLG:
        noisy_circs_per_step.append(
            generate_ising_circuit(
                num_qubits,
                steps,
                rx_angle,
                rzz_angle,
                color_to_edge,
                layer_noise_models=noise_models,
                initial_state_angle=initial_state_angle,
            )
        )
    noisy_circs.append(noisy_circs_per_step)
print(
    f"{num_steps} noiseless and {num_steps * len(target_EPLGs)} noisy Trotter circuits generated. {num_steps} different depths across {len(target_EPLGs)} different noise models"
)
print("\nBelow: Initial state and one noisy Trotter step.")
noisy_circs[0][0].draw("mpl", fold=-1)

10 noiseless and 50 noisy Trotter circuits generated. 10 different depths across 5 different noise models

Below: Initial state and one noisy Trotter step.

Zadání observabli a spuštění simulací

V tomto dema simulujeme střední hodnoty průměrného dvoumístného korelátoru:

$\langle O \rangle = \langle Z_{tot}^2(s) \rangle = \frac{1}{N^2}\sum \langle \Psi(\theta)|(\mathscr{U}^{\dagger})^sZ_jZ_k(\mathscr{U})^s|\Psi(\theta) \rangle$

kde $\Psi(\theta)$ odpovídá rovnoměrné rotaci $R_y(\theta)$ na všech qubitech, $\mathscr{U}^s$ popisuje $s$ Trotterových vrstev a $(j,k)$ indexuje všechny propojené dvojice vrcholů na mřížce.

Nakonec používáme pauli_prop k simulaci středních hodnot observabli pro každý z obvodů. V tomto devítiqubitovém dema provádíme všechny simulace přesně. Nebude provedeno žádné zkracování Pauliho propagace; rozdíly ve středních hodnotách mezi různými modely šumu lze proto zcela přičíst chybě hradel. Proces simulace probíhá ve čtyřech krocích:

• Přesuň cliffordovská hradla v obvodu na jeho začátek pomocí pauli_prop.evolve_through_cliffords
• Propaguj observabli necliffordovskou částí obvodu pomocí pauli_prop.propagate_through_circuit
  • Přesné simulace provádíme tak, že necháme observabli narůst na velikost celého Pauliho prostoru, $4^9$
• Propaguj vyvinutou observabli cliffordovskou částí obvodu pomocí SparsePauliOp.evolve z Qiskit
• Odhadni střední hodnotu vůči nulovému stavu, $|0\rangle^{\otimes N}$, sečtením koeficientů každého diagonálního členu v $O$ (členy obsahující na všech qubitech $Z$ nebo $I$)

import time

from pauli_prop import evolve_through_cliffords, propagate_through_circuit
from qiskit.quantum_info import Pauli

# Average ZZ-correlator observable
id_pauli = Pauli("I" * num_qubits)
observable = 2 * SparsePauliOp(
    [
        id_pauli.dot(Pauli("ZZ"), [i, j])
        for i in range(num_qubits)
        for j in range(i + 1, num_qubits)
    ]
)
observable /= num_qubits**2

# Pauli propagation parameters
max_terms = 4**num_qubits  # Exact propagation
atol = 1e-12

# Run simulations
exact_evs = []
noisy_evs = [[] for _ in range(len(target_EPLGs))]
st = time.perf_counter()
for i, noiseless_circ in enumerate(noiseless_circs):
    cliff, non_cliff = evolve_through_cliffords(noiseless_circ)
    evolved_obs = propagate_through_circuit(
        observable, non_cliff, max_terms=max_terms, atol=atol, frame="h"
    )[0]
    evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
    exact_evs.append(
        float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())
    )
    for j in range(len(target_EPLGs)):
        noisy_circ = noisy_circs[i][j]
        cliff, non_cliff = evolve_through_cliffords(noisy_circ)
        evolved_obs = propagate_through_circuit(
            observable, non_cliff, max_terms=max_terms, atol=1e-12, frame="h"
        )[0]
        evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame="h")
        noisy_evs[j].append(
            float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())
        )
print(
    f"Ran {len(noiseless_circs)} noiseless and {len(target_EPLGs) * num_steps} noisy simulations in {int(time.perf_counter() - st)}s."
)

Ran 10 noiseless and 50 noisy simulations in 103s.

Pozorování vlivu chyby hradel na model

Pamatuj, že protože jde o devítiqubitový experiment, rutina Pauliho propagace je přesná a veškerou chybu v grafech šumu lze přičíst chybě hradel.

import matplotlib.pyplot as plt

xs = range(1, num_steps + 1)
plt.plot(xs, exact_evs, label="Noiseless", color="black", marker="o")
colors = [".3", ".4", ".5", ".6", ".7"]
for i, evs in enumerate(noisy_evs):
    plt.plot(
        xs,
        evs,
        label=f"{target_EPLGs[i]} EPLG",
        linestyle="--",
        color=colors[i],
        marker="o",
    )
plt.xlabel("# Trotter steps")
plt.ylabel(r"$\langle Z_{tot}^2 \rangle$")
plt.legend()
plt.show()