Strong Simulation¶
Strong digital simulation evolves a matrix-product state (MPS) through a Qiskit circuit and evaluates Pauli (or custom) observables. Pass an optional NoiseModel as the fourth argument to run() for open-system tensor-jump trajectories; omit it for a single unitary path (regardless of num_traj).
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 should use representation="mps" (the default for State presets). For accuracy presets, truncation knobs, and random_seed, see Configuring Simulation Parameters. For log-normal disorder on noise strengths, see Realistic Noise Models.
1import matplotlib.pyplot as plt
2import numpy as np
3
4from mqt.yaqs import Simulator
5
6sim = Simulator(show_progress=False)
1. Minimal run: unitary vs open-system noise¶
Evolve a short Trotterized Ising circuit and compare final \(\langle Z_i\rangle\) without noise and with on-site amplitude damping:
1from mqt.yaqs import NoiseModel, Observable, State, StrongSimParams
2from mqt.yaqs.core.libraries.circuit_library import create_ising_circuit
3
4num_qubits = 3
5qc = create_ising_circuit(L=num_qubits, J=1.0, g=0.8, dt=0.1, timesteps=6)
6circuit_state = State(num_qubits, initial="zeros")
7circuit_params = StrongSimParams(
8 observables=[Observable("z", site) for site in range(num_qubits)],
9 preset="fast",
10 num_traj=32,
11)
12noise_model = NoiseModel([
13 {"name": "lowering", "sites": [site], "strength": 0.05} for site in range(num_qubits)
14])
15
16clean_result = sim.run(circuit_state, qc, circuit_params)
17noisy_result = sim.run(State(num_qubits, initial="zeros"), qc, circuit_params, noise_model)
18clean_z = np.array([float(np.real(v[0])) for v in clean_result.expectation_values])
19noisy_z = np.array([float(np.real(v[0])) for v in noisy_result.expectation_values])
20
21fig, ax = plt.subplots(figsize=(5, 3), layout="constrained")
22x = np.arange(num_qubits)
23bar_width = 0.35
24ax.bar(x - bar_width / 2, clean_z, bar_width, label="unitary", color="0.55")
25ax.bar(x + bar_width / 2, noisy_z, bar_width, label="with damping", color="C0")
26ax.set_xticks(x, [rf"$\langle Z_{i}\rangle$" for i in range(num_qubits)])
27ax.set_ylim(-1.05, 1.05)
28ax.set_ylabel("expectation value")
29ax.set_title("Optional noise model on the fourth `run` argument")
30ax.legend(frameon=False)
<matplotlib.legend.Legend at 0x7078df65d550>
2. Noise-strength sweep¶
On a longer chain, sweep a global relaxation rate \(\gamma\) and track how each qubit’s final \(\langle Z_i \rangle\) moves toward \(+1\) as damping dominates:
1num_qubits = 5
2circuit = create_ising_circuit(L=num_qubits, J=1.0, g=0.5, dt=0.1, timesteps=10)
3state = State(num_qubits, initial="zeros")
4sim_params = StrongSimParams(
5 observables=[Observable("z", site) for site in range(num_qubits)],
6 num_traj=64,
7 max_bond_dim=8,
8 svd_threshold=1e-6,
9)
10
11gammas = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1]
12heatmap = np.empty((num_qubits, len(gammas)))
13for j, gamma in enumerate(gammas):
14 damping = NoiseModel([
15 {"name": "lowering", "sites": [site], "strength": gamma} for site in range(num_qubits)
16 ])
17 result = sim.run(state, circuit, sim_params, damping)
18 for i in range(num_qubits):
19 heatmap[i, j] = float(np.real(result.expectation_values[i][0]))
20
21fig, ax = plt.subplots(figsize=(7, 4), layout="constrained")
22colors = plt.cm.viridis(np.linspace(0.15, 0.85, num_qubits))
23for i in range(num_qubits):
24 ax.semilogx(gammas, heatmap[i], "o-", color=colors[i], linewidth=1.8, markersize=5, label=rf"$q_{i}$")
25ax.set_xlabel(r"Relaxation rate $\gamma$")
26ax.set_ylabel(r"$\langle Z_i \rangle$")
27ax.set_ylim(-1.05, 1.05)
28ax.legend(ncol=num_qubits, fontsize=8, loc="lower left", frameon=False)
29ax.set_title("Final magnetization vs damping strength")
30ax.grid(alpha=0.3, which="both")
3. Mid-circuit observables¶
Note
This section uses num_traj=64 during the documentation build. Increase num_traj locally for lower-variance layer curves.
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.
The example below starts from \(\ket{+}^{\otimes n}\), applies a chain of \(R_{ZZ}\) entanglers, and tracks how amplitude damping gradually drives each \(\langle Z_i \rangle\) toward \(+1\). Only barriers labelled SAMPLE_OBSERVABLES trigger sampling; unlabelled barriers are ignored.
1from qiskit.circuit import QuantumCircuit
2
3layer_qubits = 5
4qc = QuantumCircuit(layer_qubits)
5
6for segment in range(6):
7 for i in range(layer_qubits - 1):
8 qc.rzz(0.7, i, i + 1)
9 if segment < 5:
10 qc.barrier(label="SAMPLE_OBSERVABLES")
11
12noise_factor = 0.1
13layer_noise = NoiseModel([
14 {"name": "lowering", "sites": [i], "strength": noise_factor} for i in range(layer_qubits)
15])
16
17layer_state = State(layer_qubits, initial="x+", pad=16)
18layer_params = StrongSimParams(
19 observables=[Observable("z", i) for i in range(layer_qubits)],
20 num_traj=64,
21 sample_layers=True,
22 max_bond_dim=12,
23)
24
25layer_result = sim.run(layer_state, qc, layer_params, layer_noise)
26layer_traj = np.vstack([np.real(v) for v in layer_result.expectation_values])
27
28fig, ax = plt.subplots(figsize=(8, 4), layout="constrained")
29depth = np.arange(layer_traj.shape[1])
30qubit_labels = [rf"$q_{i}$" for i in range(layer_qubits)]
31im = ax.imshow(
32 layer_traj,
33 aspect="auto",
34 origin="lower",
35 vmin=-1,
36 vmax=1,
37 extent=(-0.5, layer_traj.shape[1] - 0.5, -0.5, layer_qubits - 0.5),
38)
39ax.set_xlabel("Sampling index")
40ax.set_ylabel("Qubit")
41ax.set_xticks(depth)
42ax.set_yticks(range(layer_qubits), qubit_labels)
43ax.set_title(r"Mid-circuit $\langle Z \rangle$ under damping")
44fig.colorbar(im, ax=ax, shrink=0.9, label=r"$\langle Z \rangle$")
<matplotlib.colorbar.Colorbar at 0x7078dca43a10>
4. 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 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)
Measuring shots: 0%| | 0/128 [00:00<?, ?it/s]
Measuring shots: 62%|██████████████▊ | 79/128 [00:00<00:00, 789.82it/s]
Measuring shots: 100%|███████████████████████| 128/128 [00:00<00:00, 987.23it/s]
OpenQASM 3 requires pip install mqt-yaqs[qasm3]. EquivalenceChecker accepts the same path and string forms; see Equivalence Checking.
5. 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 Gates in YAQS for the full matrix.
Below, a long-range cx on qubits 0 and 2 is simulated noiselessly with both modes:
1lr_qc = QuantumCircuit(3)
2lr_qc.h(0)
3lr_qc.cx(0, 2)
4
5lr_state = State(3, initial="zeros")
6z0_by_mode = {}
7for mode in ("mpo", "tdvp"):
8 mode_params = StrongSimParams(
9 observables=[Observable("z", 0)],
10 num_traj=1,
11 gate_mode=mode,
12 max_bond_dim=8,
13 )
14 mode_result = sim.run(lr_state, lr_qc, mode_params)
15 z0_by_mode[mode] = float(np.real(mode_result.expectation_values[0][0]))
16
17print({mode: round(value, 4) for mode, value in z0_by_mode.items()})
{'mpo': 0.0, 'tdvp': 0.0}