CUDA Quantum Simulator

Status: not actively maintained. Left up as a reference.

A quantum state-vector simulator in CUDA C++, with noise models, density-matrix simulation, and benchmarks against both a CPU reference and NVIDIA cuStateVec.

What's included

• State-vector simulation on NVIDIA GPUs, double precision
• Gates: X, Y, Z, H, S, T, S†, T†, Rx, Ry, Rz, CNOT, CZ, SWAP, Toffoli, CRY, CRZ
• Noise channels: depolarizing, amplitude damping (T1), phase damping (T2), bit/phase flip
• Density-matrix simulation via Kraus operators (1-14 qubits)
• Batched trajectories for Monte Carlo sampling
• Single-qubit measurement with collapse, plus multi-shot sampling
• A chainable circuit API and RAII GPU memory
• GPU output checked against a CPU reference and against Qiskit/Cirq; gate conventions follow Qiskit (little-endian)

Benchmark Results

Tested on NVIDIA RTX 4070 Laptop GPU (8GB VRAM, Compute Capability 8.9)

GPU vs CPU Performance

The GPU simulator passes single-threaded CPU at 12 qubits, and the gap widens as qubit count grows (100 mixed H + CNOT gates):

Qubits	GPU (ms)	CPU (ms)	Speedup
10	0.33	0.14	0.4x (overhead dominates)
12	0.29	0.54	1.9x
14	0.29	2.14	7.5x
16	0.29	8.62	29.8x
18	0.29	35.46	123x
20	0.28	143.50	515x
22	0.28	720.39	2,551x

GPU time stays near-constant (~0.28 ms) because 100 gates saturate the GPU at any qubit count in this range; CPU time grows exponentially (O(2^n) per gate).

Scaling Characteristics

Qubits	States	Memory	Init (ms)	100 H gates (ms)
20	1M	16 MB	0.12	0.25
22	4M	64 MB	0.56	0.25
24	16M	256 MB	2.28	0.24
26	67M	1 GB	7.07	0.24

Gate Throughput (20 qubits)

Gate	Throughput
CNOT	53,200 gates/s
X	34,300 gates/s
H	24,600 gates/s
Rz	5,940 gates/s

Requirements

• CUDA Toolkit 12.0 or later
• CMake 3.18 or later
• C++17 compatible compiler
• NVIDIA GPU with Compute Capability 7.0+ (tested on RTX 4070)

Building

# Clone and build
cd cuda-quantum-simulator
mkdir build && cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
make -j$(nproc)

# Run tests
ctest --output-on-failure

# Run demo
./quantum_sim

# Run benchmarks
./benchmark_gates
./benchmark_scaling

Usage

Basic Circuit Execution

#include "Simulator.hpp"
using namespace qsim;

// Create a 2-qubit Bell state
Simulator sim(2);
Circuit circuit(2);
circuit.h(0).cnot(0, 1);  // |Φ+⟩ = (|00⟩ + |11⟩)/√2

sim.run(circuit);
auto probs = sim.getProbabilities();
// probs[0] ≈ 0.5, probs[3] ≈ 0.5

Measurement and Sampling

// Sample from the probability distribution (non-destructive)
auto samples = sim.sample(1000);

// Measure a single qubit (collapses state)
int result = sim.measureQubit(0);

Factory Circuits

// Pre-built circuits for common states
auto bell = Circuit::createBellCircuit();
auto ghz = Circuit::createGHZCircuit(5);
auto random = Circuit::createRandomCircuit(4, 20);  // 4 qubits, depth 20

Noisy Simulation

#include "NoiseModel.cuh"
using namespace qsim;

// Create a noise model with depolarizing error
NoiseModel noise;
noise.addDepolarizing(0.01);  // 1% depolarizing probability per gate

// Run noisy simulation (Monte Carlo wavefunction method)
NoisySimulator sim(3, noise, /*seed=*/42);
Circuit circuit(3);
circuit.h(0).cnot(0, 1).cnot(1, 2);

sim.run(circuit);
auto probs = sim.getProbabilities();

Batched Monte Carlo Sampling

// Run 1000 noisy trajectories in parallel on GPU
NoiseModel noise;
noise.addDepolarizing(0.005);
noise.addAmplitudeDamping(0.001);  // T1 decay

BatchedSimulator batch(/*n_qubits=*/3, /*batch_size=*/1000, noise);
Circuit circuit(3);
circuit.h(0).cnot(0, 1).cnot(1, 2);

batch.run(circuit);

// Get average probabilities across all trajectories
auto avg_probs = batch.getAverageProbabilities();

// Get histogram of measurement outcomes
auto histogram = batch.getHistogram();  // Counts for each basis state

Available Noise Channels

Noise Type	Description	Usage
Depolarizing	Random Pauli error (X, Y, or Z with equal probability)	noise.addDepolarizing(p)
Amplitude Damping	T1 decay - relaxation to ground state	noise.addAmplitudeDamping(gamma)
Phase Damping	T2 dephasing - loss of phase coherence	noise.addPhaseDamping(gamma)
Bit Flip	X error with probability p	noise.addBitFlip(p)
Phase Flip	Z error with probability p	noise.addPhaseFlip(p)
Bit-Phase Flip	Y error with probability p	noise.addBitPhaseFlip(p)

Density Matrix Simulation

For exact simulation of mixed states and noise (without Monte Carlo sampling), use the density matrix simulator:

#include "DensityMatrix.cuh"
using namespace qsim;

// Create a 3-qubit density matrix simulator
DensityMatrixSimulator sim(3);
Circuit circuit(3);
circuit.h(0).cnot(0, 1).cnot(1, 2);

sim.run(circuit);

// Apply exact noise channels (Kraus operators, not Monte Carlo)
NoiseModel noise;
noise.addDepolarizing(0.01);
sim.applyNoise(noise);

// Get probabilities and state properties
auto probs = sim.getProbabilities();
double purity = sim.getPurity();  // tr(rho^2), 1.0 for pure states

When to use density matrix vs state vector:

Approach	Memory	Use Case
StateVector	O(2^n)	Pure states, large qubit counts (20+)
NoisySimulator	O(2^n)	Noisy circuits, Monte Carlo sampling
DensityMatrix	O(4^n)	Exact mixed states, small systems (1-14 qubits)

Architecture

Core Components

include/
├── StateVector.cuh    # GPU-resident quantum state (2^n complex amplitudes)
├── Gates.cuh          # CUDA kernels for quantum gates
├── Circuit.hpp        # Circuit representation with fluent API
├── Simulator.hpp      # GPU simulator orchestration
├── NoiseModel.cuh     # Noise models and batched simulation
├── DensityMatrix.cuh  # Density matrix simulation for mixed states
├── OptimizedGates.cuh # Shared memory and coalesced access kernels
├── CudaMemory.cuh     # RAII wrapper for CUDA memory
└── Constants.hpp      # Configuration and math constants

src/
├── StateVector.cu     # State vector implementation
├── Gates.cu           # Gate kernel implementations
├── Circuit.cpp        # Circuit builder
├── Simulator.cu       # GPU and CPU simulator implementations
├── NoiseModel.cu      # Noise model and batched simulator
├── DensityMatrix.cu   # Density matrix implementation
├── OptimizedGates.cu  # Optimized kernel implementations
└── main.cpp           # Demo executable

Qubit Ordering Convention

This simulator uses little-endian ordering, matching Qiskit's default: qubit q is bit q of the state index.

• Qubit 0 is the least significant bit; qubit n−1 is the most significant
• For a 3-qubit state: index = q0×1 + q1×2 + q2×4, so X on qubit 0 takes |000⟩ (index 0) to index 1

Example for 3 qubits (kets written most-significant-qubit-first, |q2 q1 q0⟩):

Index:  0     1     2     3     4     5     6     7
State: |000⟩ |001⟩ |010⟩ |011⟩ |100⟩ |101⟩ |110⟩ |111⟩

Validation

The simulator is validated through multiple approaches:

1. GPU vs CPU Equivalence (Primary Validation)

The C++ test suite verifies GPU kernels produce identical results to a CPU reference implementation:

cd build && ./test_gpu_cpu_equivalence  # All tests pass

2. Gate Algebra Tests

Mathematical identities verify correctness (H²=I, CNOT²=I, S²=Z, etc.):

cd build && ./test_gate_algebra  # All tests pass

3. Cross-Validation Against Qiskit and Cirq

These scripts run each circuit on both this simulator (via the dump_state tool) and the reference library, then compare full state vectors up to a global phase:

cmake --build build --target dump_state          # one-time build of the tool
python validation/validate_against_qiskit.py     # Qiskit (little-endian, same as us)
python validation/validate_against_cirq.py        # Cirq (big-endian; qubit order reversed)

Cross-validation coverage:

• All single-qubit gates (X, Y, Z, H, S, T, S†, T†, Rx, Ry, Rz)
• All two-qubit gates (CNOT, CZ, SWAP, CRY, CRZ) and Toffoli
• Bell states, GHZ states, and random circuits

The C++ test_gpu_cpu_equivalence suite additionally checks random circuits up to 500 gates deep against the CPU reference.

Testing

The test suite uses Google Test and covers:

Test Suite	Tests	Description
test_warmup	4	CUDA infrastructure: vector add, shared memory, GPU properties, bandwidth
test_statevector	16	State initialization, normalization, measurement, move semantics
test_gates	26	Gate correctness for all 17 gate types
test_gate_algebra	24	Gate identities (H²=I, X²=I, S²=Z, T⁸=I, CNOT²=I, etc.)
test_gpu_cpu_equivalence	13	GPU matches CPU reference implementation
test_boundary	18	Edge cases and error handling
test_noise	25	Noise models, NoisySimulator, BatchedSimulator
test_density_matrix	28	Density matrix operations and Kraus noise channels
test_optimized_gates	8	Optimized kernels match standard kernels (incl. shared-memory path)

Total: 162 test cases across 9 test suites — all passing.

# Run all tests
cd build && ctest --output-on-failure

# Run specific test suite
./test_statevector --gtest_filter='*Measurement*'

Memory Safety (Valgrind)

All 9 test suites pass with zero memory leaks verified by Valgrind 3.22:

definitely lost: 0 bytes in 0 blocks   ← no application leaks
indirectly lost: 0 bytes in 0 blocks

CUDA programs produce known false positives (unhandled ioctl 0x30000001 from the NVIDIA kernel driver, plus context-lifetime "still reachable" memory). These are suppressed via cuda.supp. The valgrind.sh script runs all suites and exits non-zero if any definite or indirect leak is detected:

./valgrind.sh           # build + run all 9 suites under valgrind
./valgrind.sh --no-build  # skip build step

RAII ensures every GPU allocation is freed automatically:

Class	Resource	Freed by
CudaMemory<T>	cudaMalloc allocation	~CudaMemory() / move assignment
StateVector	device state array	~StateVector()
DensityMatrix	device ρ + scratch	~DensityMatrix()
NoisySimulator	state + RNG states	CudaMemory members
BatchedSimulator	batched states + RNG	CudaMemory members

Technical Details

Memory Requirements

State vector size = 2^n × 16 bytes (double precision complex):

• 20 qubits: 16 MB
• 25 qubits: 512 MB
• 28 qubits: 4 GB (practical limit for 8GB GPU)

Performance Characteristics

• CUDA initialization overhead: First kernel launch includes JIT compilation and context setup
• Memory bandwidth bound: Performance scales with state vector size at high qubit counts
• GPU advantage: Emerges at 20+ qubits where parallelism outweighs overhead

Design Decisions

1. Double precision only: Ensures numerical stability for deep circuits
2. Synchronous execution: Simpler debugging, consistent timing
3. Little-endian ordering: Matches Qiskit; Cirq comparison requires reversing the qubit order
4. RAII everywhere: No manual cudaFree calls, exception-safe

Benchmarks against cuStateVec

Multi-GPU support is the main thing not implemented.

cuStateVec Comparison Results

Benchmarked against NVIDIA cuStateVec 1.11.0 (part of cuQuantum SDK):

Gate	Qubits	Our Time	cuStateVec Time	Speedup
Hadamard	20	0.035 ms	0.067 ms	1.9x faster
Hadamard	24	2.7 ms	2.7 ms	1.0x (equal)
Hadamard	26	9.9 ms	9.7 ms	1.0x (equal)
CNOT (adj pair, 20q)	20	0.004-0.025 ms	0.033 ms	1.3-8.9x faster

Circuit Benchmark (H + CNOT layers, depth 10):

Qubits	Our Throughput	cuStateVec Throughput	Speedup
20	48,791 gates/s	20,238 gates/s	2.4x faster
24	637 gates/s	547 gates/s	1.2x faster

Key Findings:

• Our simple kernels match or exceed cuStateVec performance
• At smaller qubit counts (12-20), we're significantly faster due to lower overhead
• At larger qubit counts (24-26), performance converges as both become memory-bound
• Our CNOT implementation is particularly efficient for adjacent qubit pairs

Optimization Benchmark Results

The optimized kernels were benchmarked against the original implementations:

Optimization	Speedup	Notes
Shared memory tiling (target qubit 0)	1.5x	Best case for low-order qubits
Shared memory tiling (target qubits 1-7)	~1.0x	L2 cache already effective
Coalesced access patterns	1.0-1.1x	Marginal improvement

Finding: Modern GPU memory hierarchies (RTX 4070 L2 cache, high bandwidth) already handle random access patterns efficiently. The original naive kernels were near-optimal for this hardware, demonstrating that simple implementations can be surprisingly performant on modern GPUs.

License

MIT (see LICENSE).