Noisy Circuit Simulation¶
Strong digital simulation evolves a matrix-product state (MPS) through a Qiskit circuit while applying tensor-jump noise between gates. You choose which Pauli (or other) observables to measure and how many Monte Carlo trajectories to average.
Workflow |
Typical use |
Key settings |
|---|---|---|
Final observables |
Noise scaling, benchmarking, device studies |
|
Mid-circuit observables |
Layer-wise diagnostics, depth-dependent calibration |
|
Shot-based readout |
Hardware-like bitstring statistics |
Circuits enter YAQS as qiskit.circuit.QuantumCircuit objects (or OpenQASM strings). The initial state must use representation="mps" (the default for State presets). For accuracy presets, truncation knobs, and random_seed, see Configuring Simulation Parameters. For bell-curve (Gaussian) noise strengths, see Realistic Noise Models.
1from mqt.yaqs import Simulator
2
3sim = Simulator(show_progress=False)
1. Final observables and noise scaling¶
We build an Ising-style circuit, prepare \(\ket{0}^{\otimes n}\), and sweep a global relaxation rate \(\gamma\). Expectation values \(\langle Z_i \rangle\) at the end of the circuit are collected into a heatmap.
1.1 Circuit and initial state¶
1from mqt.yaqs import State
2from mqt.yaqs.core.libraries.circuit_library import create_ising_circuit
3
4num_qubits = 5
5circuit = create_ising_circuit(L=num_qubits, J=1, g=0.5, dt=0.1, timesteps=10)
6circuit.measure_all()
7circuit.draw(output="mpl")
8
9state = State(num_qubits, initial="zeros")
1.2 Simulation parameters¶
1from mqt.yaqs import Observable, StrongSimParams
2
3sim_params = StrongSimParams(
4 observables=[Observable("z", site) for site in range(num_qubits)],
5 num_traj=100,
6 max_bond_dim=4,
7 svd_threshold=1e-6,
8)
1.3 Noise sweep and heatmap¶
1import numpy as np
2
3from mqt.yaqs import NoiseModel
4
5gammas = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1]
6heatmap = np.empty((num_qubits, len(gammas)))
7for j, gamma in enumerate(gammas):
8 noise_model = NoiseModel([
9 {"name": "lowering", "sites": [i], "strength": gamma} for i in range(num_qubits)
10 ])
11 result = sim.run(state, circuit, sim_params, noise_model)
12 for i in range(len(result.observables)):
13 heatmap[i, j] = result.expectation_values[i][0]
1import matplotlib.pyplot as plt
2
3fig, ax = plt.subplots(figsize=(6, 3.5))
4im = ax.imshow(heatmap, aspect="auto", vmin=0.5, vmax=1)
5ax.set_ylabel("Qubit")
6ax.set_xlabel(r"Relaxation rate $\gamma$")
7ax.set_yticks(range(num_qubits))
8ax.set_xticks(range(len(gammas)))
9ax.set_xticklabels([f"$10^{{{int(np.log10(g))}}}$" for g in gammas])
10
11fig.subplots_adjust(top=0.95, right=0.88)
12cbar_ax = fig.add_axes(rect=(0.9, 0.11, 0.025, 0.8))
13cbar = fig.colorbar(im, cax=cbar_ax)
14cbar.ax.set_title(r"$\langle Z \rangle$")
15plt.tight_layout()
16plt.show()
/tmp/ipykernel_2446/2494117714.py:15: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
plt.tight_layout()
2. Mid-circuit observables¶
Note
This section uses num_traj=1000 and may take one to two minutes during a documentation build. Reduce num_traj for faster local runs.
Set sample_layers=True on StrongSimParams and insert barriers labelled SAMPLE_OBSERVABLES (case-insensitive) where you want measurements. YAQS records observables at the circuit start, after each labelled barrier, and after the final gate layer.
Other barriers and measure operations are ignored for this sampling schedule.
2.1 Circuit with sampling barriers¶
1from qiskit.circuit import QuantumCircuit
2
3layer_qubits = 3
4qc = QuantumCircuit(layer_qubits)
5
6for segment in range(5):
7 qc.rzz(0.5, 0, 1)
8 qc.rzz(0.5, 1, 2)
9 if segment < 4:
10 qc.barrier(label="SAMPLE_OBSERVABLES")
11
12qc.draw(output="mpl")
Five entangler layers with four labelled barriers yield six sampling points: initial state, four mid-circuit checkpoints, and the final layer.
2.2 Noise model and parameters¶
1import numpy as np
2
3from mqt.yaqs import Observable, StrongSimParams
4
5noise_factor = 0.01
6layer_noise = NoiseModel(
7 [{"name": "pauli_x", "sites": [i], "strength": noise_factor} for i in range(layer_qubits)]
8 + [{"name": "crosstalk_xx", "sites": [i, i + 1], "strength": noise_factor} for i in range(layer_qubits - 1)]
9)
10
11layer_state = State(layer_qubits, initial="zeros", pad=2)
12layer_observables = [Observable("z", i) for i in range(layer_qubits)]
13layer_params = StrongSimParams(layer_observables, num_traj=1000, sample_layers=True)
Higher num_traj reduces Monte Carlo variance; 100 trajectories are often enough for prototyping.
2.3 Run and validate¶
1layer_result = sim.run(layer_state, qc, layer_params, layer_noise)
2
3# Reference expectations (rows: qubits; columns: initial + 4 barriers + final)
4reference = np.array([
5 [1.0, 0.9607894391523233, 0.9231163463866354, 0.8869204367171571, 0.8521437889662108, 0.8187307530779814],
6 [1.0, 0.9231163463866359, 0.8521437889662113, 0.7866278610665535, 0.726149037073691, 0.6703200460356394],
7 [1.0, 0.9607894391523233, 0.9231163463866354, 0.8869204367171571, 0.8521437889662108, 0.8187307530779814],
8])
9
10yaqs = np.vstack([np.real(v) for v in layer_result.expectation_values])
11max_diff = float(np.abs(yaqs - reference).max())
12print(f"max|YAQS − reference| = {max_diff:.4f} (decreases with num_traj)")
1import matplotlib.pyplot as plt
2
3layers = range(yaqs.shape[1])
4labels = [f"q{i}" for i in range(layer_qubits)]
5
6fig, ax = plt.subplots(figsize=(7, 4))
7for q in range(layer_qubits):
8 ax.plot(layers, reference[q], "--", color=f"C{q}", alpha=0.6, label=f"reference {labels[q]}")
9 ax.plot(layers, yaqs[q], "o-", color=f"C{q}", label=f"YAQS {labels[q]}")
10
11ax.set_xlabel("Sampling point (initial → barriers → final)")
12ax.set_ylabel(r"$\langle Z \rangle$")
13ax.set_xticks(list(layers))
14ax.legend(ncols=2, fontsize=8)
15ax.grid(alpha=0.3)
16plt.tight_layout()
17plt.show()
3. OpenQASM inputs¶
Pass an OpenQASM 2 source string (or file path) directly to run() instead of building a qiskit.circuit.QuantumCircuit in Python. Custom gate bodies declared in the program are translated like any other Qiskit operation.
1from mqt.yaqs import State, WeakSimParams
2
3qasm = """
4OPENQASM 2.0;
5include "qelib1.inc";
6
7gate entangle a,b {
8 h a;
9 cx a,b;
10}
11
12qreg q[2];
13entangle q[0], q[1];
14"""
15
16qasm_state = State(2, initial="zeros")
17qasm_result = sim.run(
18 qasm_state,
19 qasm,
20 WeakSimParams(shots=128, max_bond_dim=4),
21)
22print(f"Top bitstrings: {list(qasm_result.counts.items())[:3]}")
Measuring shots: 0%| | 0/128 [00:00<?, ?it/s]
Measuring shots: 77%|██████████████████▌ | 99/128 [00:00<00:00, 984.59it/s]
Measuring shots: 100%|██████████████████████| 128/128 [00:00<00:00, 1019.22it/s]
Top bitstrings: [(0, 58), (3, 70)]
OpenQASM 3 requires pip install mqt-yaqs[qasm3]. EquivalenceChecker accepts the same path and string forms; see Equivalence Checking.
4. Gate application modes¶
StrongSimParams.gate_mode (and WeakSimParams.gate_mode) selects how two-qubit gates are applied to the MPS. The default "mpo" uses extended gate MPOs for long-range pairs; "tdvp" uses a local TDVP window when an analytic generator is available. See Configuring Simulation Parameters and Custom and Qiskit Gates in YAQS for the full matrix.
Below, a long-range cx on qubits 0 and 2 is simulated noiselessly with both modes:
1from qiskit.circuit import QuantumCircuit
2
3lr_qc = QuantumCircuit(3)
4lr_qc.h(0)
5lr_qc.cx(0, 2)
6
7lr_state = State(3, initial="zeros")
8for mode in ("mpo", "tdvp"):
9 mode_params = StrongSimParams(
10 observables=[Observable("z", 0)],
11 num_traj=1,
12 gate_mode=mode,
13 max_bond_dim=8,
14 )
15 mode_result = sim.run(lr_state, lr_qc, mode_params)
16 z0 = float(mode_result.expectation_values[0][0])
17 print(f"gate_mode={mode!r}: ⟨Z₀⟩ = {z0:.6f}")
gate_mode='mpo': ⟨Z₀⟩ = 0.000000
gate_mode='tdvp': ⟨Z₀⟩ = 0.000000
/tmp/ipykernel_2446/2276085960.py:16: ComplexWarning: Casting complex values to real discards the imaginary part
z0 = float(mode_result.expectation_values[0][0])