Representation Comparison

YAQS supports multiple state representations for analog evolution. Each path targets a different scaling regime; the table below summarizes when each is appropriate.

For how to set representation on State, see Initializing Quantum States.

Choosing a representation

Path

When to use

Notes

"mps" (default)

Larger systems and tensor-network-friendly Hamiltonians

TJM trajectories; tune num_traj, max_bond_dim, and accuracy presets

"vector"

MCWF / state-vector quantum trajectories

Exponential memory in qubits; single-trajectory wavefunction dynamics

"density_matrix"

Lindblad master-equation evolution

Exponential memory; deterministic ensemble average without trajectory sampling

Practical guidance:

  • Start with preset="balanced" (or "fast" while exploring) on AnalogSimParams and increase num_traj until observables stabilize.

  • Tighten max_bond_dim / svd_threshold when entanglement growth demands it.

  • For trade-offs between unravellings and trajectory cost, see [] (references).

The sections below run the same noisy benchmark on all three paths so you can validate agreement on small systems. We use a product \(|{+}\rangle^{\otimes L}\) initial state so the MPS and dense backends encode the same physical state (Haar-random MPS states can disagree with unfolded vectors on non-local observables).

1. Noisy open-system benchmark

 1import matplotlib.pyplot as plt
 2import numpy as np
 3
 4from mqt.yaqs import AnalogSimParams, Hamiltonian, NoiseModel, Observable, Simulator, State
 5
 6sim = Simulator(show_progress=False)
 7
 8L = 3
 9H = Hamiltonian.ising(L, J=1.0, g=0.5)
10noise = NoiseModel([{"name": "pauli_z", "sites": [i], "strength": 0.2} for i in range(L)])
11obs = Observable("x", sites=[0])
12
13init_ref = State(L, initial="x+")
14psi0 = init_ref.mps.to_vec()
15rho0 = np.outer(psi0, psi0.conj())
16mps_tensors0 = [np.asarray(t, dtype=np.complex128).copy() for t in init_ref.mps.tensors]
17
18# Doc-build-friendly settings; increase t_max / num_traj for production runs.
19t_max = 1.0
20dt = 0.1
21num_traj = 32
22seed = 7
23
24params_rho = AnalogSimParams(observables=[obs], elapsed_time=t_max, dt=dt)
25result_rho = sim.run(State(density_matrix=rho0), H, params_rho, noise)
26res_rho = result_rho.expectation_values[0].flatten()
27times = params_rho.times
28
29params_vector = AnalogSimParams(
30    observables=[obs], elapsed_time=t_max, dt=dt, num_traj=num_traj, random_seed=seed,
31)
32result_vector = sim.run(State(vector=psi0), H, params_vector, noise)
33res_vector = result_vector.expectation_values[0].flatten()
34
35params_mps = AnalogSimParams(
36    observables=[obs], elapsed_time=t_max, dt=dt, num_traj=num_traj, max_bond_dim=16, random_seed=seed,
37)
38result_mps = sim.run(State(L, tensors=[t.copy() for t in mps_tensors0]), H, params_mps, noise)
39res_mps = result_mps.expectation_values[0].flatten()

 1fig, ax = plt.subplots(figsize=(6, 3.5), layout="constrained")
 2ax.plot(times, res_rho, label="density_matrix (exact)", linewidth=2, color="black")
 3ax.plot(times, res_vector, label=f"vector ({num_traj} traj)", linestyle="--")
 4ax.plot(times, res_mps, label=f"mps ({num_traj} traj)", linestyle=":")
 5ax.set_xlabel("Time")
 6ax.set_ylabel(r"$\langle X_0 \rangle$")
 7ax.set_ylim(-1.05, 1.05)
 8ax.legend()
 9ax.set_title(r"Open-system evolution across representations ($|+\rangle^{\otimes L}$ init)")
10ax.grid(alpha=0.3)
11plt.show()
../_images/0b47b1f14d9e2897a6f1a6d03feeb06cfa591ddab71f69eebf642240d81d899d.svg

Note

vector and mps curves are Monte Carlo means over num_traj trajectories; statistical error scales as \(1/\sqrt{N_{\mathrm{traj}}}\). The density_matrix path returns the deterministic ensemble average directly. Increase num_traj and t_max for smoother curves.

2. Noiseless cross-check

With noise_model=None, all three representations should agree on unitary observables (single trajectory for mps and vector). We use a product state here so the MPS path is exact at modest bond dimension.

 1obs_z = Observable("z", sites=[0])
 2params_mps_u = AnalogSimParams(observables=[obs_z], elapsed_time=0.5, dt=0.1, max_bond_dim=16)
 3params_rho_u = AnalogSimParams(observables=[obs_z], elapsed_time=0.5, dt=0.1)
 4
 5init_product = State(L, initial="x+")
 6psi_prod = init_product.mps.to_vec()
 7rho_prod = np.outer(psi_prod, psi_prod.conj())
 8mps_prod = [np.asarray(t, dtype=np.complex128).copy() for t in init_product.mps.tensors]
 9
10z_mps = sim.run(State(L, tensors=mps_prod), H, params_mps_u, None).expectation_values[0][-1]
11z_vec = sim.run(State(vector=psi_prod), H, params_mps_u, None).expectation_values[0][-1]
12z_rho = sim.run(State(density_matrix=rho_prod), H, params_rho_u, None).expectation_values[0][-1]