BQSim: GPU-accelerated Batch Quantum Circuit Simulation using Decision Diagram

Simulating a quantum circuit on an ordinary computer is how researchers check that a quantum algorithm does what they think, and every simulator is built to push one input through one circuit as fast as possible — but testing and verification need hundreds of different inputs through the same circuit, which is a different shape of problem with data parallelism a graphics processor could exploit; BQSim (Chinese University of Hong Kong and Wisconsin-Madison) rebuilds the pipeline around that, converting a compressed matrix representation into a form a GPU handles well and describing the whole simulation as a task graph.

The argument in one glance

Problem

Checking a quantum algorithm means simulating it on an ordinary computer

Quantum hardware is scarce and expensive, so researchers test their algorithms by simulating them classically. That is costly by construction: representing the state of a machine with n quantum bits takes two to the power of n numbers, and serious circuits contain thousands of operations.
What people did

Compress the matrices and fuse the operations

The matrices describing quantum operations are full of zeros and repeated blocks, so a decision diagram can store them as a small graph rather than a large array. That compression also lets many consecutive operations be merged into one, which is the main way simulators reduce work.
Why it fell short

All of it is built for one input at a time

Testing and verifying an algorithm means running many different inputs through the same circuit to see how it behaves in different scenarios. Every existing simulator optimises the single-input case, so the usual answer is to launch one process per input — trivially parallel, and blind to the fact that all those inputs are doing identical work.
This paper's idea

Treat the batch as the unit of work and move it to a graphics processor

Many inputs through one circuit is a large, regular, data-parallel computation. Price each operation by what it will cost in that setting, convert the compressed representation into a layout the hardware handles well, and describe the whole simulation as a graph of tasks so the runtime can overlap transfers with computation.

Why this problem is worth solving

The cost of the thing being simulated grows in the worst possible way.

Computing what a quantum algorithm would produce, using ordinary hardware. The state of a machine with n quantum bits is described by two to the power of n complex numbers, and applying an operation means multiplying that vector by a matrix. Full-state simulation of this kind gives the complete set of ideal amplitudes after every operation, which is exactly what makes it useful for understanding an algorithm and exactly what makes it expensive. A circuit with 21 quantum bits already requires over two million amplitudes to be updated by each of more than a thousand gates, and that is before anyone asks for it to be done for many different starting states.

is why so much work goes into avoiding the arithmetic rather than doing it faster.

The structure that makes it tractable is a compression scheme with a useful side effect.

A way of storing a matrix as a small graph instead of a large array. Each node stands for a sub-matrix, and its outgoing edges point to the four quarters of that sub-matrix, each edge carrying a weight. Wherever two quarters are identical, they share a node; wherever a quarter is entirely zero, the edge simply terminates. The value at any position is the product of the edge weights along the path that reaches it. An eight-by-eight matrix that would occupy 64 array entries can be stored in six nodes and 26 edges, and because quantum operations are highly regular this compression is dramatic in practice. Its side effect is that operations on the compressed form work directly on the shared structure, so repeated sub-computations are done once.

is the state of the art, and it was built entirely for processors that branch cheaply.

What people did before, and where it broke

Partition the state across threads. The standard parallelisation: split the amplitude vector, have each thread update its share. Where it broke. Nothing, but it is strong scaling on one input. Adding cores makes a single simulation faster and does nothing for the case where many simulations are wanted.
Fuse consecutive operations into one bigger matrix. Well established and effective, because it turns several passes over the state into one. Where it broke. Done on plain arrays, the fused matrix grows as two to the power of the number of quantum bits involved, so a ten-bit fused operation is a 1024 by 1024 array and further fusion becomes impractical.
Fuse them in the compressed representation instead. Recent work does exactly this, avoiding the stored zeros and the repeated blocks, which permits much more aggressive fusion. Where it broke. It is built for processors, and its cost model is about a single input. The measured gain from redoing that fusion for the batch setting turns out to be real but modest.
Use a graphics processor for the single-input case. Several mature simulators do this, and the vendor library even exposes an interface for applying one operation to a batch of states. Where it broke. That interface takes only dense matrices, so it cannot accept a heavily fused operation without exceeding memory — in the comparison here it fails to run on most circuits when given aggressively fused input.
Fork one process per input. The obvious way to get a batch done, and embarrassingly parallel. Where it broke. It treats inputs as unrelated when in fact every one of them is applying the same sequence of identical matrices, which is precisely the regularity a graphics processor is built to exploit.
The gap. The batch case is not a harder version of the single-input case; it is a differently shaped one, with abundant data parallelism and a fixed circuit that can be optimised once and reused. Nobody had built for that shape, and the compressed representation everyone relies on stands in the way, because a structure full of pointer chasing and branching is the opposite of what this hardware wants.

How a batch simulation gets built

From one input per process to a compiled task graph over many batches

↓ →
new
↓ →
new
↓ →
new
↓ →
new
↓ →
new

Sixteen circuits from a public benchmark suite up to 21 quantum bits, 200 batches of 256 inputs each, on a single 48 gigabyte graphics processor.

The idea

The batch case is a different problem, not a bigger one. Every existing simulator improves runtime by adding cores to a single input, which is the right target if you want one answer quickly. But testing, verification and state analysis want to know how a circuit behaves across many scenarios, so they push hundreds or thousands of different starting states through one fixed circuit. That workload has properties the single-input case does not: the arithmetic is identical across inputs, the circuit never changes so any preprocessing is paid once and amortised, and the whole thing becomes a large regular matrix product rather than a sequence of vector updates. Running one process per input satisfies the requirement while discarding all of that. Noticing that the parallelism a workload actually offers is not the parallelism the standard execution model reaches for is the same observation the hardware-accelerated circuit simulation work covered elsewhere on this site starts from, where a design is manifestly full of fine-grained parallelism that ordinary processors cannot use.

Use the compressed structure to reason about itself. The pivot of the fusion algorithm is a cost model: since the eventual computation is a sparse matrix product, an operation costs the number of multiply-accumulate steps needed per output amplitude, which is the largest number of non-zeros in any row of its matrix. That is the right quantity and it looks impossible to obtain, because the matrix has two to the power of n rows and nobody can afford to look at them. The resolution is that the row populations inherit the same regularity as the matrix — identical sub-blocks have identical population patterns — so the populations can themselves be stored as a decision diagram and computed recursively using the structure’s own concatenate and add operations, memoised per node. The compressed representation is used to answer a question about the thing it compresses, without ever expanding it. Refusing to materialise a regular structure and instead computing over its compact description is the same instinct as the bounded-repetition automata work covered elsewhere on this site, which stores a repetition count rather than unrolling it.

Once the arithmetic is cheap, launching the work costs more than doing it. The third stage addresses a bottleneck that only exists because the first two succeeded. A circuit of a thousand operations run over two hundred batches means hundreds of thousands of individual kernel launches, each cheap and collectively not, plus a transfer to and from the host per batch. So the kernels and the transfers are described as nodes in a dependency graph and handed to the hardware runtime to schedule as a unit, with four rotating buffers arranged so that one batch’s transfers overlap another’s computation. The ablation puts this at between 1.46 and 1.73 times, which is smaller than the other two stages and is a pure scheduling gain over identical arithmetic. Treating the overlap of communication and computation as a scheduling problem to be described declaratively rather than hand-orchestrated is the same approach the collective communication scheduling work covered elsewhere on this site takes for distributed training.

Following one gate matrix to the hardware

From a compressed graph to rows a thread block can process

Everything here happens once per circuit, then serves every batch.

Step 1 of 5

Work out what the operation will cost

Traverse the compressed graph depth-first, computing for each node the pattern of row populations of the sub-matrix it represents, combining children with the structure own concatenate and add operations, and caching the answer per node. The maximum over the result is the cost.

This is what makes the cost affordable to compute at all. Two sub-matrices that share a node share a population pattern, so the work is proportional to the size of the compressed graph rather than to the matrix.

Why it works

The right cost model is about memory traffic as much as arithmetic. In the flat sparse layout every output amplitude costs exactly the maximum row population in multiply-accumulate steps → so that single number prices an operation for the batch setting → and it explains why merging two cost-two operations into one cost-four operation is worth doing despite the arithmetic being identical, since it halves the passes over the state.
The compression that stores the matrix also prices it. Finding the largest row population means examining two to the power of n rows → but identical sub-blocks share a node and therefore share their population pattern → so the populations can be computed recursively over the compressed graph and memoised, costing time proportional to the graph rather than the matrix.
The chosen sparse format matches a measured property of these matrices. Padded formats waste work when row populations vary → but the measured coefficient of variation is exactly zero for three of the four circuits examined and 0.0328 for the fourth → so padding costs almost nothing here, and every thread gets equal work with predictable addresses.
Parallelising the conversion is not automatically right. Enumerating rows in parallel is the obvious use of the hardware → but each row is a graph traversal, and more edges means more branching and more divergence between threads in a group → so past a threshold the sequential version is faster, and choosing between them is worth doing.
Preprocessing is only free if the batch is large. Fusion and conversion happen once per circuit and cost real time → at ten batches they are more than half the runtime for a large circuit → at two hundred batches they are seven percent, which is why the whole design presumes the many-input workload it was built for.

What it costs (trade-offs and scope)

Two of the three headline comparisons measure a missing feature. Neither of the two largest speedups is against a system that supports batch simulation; both were run as eight parallel single-input processes because that is all they can do. The meaningful comparison is the 3.25 times over the one baseline that does support batching.
Against the closest prior work, the algorithmic contribution is modest. On the arithmetic actually eliminated, the new fusion beats the existing decision-diagram fusion by 1.23 times. It beats the array-based fusion by 3.85 times and the baseline with no fusion by 10.76 times, but the interesting comparison is the first one.
The biggest single contributor is a data structure change. The ablation attributes between 5.55 and 35.08 times to converting into the flat sparse layout, against 1.39 to 6.73 for the new fusion and 1.46 to 1.73 for the task graph. Most of the result is that the compressed graph is a bad fit for this hardware and a flat array is a good one.
One reported speedup is against a configuration that mostly does not run. The 407 times figure compares against the vendor library given aggressively fused operations, which its batch interface cannot accept because it takes only dense matrices — it fails on 11 of the 16 circuits. The comparable number is 3.62 times.
Some baseline runtimes look pathological. One simulator takes over six minutes on a 39-gate, six-qubit circuit, which produces a speedup figure above eleven thousand. That reflects fixed overhead in the baseline rather than simulation cost.
The largest circuit is 21 quantum bits. Scaling further is described as achievable by partitioning batches across several graphics processors, which is a plausible argument and is not demonstrated.
The choice of sparse format rests on uniformity that may not hold. Three circuits show perfectly uniform row populations and the fourth already does not. A circuit whose operations are irregular would pay for the padding.
The conversion threshold has no principled value. The authors state there is no universally optimal choice and parameterise it, setting it empirically for their machine.
Two baseline entries are timeouts. Runs exceeding 24 hours were terminated, so two of the reported speedups are lower bounds rather than measurements.
Preparation dominates at small batch counts, at 57 percent of runtime for the largest circuit at ten batches. The design is explicitly for workloads wanting hundreds.

The evaluation

What the three stages contribute, and against what

Property	One process per input	Vendor library with batch support
Overall speedup
Arithmetic eliminated by fusion		No fusion at all
Contribution of each stage, from the ablation	Not applicable	Not applicable
Behaviour with batch size	Linear in the number of processes	Baseline
Power		Baseline

Sixteen circuits from a public benchmark suite spanning machine learning, optimisation and supremacy workloads, validated by comparing output amplitudes against every baseline.Underlined values explain themselves when clicked.

Where this stops working

What to keep in mind

The two arresting numbers, 159 and 311 times, are largely measuring the absence of a feature: neither of those baselines supports batch simulation, so both were run as eight parallel single-input processes. The comparison that carries weight is 3.25 times over the one baseline that does batch. Within that, the new contribution is smaller than the headline suggests — measured on the arithmetic actually eliminated, the batch-aware fusion beats the closest prior fusion by 1.23 times, and the ablation attributes the largest share, between 5.55 and 35.08 times, to converting the compressed graph into a flat sparse layout, which is a data structure change rather than an algorithmic one. A separately reported 407 times figure compares against a configuration that fails on 11 of 16 circuits, because the vendor batch interface accepts only dense matrices and cannot hold aggressively fused operations; the runnable comparison there is 3.62 times. Some baseline runtimes are pathological enough to produce a speedup above eleven thousand on a six-qubit circuit, which reflects fixed overhead rather than simulation cost, and two entries are 24-hour timeouts, making those speedups lower bounds. The largest circuit tested has 21 quantum bits, with scaling beyond that argued rather than shown. The choice of sparse format depends on row populations being uniform, which holds perfectly for three of four circuits examined and already does not for the fourth. The conversion threshold is acknowledged to have no universally optimal value. And preparation costs 57 percent of runtime at ten batches, falling to seven percent only at two hundred.

Background

Why simulate a quantum computer on a classical one

It can seem perverse to spend so much effort simulating machines whose entire justification is doing things classical computers cannot. Three reasons make it necessary. The first is availability: quantum hardware is scarce, expensive and heavily queued, so an algorithm designer who wants to know whether their circuit is correct cannot simply run it. The second is that real devices are noisy, so a run on hardware conflates algorithmic errors with physical ones — a classical simulation gives the ideal answer, which is the only way to tell which is which. The third is that simulation gives access to quantities the hardware cannot provide at all: measuring a quantum state destroys it and yields one sample, whereas a simulator can hand over every amplitude after every operation. That last property is what makes simulation indispensable for testing and verification, and it is also why the exponential cost is unavoidable rather than an artefact of clumsy implementation. The interesting corollary is that classical simulability becomes a design resource in its own right — certain restricted families of circuits can be simulated efficiently, and the variational algorithm work covered elsewhere on this site exploits exactly that by confining a search to a simulable subspace to find good starting parameters before ever touching hardware.

If you remember only three things

1. Quantum circuit simulators are all built to push one input through one circuit as fast as possible, but testing, verification and state analysis need hundreds of inputs through a fixed circuit — a differently shaped problem where the arithmetic is identical across inputs, preprocessing is paid once, and the whole thing becomes a large regular matrix product that suits a graphics processor, which running one process per input throws away entirely.
2. The fusion algorithm rests on pricing each operation by the largest number of non-zeros in any row of its matrix, since that is the multiply-accumulate count per output amplitude — a quantity that appears to require examining two to the power of n rows, but which inherits the same regularity as the matrix, so it can be computed recursively over the compressed representation itself and memoised, never expanding anything.
3. The honest accounting is that the striking speedups of 159 and 311 times are against systems that cannot do batch simulation at all, the meaningful figure is 3.25 times against the one that can, the batch-aware fusion beats the closest prior fusion by 1.23 times on arithmetic eliminated, and the largest single contributor in the ablation is converting the compressed graph into a flat sparse layout that the hardware handles well.

Check yourself

Check your understanding

1. Why is running one process per input a poor way to simulate a batch?
2. How is the cost of an operation computed without expanding its matrix?
3. Why is merging two operations that each cost two into one that costs four still worthwhile?
4. Which comparison in the evaluation carries the most weight, and why?

Sources

The paper is “BQSim: GPU-accelerated Batch Quantum Circuit Simulation using Decision Diagram” by Shui Jiang, Yi-Hua Chung, Chih-Chun Chang, Tsung-Yi Ho and Tsung-Wei Huang, published at ASPLOS 2025, pages 79-94, DOI 10.1145/3676641.3715984. It is open access under a Creative Commons licence, the source code is public, and the paper includes an artifact appendix.

This explainer is written from the full 16-page text. Everything above comes from the paper’s own sections: the background on full-state simulation and how single-qubit and controlled operations update amplitudes, the worked construction of decision diagrams for both a gate matrix and a state vector including how values are recovered as products of edge weights, and the account of why array-based fusion becomes impractical as fused matrices grow; the algorithm covering the definition of batch cost as maximum row population, the recursive computation of the row-population vector over the diagram using its native concatenate and add operations with memoisation, and the three-step fusion procedure with its worked example; the conversion stage covering the measured coefficient of variation of row populations across four circuits, the reason for choosing a padded format over alternatives, the representation of the diagram on the device as node and edge arrays, the conversion kernel with its per-block depth-first traversal using an explicit stack in shared memory, and the finding that the parallel conversion loses to the sequential one as edge counts grow together with the resulting hybrid selection; the task graph stage covering kernel launch and transfer overheads, the four-buffer rotation with its indexing scheme, and the worked dependency example across four batches; and the evaluation covering the machine and benchmark configuration, overall runtimes against three baselines with their differing batch capabilities, the comparison of arithmetic eliminated, conversion timing, task graph results including the configuration that fails on most circuits, batch size scaling to saturation, power measurements, the runtime breakdown showing how preparation amortises, and the three-way ablation.