{
"cells": [
{
"cell_type": "markdown",
"id": "frontmatter",
"metadata": {},
"source": [
"---\n",
"title: \"Simulate noisy quantum systems with Pauli propagation\"\n",
"description: \"Simulate noisy quantum systems with Pauli propagation for the latest version of Pauli propagation\"\n",
"---\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "91eb1163-e828-4c59-a06d-e4f0030f8071",
"metadata": {},
"source": [
"# Simulate noisy quantum systems with Pauli propagation\n",
"\n",
"{/* cspell:ignore mathscr, mapsto */}\n",
"\n",
"In this guide we use the `pauli-prop` package to classically simulate the time dynamics of a noisy nine-qubit transverse-field Ising model (TFIM) on a 3x3 square lattice. We use [PauliLindbladError](https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.PauliLindbladError.html) instructions to define a noise channel, $\\Lambda$, acting on a set of entangling layers, $\\mathcal{U}$. We then propagate the observable, $O$, backward through the noisy circuit and estimate expectation values for a variety of noise models, as well as the noiseless case.\n",
"\n",
"\n",
"\n",
"As the observable is propagated backwards through the circuit, each noise channel, $\\Lambda_k$, associated with entangling layer, $\\mathcal{U}_k$, damps the Pauli terms in $O$ which anti-commute with its Pauli-Lindblad generators. Specifically, if $G_{k,i}$ is a Pauli generator of $\\Lambda_k$ with rate, $\\gamma_{k,i}$, then a Pauli term, $P$, in $O$ transforms as: $c_P \\mapsto c_P e^{-2\\gamma_{k,i}} \\quad \\text{if } \\{P, G_{k,i}\\}=0$, where $c_P$ is the coefficient of $P$. Once $O$ has been propagated to the beginning of the circuit, the expectation value with respect to the zero state, $|0\\rangle^{\\otimes N}$, can be trivially calculated by summing the coefficients of each diagonal term in $O$ (terms containing $Z$ or $I$ on all qubits).\n",
"\n",
"Workflow:\n",
"\n",
"* Specify the TFIM lattice, and use edge coloring to identify a minimal set of entangling layers\n",
"* Generate synthetic noise models, $\\Lambda_k$, for each unique entangling layer, $U_k$\n",
" * Create noise models of various scales to study the impact of gate noise on the system\n",
"* Create noiseless and noisy quantum circuits for the various depths and noise scales of interest\n",
" * In noisy circuits, `PauliLindbladError` instructions are inserted before each entangling layer\n",
"* Use Pauli propagation to simulate exact expectation values of the system at various depths\n",
" * For nine qubits, this is done by letting $O$ grow to $4^9$ terms, covering the full Pauli space\n",
"* Use Pauli propagation to simulate noisy expectation values\n",
"* Observe how increasing gate noise degrades the accuracy of the quantum model\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b70d9f4d-2771-45fb-98ac-3c2c4ce978b5",
"metadata": {},
"source": [
"## Generate a 3x3 square lattice and find a 4-coloring on the edges\n",
"\n",
"The vertices in the graph represent qubits, and the edges represent a connection between two qubits. The edge coloring corresponds to unique entangling layers in the quantum circuit such that gates on connections associated with differing colors cannot be applied simultaneously.\n",
"\n",
"Identifying a minimal set of unique entangling layers is often important for implementing efficient noise-learning protocols, as the noise for each layer must be learned independently. The more layers we must learn, the more shots we need to take from the QPU. For this demo, we use the layer information to build up noisy circuits and inject [PauliLindbladError](https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.PauliLindbladError.html) instructions from `qiskit-aer` before each entangling layer to model QPU gate noise.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "58bcd725-3447-4bfc-8d63-12fd58082ff7",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"from collections import defaultdict\n",
"\n",
"import numpy as np\n",
"from qiskit.transpiler import CouplingMap\n",
"from qiskit_addon_utils.coloring import auto_color_edges\n",
"\n",
"# Define rectangular square-lattice on 20 qubits\n",
"num_rows = 3\n",
"num_cols = 3\n",
"num_qubits = num_rows * num_cols\n",
"\n",
"coupling_map = CouplingMap.from_grid(\n",
" num_rows=num_rows, num_columns=num_cols, bidirectional=False\n",
")\n",
"\n",
"# Create mapping from color to edge list\n",
"coloring = auto_color_edges(coupling_map.get_edges())\n",
"color_to_edge = defaultdict(list)\n",
"for edge, color in coloring.items():\n",
" color_to_edge[color].append(edge)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "adbceb93-6344-465c-9c4a-1c6e98ebc3d6",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"remove-input"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The circuit will have 9 qubits and 4 unique entangling layers.\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from rustworkx import PyDiGraph\n",
"from rustworkx.visualization import graphviz_draw\n",
"\n",
"# Inspect graph coupling and unique entangling layers\n",
"print(\n",
" f\"The circuit will have {num_qubits} qubits and {len(color_to_edge)} unique entangling layers.\"\n",
")\n",
"sq_lattice = PyDiGraph()\n",
"sq_lattice.extend_from_weighted_edge_list(\n",
" [\n",
" (source, target, color)\n",
" for ((source, target), color) in coloring.items()\n",
" ]\n",
")\n",
"\n",
"\n",
"def color_edge_4color(edge):\n",
" color_dict = {0: \"red\", 1: \"green\", 2: \"blue\", 3: \"orange\"}\n",
" return {\"color\": color_dict[edge]}\n",
"\n",
"\n",
"graphviz_draw(sq_lattice, edge_attr_fn=color_edge_4color, method=\"neato\")"
]
},
{
"cell_type": "markdown",
"id": "538ccd51-897b-4794-aaff-864d0cea02b3",
"metadata": {},
"source": [
"## Generate synthetic noise models\n",
"\n",
"Before creating the quantum circuits, we generate a noise model ([PauliLindbladError](https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.PauliLindbladError.html) instance) for each of the entangling layers. We will embed them as instructions in our quantum circuits later. For each layer, we generate noise channels of varying scales. Specifically, we generate noise models with [Error Per Layered Gate (EPLG)](https://www.ibm.com/quantum/blog/quantum-metric-layer-fidelity) of approximately `.0004, .0008, .0012, .0016,` and `.002`.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "3b3bb04e-d8f7-429f-b396-258ccd7197da",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"from qiskit.quantum_info import SparsePauliOp, pauli_basis\n",
"from qiskit_aer.noise import PauliLindbladError\n",
"\n",
"# Pauli-Lindblad noise parameters\n",
"seed = 1764\n",
"target_EPLGs = [0.0004, 0.0008, 0.0012, 0.0016, 0.002]\n",
"\n",
"\n",
"def generate_random_pauli_lindblad_noise(\n",
" edges,\n",
" num_qubits: int | None = None,\n",
" noise_scale: float = 1e-3,\n",
" seed: int | None = None,\n",
") -> PauliLindbladError:\n",
" \"\"\"Generate random Pauli-Lindblad noise over the full Pauli basis.\"\"\"\n",
" if num_qubits is None:\n",
" num_qubits = np.max(edges)\n",
"\n",
" basis_paulis = [p for p in pauli_basis(2) if np.sum(p.x + p.z)]\n",
" basis_paulis = SparsePauliOp.from_sparse_list(\n",
" [\n",
" (pauli.to_label(), edge, 1)\n",
" for pauli in basis_paulis\n",
" for edge in edges\n",
" ],\n",
" num_qubits=num_qubits,\n",
" )\n",
" basis_paulis = basis_paulis.simplify()\n",
" basis_paulis = basis_paulis.paulis\n",
"\n",
" rng = np.random.default_rng(seed=seed)\n",
" rates = rng.random(len(basis_paulis)) * noise_scale\n",
"\n",
" return PauliLindbladError(generators=basis_paulis, rates=rates)\n",
"\n",
"\n",
"num_generators = (\n",
" (num_rows * num_cols)\n",
" + (num_rows - 1) * num_cols\n",
" + num_rows * (num_cols - 1)\n",
")\n",
"noise_scales = [\n",
" EPLG * (num_rows * num_cols) / num_generators for EPLG in target_EPLGs\n",
"]\n",
"noise_models_per_EPLG = [\n",
" [\n",
" generate_random_pauli_lindblad_noise(\n",
" color_to_edge[color],\n",
" num_qubits=num_qubits,\n",
" noise_scale=noise_scale,\n",
" seed=seed,\n",
" )\n",
" for color in range(len(color_to_edge))\n",
" ]\n",
" for noise_scale in noise_scales\n",
"]"
]
},
{
"cell_type": "markdown",
"id": "5e2044c3-d265-4de3-b489-acfb4b4666de",
"metadata": {},
"source": [
"## Create the quantum circuits\n",
"\n",
"For this demo, we simulate the time dynamics of a transverse-field Ising model (TFIM) for increasing numbers of Trotter steps (1-10 steps). For each of the 10 circuit depths, we simulate the effect of gate noise, given noise models of varying scales (`EPLGs = .0004, .0008, .0012, .0016, .002`). The noise is inserted into the `QuantumCircuit` as a [PauliLindbladError](https://qiskit.github.io/qiskit-aer/stubs/qiskit_aer.noise.PauliLindbladError.html) instruction from Qiskit Aer. The Hamiltonian considered is:\n",
"\n",
"$H = -J\\sum\\limits_{\\langle i,j \\rangle} Z_iZ_j + h\\sum\\limits_iX_i$\n",
"\n",
"where $J>0$ describes the coupling of nearest-neighbor spins, $i"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from typing import Any\n",
"\n",
"from qiskit import QuantumCircuit\n",
"\n",
"# Ising model parameters\n",
"num_steps = 10\n",
"J = -1.0\n",
"dt = 0.25 / abs(J)\n",
"h = 2.0 * abs(J)\n",
"initial_state_angle = np.pi / 18.0\n",
"rx_angle = 2.0 * h * dt\n",
"rzz_angle = 2.0 * J * dt\n",
"\n",
"\n",
"def generate_ising_circuit(\n",
" num_qubits: int,\n",
" num_steps: int,\n",
" rx_angle: float,\n",
" rzz_angle: float,\n",
" coloring: dict[Any, list[tuple[int, int]]],\n",
" layer_noise_models: list[PauliLindbladError] | None = None,\n",
" initial_state_angle: float | None = None,\n",
") -> QuantumCircuit:\n",
" \"\"\"Generate a quantum circuit implementing a transverse-field Ising model\"\"\"\n",
" qc = QuantumCircuit(num_qubits)\n",
" if initial_state_angle:\n",
" qc.ry(initial_state_angle, range(num_qubits))\n",
" qc.rx(rx_angle / 2, range(num_qubits))\n",
" for i in range(num_steps):\n",
" for j, layer in enumerate(coloring):\n",
" edges = coloring[layer]\n",
" if layer_noise_models:\n",
" qc.append(layer_noise_models[j], qargs=range(num_qubits))\n",
" for edge in edges:\n",
" qc.rzz(rzz_angle, *edge)\n",
" if i == num_steps - 1:\n",
" qc.rx(rx_angle / 2, range(num_qubits))\n",
" else:\n",
" qc.rx(rx_angle, range(num_qubits))\n",
" return qc\n",
"\n",
"\n",
"# Create the noiseless and noisy circuits\n",
"noiseless_circs = []\n",
"noisy_circs = []\n",
"for steps in range(1, num_steps + 1):\n",
" noiseless_circs.append(\n",
" generate_ising_circuit(\n",
" num_qubits,\n",
" steps,\n",
" rx_angle,\n",
" rzz_angle,\n",
" color_to_edge,\n",
" initial_state_angle=initial_state_angle,\n",
" )\n",
" )\n",
" noisy_circs_per_step = []\n",
" for noise_models in noise_models_per_EPLG:\n",
" noisy_circs_per_step.append(\n",
" generate_ising_circuit(\n",
" num_qubits,\n",
" steps,\n",
" rx_angle,\n",
" rzz_angle,\n",
" color_to_edge,\n",
" layer_noise_models=noise_models,\n",
" initial_state_angle=initial_state_angle,\n",
" )\n",
" )\n",
" noisy_circs.append(noisy_circs_per_step)\n",
"print(\n",
" 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\"\n",
")\n",
"print(\"\\nBelow: Initial state and one noisy Trotter step.\")\n",
"noisy_circs[0][0].draw(\"mpl\", fold=-1)"
]
},
{
"cell_type": "markdown",
"id": "007a0560-e1bf-4715-a2fd-a8a5f1e5ac24",
"metadata": {},
"source": [
"## Specify observable and run simulations\n",
"\n",
"For this demo, we simulate expectation values of the average two-site correlator:\n",
"\n",
"$\\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$\n",
"\n",
"where $\\Psi(\\theta)$ corresponds to a uniform $R_y(\\theta)$ rotation on all qubits, $\\mathscr{U}^s$ describes $s$ Trotter layers, and $(j,k)$ index all connected pairs of vertices on the lattice.\n",
"\n",
"Finally, we use `pauli_prop` to simulate observable expectation values for each of the circuits. For this nine-qubit demo, we perform all simulations **exactly**. **No Pauli propagation truncation will be performed; therefore, the differences in expectation values across the different noise models can be entirely attributed to the gate error**. The simulation process is handled in four steps:\n",
"\n",
"* Evolve the Clifford gates in the circuit to the front of the circuit using `pauli_prop.evolve_through_cliffords`\n",
"* Propagate the observable through the non-Clifford part of the circuit using `pauli_prop.propagate_through_circuit`\n",
" * We perform exact simulations by allowing the observable to grow to the size of the full Pauli space, $4^9$\n",
"* Propagate the evolved observable through the Clifford part of the circuit using Qiskit's `SparsePauliOp.evolve`\n",
"* Estimate the expectation value with respect to the zero state, $|0\\rangle^{\\otimes N}$, by summing the coefficients of each diagonal term in $O$ (terms containing $Z$ or $I$ on all qubits)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "7ceb7937-ef7d-483a-b49a-d3aa2a3e6c43",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ran 10 noiseless and 50 noisy simulations in 103s.\n"
]
}
],
"source": [
"import time\n",
"\n",
"from pauli_prop import evolve_through_cliffords, propagate_through_circuit\n",
"from qiskit.quantum_info import Pauli\n",
"\n",
"# Average ZZ-correlator observable\n",
"id_pauli = Pauli(\"I\" * num_qubits)\n",
"observable = 2 * SparsePauliOp(\n",
" [\n",
" id_pauli.dot(Pauli(\"ZZ\"), [i, j])\n",
" for i in range(num_qubits)\n",
" for j in range(i + 1, num_qubits)\n",
" ]\n",
")\n",
"observable /= num_qubits**2\n",
"\n",
"# Pauli propagation parameters\n",
"max_terms = 4**num_qubits # Exact propagation\n",
"atol = 1e-12\n",
"\n",
"# Run simulations\n",
"exact_evs = []\n",
"noisy_evs = [[] for _ in range(len(target_EPLGs))]\n",
"st = time.perf_counter()\n",
"for i, noiseless_circ in enumerate(noiseless_circs):\n",
" cliff, non_cliff = evolve_through_cliffords(noiseless_circ)\n",
" evolved_obs = propagate_through_circuit(\n",
" observable, non_cliff, max_terms=max_terms, atol=atol, frame=\"h\"\n",
" )[0]\n",
" evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame=\"h\")\n",
" exact_evs.append(\n",
" float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())\n",
" )\n",
" for j in range(len(target_EPLGs)):\n",
" noisy_circ = noisy_circs[i][j]\n",
" cliff, non_cliff = evolve_through_cliffords(noisy_circ)\n",
" evolved_obs = propagate_through_circuit(\n",
" observable, non_cliff, max_terms=max_terms, atol=1e-12, frame=\"h\"\n",
" )[0]\n",
" evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame=\"h\")\n",
" noisy_evs[j].append(\n",
" float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())\n",
" )\n",
"print(\n",
" f\"Ran {len(noiseless_circs)} noiseless and {len(target_EPLGs) * num_steps} noisy simulations in {int(time.perf_counter() - st)}s.\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "d138fc46-5a47-46d8-9074-6aa30b273a86",
"metadata": {},
"source": [
"## Observe effect of gate error on the model\n",
"\n",
"Remember that, since this is a nine-qubit experiment, the Pauli propagation routine is exact, and all of the error in the noise plots can be attributed to gate error.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "73b00aa9-ee66-4b65-ad41-f6278555e5c5",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"remove-input"
]
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"xs = range(1, num_steps + 1)\n",
"plt.plot(xs, exact_evs, label=\"Noiseless\", color=\"black\", marker=\"o\")\n",
"colors = [\".3\", \".4\", \".5\", \".6\", \".7\"]\n",
"for i, evs in enumerate(noisy_evs):\n",
" plt.plot(\n",
" xs,\n",
" evs,\n",
" label=f\"{target_EPLGs[i]} EPLG\",\n",
" linestyle=\"--\",\n",
" color=colors[i],\n",
" marker=\"o\",\n",
" )\n",
"plt.xlabel(\"# Trotter steps\")\n",
"plt.ylabel(r\"$\\langle Z_{tot}^2 \\rangle$\")\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}