Quantum simulation underpins nearly every stage of quantum algorithm development — from prototyping circuits to validating results before hardware runs.
But once circuits move beyond ~20–30 qubits, no single classical method works universally. Statevector methods scale exponentially. Tensor networks like MPS depend strongly on entanglement structure and bond dimension. Sampling-based approaches trade accuracy for feasibility. GPU acceleration helps — but only when applied to the right workload.
The real bottleneck isn’t just scaling. It’s method selection and orchestration.
Today, developers are forced to manually choose simulators, tune parameters, manage hardware resources, and rewrite workflows when moving between simulation engines. The complexity of the tooling becomes the constraint.
At Qoro Quantum, we built Maestro to abstract the engineering complexity away — not the physics. One QuantumCircuit. One estimate() or execute() call, no matter the engine.
Maestro dynamically manages backend selection, CPU/GPU routing, multi-process and multi-thread execution, and parameter configuration behind a unified interface. It integrates multiple simulation methods from different providers and makes it straightforward to switch between them — or combine them — to optimize runtime, accuracy, and cost.
With a single line of code, a user can move from a CPU-based simulator to our GPU-accelerated MPS backend, reducing runtime from hours to minutes in the right regimes — and avoiding unnecessary GPU spend where acceleration offers no advantage.
To demonstrate this in practice, we’ve published an example tackling a real physics workload: strongly correlated electron dynamics in the Fermi–Hubbard model. The example shows where GPU acceleration delivers measurable speedups — and where intelligent backend orchestration prevents wasted compute on steps that don’t benefit from it.
Why GPU Matters for Quantum Simulation
Before diving into the example, it’s worth understanding why GPU acceleration matters specifically for MPS-based quantum simulation — and why it isn’t a blanket “everything runs faster” switch.
MPS simulation represents a quantum state as a chain of tensors, each with a bond dimension χ that controls the amount of entanglement the representation can capture. The core operation — contracting two tensors during a gate application — costs O(χ³). This means:
| Bond Dimension | Relative Cost | Typical Use |
| 16 | 1x | Quick exploration |
| 64 | 64x | Production accuracy |
| 256 | 4096x | Precision convergence |
At χ=16, tensor contractions are small matrix multiplications — a CPU handles them with negligible overhead. At χ=256, each contraction involves large dense matrices where GPU parallelism can provide 10–100× speedups over CPU. The data transfer overhead to the GPU (moving tensors to device memory and back) is amortized across the much larger compute workload.
This is why one should not blindly route everything to GPU. One should apply GPU where the workload justifies it — high bond dimensions, many repeated circuit evaluations — and stays on CPU for cheap steps. In code, this is a one-parameter change:
# CPU — fast for low bond dimension
result = qc.estimate(
simulator_type=maestro.SimulatorType.QCSim,
simulation_type=maestro.SimulationType.MatrixProductState,
max_bond_dimension=16,
observables=observables,
)
# GPU — essential for high bond dimension
result = qc.estimate(
simulator_type=maestro.SimulatorType.Gpu,
simulation_type=maestro.SimulationType.MatrixProductState,
max_bond_dimension=256,
observables=observables,
)Same circuit. Same observables. One parameter flips the backend.
Fermi-Hubbard Model: 40 Qubits from 200
What is it?
The Fermi-Hubbard model is one of the most important unsolved problems in condensed matter physics. It describes electrons hopping on a lattice with on-site Coulomb interactions:
Hopping term (t): electrons move between neighboring lattice sites.
Interaction term (U): two electrons occupying the same site pay an energy penalty.
The ratio U/t controls whether the system is metallic or insulating (the Mott transition).
Despite its simple structure — just hopping plus interaction — the model is believed to capture key mechanisms behind high-temperature superconductivity. Exact solutions exist only in one dimension, and even there, dynamics quickly become intractable. Simulating its real-time evolution is therefore a prime target for quantum computers.
The reason it makes an excellent benchmark: it maps exactly onto quantum hardware via the Jordan-Wigner transformation. Each electron spin state (occupied or empty) corresponds to a single qubit. A 100-site chain (100 sites × 2 spins) becomes a 200-qubit system. That’s genuinely at the frontier of classical tractability. We’re not using it because it’s a textbook example. We’re using it because it’s physically motivated, has known limiting behaviors to validate against, and sits precisely at the scale where the classical-to-quantum crossover happens.
Why does it matter computationally?
Physically relevant system sizes are in the hundreds of lattice sites. Classical exact diagonalization scales exponentially and becomes impossible well below that. Even tensor-network methods like MPS — which exploit 1D structure — become expensive when entanglement (bond dimension) grows.
But there is a crucial physical insight: The Lieb-Robinson Bound.
After a local quench — for example, preparing a domain wall state (left half filled, right half empty) — correlations do not spread instantaneously. Information propagates at a finite velocity, forming a causal light cone. Outside this light cone, the system remains essentially unchanged from its initial state.
For a 200-qubit chain evolved for short-to-moderate times:
- Only ~40 qubits near the domain wall exhibit non-trivial dynamics
- The remaining ~160 qubits remain effectively frozen
- Simulating all 200 uniformly wastes compute
The physics tells us that the effective dynamical system is much smaller than the nominal Hilbert space.
How Maestro runs it
Our Fermi-Hubbard example implements a 3-tier adaptive pipeline that exploits this structure:
Tier 1 — Scout
Engine: Pauli Propagator (Clifford proxy)
Scope: All 200 qubits
Hardware: CPU
We first construct a Clifford-only proxy of the real circuit and execute it using a Pauli propagation method. The proxy preserves circuit connectivity — which lattice sites influence which — while remaining efficiently simulable. We measure ⟨Zᵢ⟩ across all sites and compare against the initial domain-wall configuration. Sites whose expectation values deviate beyond a chosen threshold are marked as dynamically active.
Because this tier avoids non-Clifford operations and high bond dimensions, it completes in seconds, independent of the nominal system size. The result is a data-driven identification of the effective light cone.
This works because of the Lieb-Robinson bound. In a quantum lattice system with local interactions, information cannot propagate faster than a limit set by the coupling strength. After the release of the boundary wall, change starts to propagate. The Clifford circuit represents a worst case scenario for information propagation: it substitutes the actual parameterized rotation gates from the Trotterized time evolution, which are less than maximally entangling, by CNOT gates that maximally entangle sites.
Using a simple Maestro execution, we can quickly get the results:
result = scout_circuit.estimate(
observables=obs_list,
simulator_type=maestro.SimulatorType.QCSim,
simulation_type=maestro.SimulationType.PauliPropagator,
)This runs on the full 200-qubit system in under a second. We measure ⟨Z⟩ on each site and compare to the initial state. Sites that haven’t moved stay at ±1; sites in the active region deviate. The result: active region [38, 62), 46 qubits.
Tier 2 — Sniper
Engine: MPS (χ = 64)
Scope: ~40 qubits
Hardware: CPU
Next, we execute the real Trotterized circuit — including the non-Clifford Rz rotations arising from the Jordan–Wigner mapping — but only on the active subregion identified by the scout phase.
Instead of simulating 200 qubits, we simulate ~40.
This yields roughly a 5× reduction in simulated system width, dramatically lowering tensor contraction cost while retaining physically relevant dynamics. The moderate bond dimension keeps CPU execution efficient while capturing the dominant entanglement structure.
Tier 3 — Precision
Engine: MPS (χ = 256)
Scope: Same active region
Hardware: GPU
Finally, we re-run the same reduced subregion at higher bond dimension for converged, high-accuracy results.
This is where GPU acceleration becomes critical. Large-χ tensor contractions benefit strongly from parallel linear algebra, and the performance gap between CPU and GPU becomes substantial. χ=256 gives 4× more entanglement capacity than χ=64. Entanglement grows during evolution, and lower bond dimension truncates it. The GPU tier produces the production-quality result.
Because the active region has already been localized, GPU resources are applied only where they meaningfully improve accuracy — not wasted on frozen qubits outside the light cone.
Where GPU helps
The precision tier dominates total pipeline time, and this is by design — the scout and sniper tiers are deliberately cheap. At χ=256, each tensor contraction involves dense matrix operations that scale as O(χ³) = O(16.7 million) floating-point operations per contraction. GPU acceleration delivers a ~10× speedup on this tier.
Because the precision tier accounts for the majority of the total wall-clock time, GPU acceleration on this single step translates directly to a faster end-to-end result. The scout (Tier 1) remains on CPU — it’s already sub-second. The sniper (Tier 2) stays on CPU — at χ=64, the contractions aren’t large enough for GPU to outrun the data transfer overhead. GPU resources are applied precisely where the payoff is highest.
The Pattern: Scout, Snipe, Accelerate
This example follows a three-phase pattern:
- Scout cheaply: Pauli Propagator identifies where the interesting physics lives, at negligible cost.
- Snipe precisely: MPS runs the full simulation only on the region that matters, at moderate bond dimension.
- Accelerate the bottleneck: GPU resources are applied to the high-bond-dimension step that dominates runtime.
This isn’t a theoretical framework, it’s a practical workflow that runs through Maestro’s unified API and can be applied in other scenarios. The same Maestro QuantumCircuit flows through all three stages. Backend selection, bond dimension, and CPU/GPU routing are parameters, not architectural decisions. With very few lines of code, we can easily switch between methods, saving significant runtime by using the right tool for the job.
Try It
The example is available in our maestro-examples repository, each with interactive Jupyter notebooks:
pip install qoro-maestro
# Fermi-Hubbard: 200-qubit adaptive pipeline
python fermi_hubbard_demo.pyMaestro is open source and available on GitHub. For production workloads, it integrates with SLURM-managed HPC clusters via the Maestro SPANK plugin, or runs as a managed service through Qoro Cloud.
