QSimAdv: A Late-Bound, Vendor-Agnostic Architecture for High-Performance Quantum-Circuit Simulation

AI summary

Overview

Research area: High-performance quantum-circuit simulation — specifically the software architecture that connects quantum-simulation semantics (state vectors, stabiliser tableaux, density matrices) to vendor hardware libraries (NVIDIA CUDA/cuBLAS, AMD HIP/hipBLAS, CPU/SIMD, MPI).

Technical level: Advanced. The paper assumes familiarity with state-vector simulation, Clifford/stabiliser formalism, tensor-axis permutations, GEMM, and distributed-memory MPI decomposition. The architectural argument is accessible to an intermediate reader; the cost models and correctness invariants are not.

Scope (one sentence): The paper proposes and evaluates QSimAdv, a simulator whose vendor independence comes from late binding — deferring every representation, kernel, and data-placement decision until the inputs that constrain it are known — rather than from a shared cross-vendor kernel, and it benchmarks this design on NVIDIA GH200 and AMD MI250X/EPYC systems from one to 256 ranks.

What This Paper Is About

The fastest NVIDIA GPU simulators (Qiskit Aer's GPU backend, PennyLane Lightning-GPU, Google qsim-GPU, CUDA-Q) get their speed from NVIDIA's closed cuQuantum SDK, particularly the cuStateVec state-vector library. AMD's ROCm stack on the evaluated Setonix system provides HIP, hipBLAS/hipBLASLt, hipTensor, rocRAND and GPU-aware MPI, but no analogue of the gate, index-swap and measurement operations that cuStateVec exposes — so simply translating CUDA syntax to HIP leaves the simulator algorithms themselves unwritten.

QSimAdv's answer is to move the portability boundary upward: the layers above the backend speak only of logical qubits, tensor axes, matrices and layout maps, while each vendor backend keeps its own native libraries, kernels and feasibility limits. The goal is portability that does not cost native execution speed, and that extends from a single device to distributed memory.

Key Contributions

1. Semantic routing before allocation (R1). For generic dense and gpu requests that ask for sampled counts, the simulator inspects the circuit, noise model and requested output before allocating a 2^N state vector. Eligible Clifford circuits with Pauli noise are routed to a stabiliser tableau; explicit density-matrix, MPS and tensor-network requests keep their chosen representation and are never overridden.

2. A fused-operator interface shared across backends (R2). The CPU state-vector path, the CUDA and HIP GPU paths, rank-local execution inside MPI, and the density-matrix wrapper all accept the same ordered targets, operator and layout map. The generic lowering combines an axis permutation with GEMM, while diagonal, controlled and small-width gates use specialised kernels.

3. Distributed layout as explicit, inspectable state (R3). A proved bijective logical-to-physical map records non-canonical qubit order across both local index bits and rank-address bits, so the dispatcher performs a cross-rank exchange only when a gate actually targets a rank bit. The authors show that, under an ideal all-frame, full-future demand model with uniform promotion costs, this reduces to offline paging, where Bélády's rule minimises the number of promotions; the production planner uses bounded lookahead instead.

4. Backend lowering through a shared primitive contract (R4). CUDA, HIP, Arm, x86 and MPI implementations each supply four primitives — matrix multiplication, index permutation, sampling, and distributed exchange — plus specialised gate kernels, while preserving the same representation and layout contracts. Evaluation characterises the NVIDIA lowerings with a kernel decomposition and the AMD lowerings with a bounding host-region profile.

Main Findings

• Routing pays off enormously at the Clifford boundary: at N = 30, routing the noisy GHZ benchmark through the stabiliser path is about 78,000 times faster than general execution, and the tableau path reaches N = 1000 in under a second.

• Batched Pauli frames remove the per-shot trajectory loop: for terminal computational-basis readout, up to 64 independent Pauli frames are propagated together, one bit per 64-bit integer, so S shots need only ⌈S/64⌉ batched frame walks plus S noiseless samples. The construction is exact because inserting sampled Pauli errors into a Clifford circuit leaves each trajectory Clifford.

• Full-state wins at large N on Hopper: on the random-circuit sweep, QSimAdv overtakes Aer without cuStateVec from N = 28, reaching an Aer-to-QSimAdv time ratio of 2.2 at N = 32.

• Modest but consistent wins on AMD: on the four MI250X random-brickwork cells shared by all GPU tools, at N in {24, 26, 28, 30}, the Aer-to-QSimAdv ratio ranges from 1.15 to 2.16. With matched complex 32-bit floating-point state storage, the paper reports QSimAdv leading both Aer Hopper configurations at N = 32 and Aer's HIP backend at these four shared MI250X sizes.

• Strong scaling is platform-dependent: on Setonix, QSimAdv leads both the GPU and CPU comparisons at every measured rank, with 3.4× and 2.8× speedups respectively from one to eight ranks — while neither GH200 path speeds up at eight ranks.

• Weak scaling extends to distributed memory: the distributed interface reaches 256 ranks with 2 TiB GPU and 1 TiB CPU states.

• The diagonal-aware fusion score is not what produces the diagonal gain: changing the dynamic-programming score to be diagonal-aware moves the measured Hopper intervals by less than 2%; the benefit comes from the streaming elementwise kernel, not the partition heuristic.

• Lazy permutation is free in the accounting: the trailing inverse transpose of the generic lowering is never executed — it is recorded in the layout map instead — so a materialised permutation moves Θ(2^N) complex amplitudes (the whole state, once), dense GEMM costs Θ(2^(N+k)) complex multiply-adds, and a diagonal fused operator collapses to Θ(2^N) streaming work.

Methodology in Plain English

The design principle is "decide late, decide once, and decide with the information you need."

A run starts with a parsed OpenQASM 3.0 circuit, a noise model, a requested entry point and a shot count. Parsing first throws away operations outside the causal cone of the measured qubits. Then, if the request only asked for sampled counts through a generic entry point, the router checks whether every operation — gates, rotations equivalent to Clifford gates, Pauli measurements, and supported Pauli noise — stays inside Clifford dynamics. If so, an Aaronson–Gottesman tableau is used and no exponential state vector is ever allocated; the tableau is CPU-resident and its arrays of 64-bit integers are updated with AVX-512, AVX2, SVE/SVE2, NEON or scalar kernels.

If the full-state path is retained, the backend family is chosen next, which fixes the feasibility limits (fusion cap, distributed chunk width). Only then does fusion run: gates are reordered across commuting or disjoint neighbours, and a deterministic dynamic program over windows of at most 64 gates chooses block boundaries, scoring a dense width-k block as 4^k and a diagonal block as 1. Each fused block then presents its ordered logical targets to whichever backend is active, which resolves them to physical bit positions through the layout map and picks a concrete kernel or, for cross-rank targets, a transport.

For distributed execution, the state is viewed as an N-axis tensor where the r high bits select the MPI rank among P = 2^r ranks and Q_loc = N − r bits are local. A gate targeting a rank bit triggers a promotion: an exchange of that bit with a free local bit, recorded in the map so a later reverse exchange is unnecessary. The layout map is kept non-canonical until something — an output routine, for instance — explicitly demands canonical order.

Hardware: HPC Ella (Grace Hopper GH200, CUDA, HPC-X OpenMPI) and HPC Setonix (MI250X/EPYC, ROCm, Cray MPICH).

Why This Matters

Impact on research. The paper argues that portability need not begin at the kernel. If the semantics of a fused operator and the logical-to-physical layout live above the vendor API, then a simulator can adopt a new accelerator by writing four primitives rather than reimplementing the algorithms. That matters for reproducibility: the same logical semantics are preserved across NVIDIA, AMD, Arm, x86 and MPI, so a result obtained on one vendor's stack is, in the intended sense, the same computation on another's. The explicit, provably bijective layout map also turns an informal optimisation (deferring data movement) into inspectable state that can be reasoned about, and connects distributed qubit placement to a well-studied optimality theory (offline paging and Bélády's rule).

Real-world applications:

• Cross-vendor verification of quantum algorithms: running the same circuit and noise model on NVIDIA and AMD hardware without rewriting code, to check that a discrepancy is physics rather than a porting artefact.
• Noisy intermediate-scale benchmarking: the noisy GHZ routing result makes it practical to sample Clifford-with-Pauli-noise circuits at N = 30 and up to N = 1000, which is where error-correction and noise-modelling studies live.
• Procurement and portability planning for national HPC facilities, which increasingly host mixed NVIDIA/AMD fleets and cannot afford simulator lock-in.
• Large-state reference data: weak scaling to 256 ranks with 2 TiB GPU and 1 TiB CPU states gives a route to full-state references beyond single-device memory for validating quantum hardware and approximate methods.

Industry relevance. Supercomputing centres, quantum-hardware vendors and cloud providers all face the same asymmetry: one vendor ships a quantum-simulation library, the other ships only a general-purpose toolkit. QSimAdv's separation — logical semantics above, four primitives below — is a template any accelerator vendor could implement against, and its strong-scaling result on Setonix (3.4× GPU, 2.8× CPU from one to eight ranks, against no speedup at eight ranks on either GH200 path) is a concrete data point that vendor choice affects distributed scaling behaviour, not just single-device throughput.

Future Directions

• Broadening the tableau-eligible noise set. Amplitude damping and generic custom Kraus channels currently keep a request on the full-state path, and the router deliberately does not inspect a custom Kraus list to infer that it happens to define a Pauli mixture. Extending eligibility — or safely detecting Pauli-mixture structure in user-supplied channels — would widen the fast path.

• Removing the terminal-readout restriction. The batched Pauli-frame fast path assumes one terminal joint computational-basis readout, with all sampled noise preceding the first recorded measurement and no reset, adaptive classical control or mid-circuit dependence on a measurement outcome. Circuits outside that regime fall back to the per-shot Algorithm 1, so mid-circuit measurement and feed-forward are the natural next target.

• GPU-resident tableaux. The stabiliser tableau is CPU-resident even when the request came through a CUDA or HIP entry point (and therefore allocates no GPU state vector). Moving the tableau to the accelerator, and reconciling that with the AVX-512/AVX2/SVE/NEON vectorisation already in place, is an open engineering question.

• From bounded lookahead to stronger placement policies. The paper proves that an ideal all-frame, full-future demand model with uniform promotion costs reduces to offline paging, where Bélády's rule is optimal, but the production planner uses bounded lookahead. Closing that gap — with better lookahead, or with the non-uniform promotion costs the model abstracts away — is a clear direction, as is characterising the MPS engine and the NVIDIA-only cuTensorNet path, which the paper includes but does not evaluate.

Target Audience

Quantum-computing software architects and HPC performance engineers are the primary audience, particularly those maintaining or extending a simulator that must run on more than one accelerator vendor. It is also valuable to facility staff who plan mixed NVIDIA/AMD deployments and need evidence about distributed scaling behaviour, and to researchers who need large full-state or Clifford-with-noise reference data at N = 30 and beyond. Readers wanting only the algorithmic results can take the routing and Pauli-frame findings without the cost models; readers wanting the architecture can take R1–R4 without the benchmarks. A beginner will find the introduction and §II readable but will need prior exposure to state-vector simulation and MPI to follow §III–IV.

Authors’ abstract

Portability in high-performance quantum-circuit simulation need not begin at the kernel. We present QSimAdv, which makes late binding, rather than a common kernel, the basis of vendor independence. Representation, operator lowering, and data placement are bound only when their required inputs become available. Before full-state allocation, circuit, noise, and output inspection can route eligible generic sampled-count requests to a stabiliser tableau; explicitly requested representations remain fixed. For full-state execution, backend constraints shape fusion; an ordered fused operator binds to a native lowering only after its physical targets are known. A first-class logical-to-physical layout map records non-canonical order across local and rank-address bits, so the dispatcher moves nonlocal targets only on demand. GPU, CPU, and Message Passing Interface (MPI) backends share these semantics while retaining native execution paths. We realize this design on NVIDIA GH200 and AMD MI250X/EPYC systems across local and distributed execution. With matched complex 32-bit floating-point state storage, QSimAdv leads both Aer Hopper configurations at $N=32$ and Aer's HIP backend at four shared MI250X sizes from $N=24$ to 30. Strong scaling exposes platform dependence: on setonix, QSimAdv leads both GPU and CPU comparisons at every measured rank, achieving $3.4\times$ and $2.8\times$ speedups, respectively, from one to eight ranks; neither the GH200 path nor the CPU path speeds up at eight ranks. Weak scaling reaches 256 ranks with 2 TiB GPU and 1 TiB CPU states. Together, these results support that portability can reside above the kernel boundary while execution remains native and extends across distributed memory.