Zasymuluj zaszumione układy kwantowe za pomocą propagacji Pauliego

W tym przewodniku używamy pakietu pauli-prop, aby klasycznie zasymulować dynamikę czasową zaszumionego dziewięciokubitowego modelu Isinga z polem poprzecznym (TFIM) na kwadratowej sieci 3x3. Używamy instrukcji PauliLindbladError, aby zdefiniować kanał szumu, $\Lambda$, działający na zbiór warstw splątujących, $\mathcal{U}$. Następnie propagujemy obserwablę, $O$, wstecz przez zaszumiony obwód i szacujemy wartości oczekiwane dla różnych modeli szumu, a także dla przypadku bezszumowego.

Gdy obserwabla jest propagowana wstecz przez obwód, każdy kanał szumu, $\Lambda_k$, powiązany z warstwą splątującą, $\mathcal{U}_k$, tłumi te wyrazy Pauliego w $O$, które antykomutują z jego generatorami Pauliego-Lindblada. Dokładniej, jeśli $G_{k,i}$ jest generatorem Pauliego $\Lambda_k$ o szybkości, $\gamma_{k,i}$, to wyraz Pauliego, $P$, w $O$ przekształca się jako: $c_P \mapsto c_P e^{-2\gamma_{k,i}} \quad \text{if } \{P, G_{k,i}\}=0$, gdzie $c_P$ jest współczynnikiem $P$. Gdy $O$ zostanie spropagowana na początek obwodu, wartość oczekiwaną względem stanu zerowego, $|0\rangle^{\otimes N}$, można trywialnie obliczyć, sumując współczynniki każdego wyrazu diagonalnego w $O$ (wyrazów zawierających $Z$ lub $I$ na wszystkich kubitach).

Przebieg pracy:

• Określ sieć TFIM i użyj kolorowania krawędzi, aby wyznaczyć minimalny zbiór warstw splątujących
• Wygeneruj syntetyczne modele szumu, $\Lambda_k$, dla każdej unikalnej warstwy splątującej, $U_k$
  • Utwórz modele szumu o różnych skalach, aby zbadać wpływ szumu bramek na układ
• Utwórz bezszumowe i zaszumione obwody kwantowe dla interesujących Cię głębokości i skal szumu
  • W zaszumionych obwodach instrukcje PauliLindbladError są wstawiane przed każdą warstwą splątującą
• Użyj propagacji Pauliego, aby zasymulować dokładne wartości oczekiwane układu przy różnych głębokościach
  • Dla dziewięciu kubitów robi się to, pozwalając $O$ urosnąć do $4^9$ wyrazów, co obejmuje całą przestrzeń Pauliego
• Użyj propagacji Pauliego, aby zasymulować zaszumione wartości oczekiwane
• Zaobserwuj, jak rosnący szum bramek pogarsza dokładność modelu kwantowego

Wygeneruj kwadratową sieć 3x3 i znajdź 4-kolorowanie krawędzi

Wierzchołki grafu reprezentują kubity, a krawędzie reprezentują połączenie między dwoma kubitami. Kolorowanie krawędzi odpowiada unikalnym warstwom splątującym w obwodzie kwantowym w taki sposób, że bramek na połączeniach o różnych kolorach nie można stosować jednocześnie.

Wyznaczenie minimalnego zbioru unikalnych warstw splątujących jest często ważne dla wydajnych protokołów uczenia szumu, ponieważ szum każdej warstwy trzeba poznać niezależnie. Im więcej warstw musimy poznać, tym więcej strzałów musimy pobrać z QPU. W tym demo używamy informacji o warstwach, aby zbudować zaszumione obwody i wstawić instrukcje PauliLindbladError z qiskit-aer przed każdą warstwą splątującą, modelując szum bramek 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.

Wygeneruj syntetyczne modele szumu

Przed utworzeniem obwodów kwantowych generujemy model szumu (instancję PauliLindbladError) dla każdej z warstw splątujących. Później osadzimy je jako instrukcje w naszych obwodach kwantowych. Dla każdej warstwy generujemy kanały szumu o różnych skalach. Konkretnie generujemy modele szumu o błędzie na warstwową bramkę (Error Per Layered Gate, EPLG) wynoszącym w przybliżeniu .0004, .0008, .0012, .0016, i .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
]

Utwórz obwody kwantowe

W tym demo symulujemy dynamikę czasową modelu Isinga z polem poprzecznym (TFIM) dla rosnącej liczby kroków Trottera (1-10 kroków). Dla każdej z 10 głębokości obwodu symulujemy wpływ szumu bramek, dla modeli szumu o różnych skalach (EPLGs = .0004, .0008, .0012, .0016, .002). Szum jest wstawiany do QuantumCircuit jako instrukcja PauliLindbladError z Qiskit Aer. Rozważany hamiltonian ma postać:

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

gdzie $J>0$ opisuje sprzężenie spinów najbliższych sąsiadów, $i<j$, a $h$ jest globalnym polem poprzecznym.

Tutaj implementujemy hamiltonian ewoluujący w czasie dla różnych skal czasu i szumu. Tworzymy łącznie 60 obwodów: 10 bezszumowych obwodów różniących się głębokością Trottera oraz 50 zaszumionych obwodów dla 10 głębokości Trottera w pięciu skalach szumu. Dla danego grafu połączeń model jest sparametryzowany kilkoma zmiennymi:

• num_steps: liczba kroków Trottera
• J: siła sprzężenia połączonych węzłów
• h: siła zewnętrznego pola magnetycznego
• dt: zmiana czasu w jednym kroku Trottera
• initial_state_angle: początkowe wzbudzenie, $R_y(\theta)$, stosowane jednorodnie do wszystkich kubitów

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.

Określ obserwablę i uruchom symulacje

W tym demo symulujemy wartości oczekiwane średniego korelatora dwuwęzłowego:

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

gdzie $\Psi(\theta)$ odpowiada jednorodnemu obrotowi $R_y(\theta)$ na wszystkich kubitach, $\mathscr{U}^s$ opisuje $s$ warstw Trottera, a $(j,k)$ indeksuje wszystkie połączone pary wierzchołków sieci.

Na koniec używamy pauli_prop, aby zasymulować wartości oczekiwane obserwabli dla każdego z obwodów. W tym dziewięciokubitowym demo wszystkie symulacje wykonujemy dokładnie. Nie będzie stosowane żadne obcinanie w propagacji Pauliego; dlatego różnice wartości oczekiwanych między różnymi modelami szumu można w całości przypisać błędowi bramek. Proces symulacji składa się z czterech kroków:

• Przesuń bramki Clifforda w obwodzie na jego początek za pomocą pauli_prop.evolve_through_cliffords
• Propaguj obserwablę przez niecliffordowską część obwodu za pomocą pauli_prop.propagate_through_circuit
  • Wykonujemy dokładne symulacje, pozwalając obserwablej urosnąć do rozmiaru całej przestrzeni Pauliego, $4^9$
• Propaguj ewoluowaną obserwablę przez cliffordowską część obwodu za pomocą SparsePauliOp.evolve z Qiskit
• Oszacuj wartość oczekiwaną względem stanu zerowego, $|0\rangle^{\otimes N}$, sumując współczynniki każdego wyrazu diagonalnego w $O$ (wyrazów zawierających $Z$ lub $I$ na wszystkich kubitach)

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.

Zaobserwuj wpływ błędu bramek na model

Pamiętaj, że ponieważ jest to eksperyment dziewięciokubitowy, procedura propagacji Pauliego jest dokładna, a cały błąd na wykresach szumu można przypisać błędowi bramek.

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