New StateVector simulator backend based on cuQuantum for Graphix

[ljcamargo]

Closes TeamGraphix/graphix#496 for UnitaryHack 2026

Hello @matulni here is the PR for the cuQuantum StateVector backend.
Some comments:

• Naming convention for this new backend, currently still uses template offered name
• Tested on T4 GPU.
• Log file included which shows sanity checks, tests and benchmarks passed.
• Mind I didn't remove the helper bash script used to speed up the tests on GPU server (run_tests.sh)

Copy of the backend README:

graphix-statevec-cuquantum

GPU-accelerated statevector backend for Graphix pattern simulation, built on NVIDIA cuQuantum (cuStateVec) and CuPy.

Requirements

• NVIDIA GPU with CUDA 12.x
• Python ≥ 3.13
• cuquantum-python-cu12 (≥ 24.03) — pre-built wheels, no CUDA toolkit needed at build time
• cupy-cuda12x (≥ 13.0)

Installation

uv lock
uv sync --extra cuquantum --dev

Usage

from graphix_statevec_cuquantum import Statevec, StatevectorBackend
from graphix.transpiler import Circuit

qc = Circuit(2)
qc.cz(0, 1)
qc.h(0)
pattern = qc.transpile().pattern

backend = StatevectorBackend()
sv = pattern.simulate_pattern(backend=backend)
print(sv.flatten())

Design

State representation

The quantum state is stored as a flat CuPy array of shape (2**max_space,) with dtype complex128. Only the first 2**nqubit entries are meaningful; the tail is padding to avoid GPU reallocation as qubits come and go during MBQC pattern simulation.

Qubit ordering

Graphix uses MSB convention (qubit 0 = most significant bit = axis 0 of the tensor). cuQuantum uses LSB convention (qubit 0 = least significant bit). All calls to cuQuantum APIs convert indices via _msb_to_lsb().

cuQuantum API used

Operation	cuQuantum function	Notes
Single/multi-qubit gate	custatevec.apply_matrix	Workspace size queried via apply_matrix_get_workspace_size
Expectation value	custatevec.compute_expectation	Result written to host-side buffer, returned as complex
Qubit swap	cp.swapaxes (CuPy)	cuQuantum swap_index_bits available but CuPy tensor axis swap is simpler and equally GPU-accelerated
Entangle (CZ)	custatevec.apply_matrix with 4×4 CZ matrix

Buffer growth

When tensor() needs more space than max_space, the buffer is reallocated to exactly the needed size. No doubling strategy is used, since MBQC pattern qubit counts fluctuate and doubling can overshoot GPU memory.

Project structure

src/graphix_statevec_cuquantum/
    __init__.py                        # Exports: Statevec, StatevectorBackend
    graphix_statevec_cuquantum.py       # Main implementation
    README.md                           # This file

tests/
    test_statevec_cuquantum.py          # Unit tests + legacy comparison + pattern sim

benchmarks/
    bench_statevec_cuquantum.py         # Performance benchmarks

run_tests.sh                            # Full test suite script

CI notes

Tests and benchmarks are decorated with @pytest.mark.skipif(not _gpu_available(), reason="GPU not available"). On CI runners without GPU they will be skipped. Formatting, linting, and type-checking still run normally.

Benchmark results (Tesla T4)

Approximate per-operation latency (lower is better):

Operation	4 qubits	8 qubits	12 qubits	16 qubits
evolve_single (H gate)	45 µs	46 µs	45 µs	46 µs
entangle (CZ)	52 µs	52 µs	52 µs	52 µs
expectation_single (Z)	85 µs	80 µs	81 µs	106 µs
add_nodes	354 µs	423 µs	569 µs	1.6 ms
remove_qubit	389 µs	477 µs	614 µs	1.6 ms

Gate application times are near-constant across qubit counts due to cuQuantum's optimized GPU kernels. Structural operations (add_nodes, remove_qubit) scale with state size as they involve GPU memory transfers.

Caveats

• cuQuantum API version: This implementation targets cuQuantum Python 26.x (cuquantum.bindings.custatevec). Earlier versions (24.x) with cuquantum.custatevec and apply_gate are not compatible.
• compute_expectation: cuQuantum's own API writes the scalar result to a host-side (CPU) buffer. The computation is fully GPU-accelerated; only the final 16-byte complex value is transferred to CPU.
• Buffer reallocation: Each tensor() growth triggers a new GPU allocation. For patterns with predictable max_space, initializing Statevec with the correct max_space avoids this.
• No noise model: The backend does not implement apply_noise.
• Naming: Statevec and StatevectorBackend follow the template's convention. The import path (graphix_statevec_cuquantum) disambiguates from the NumPy CPU backend (graphix.sim.statevec) and the template (graphix_statevec_template). Subject to change based on Graphix maintainer guidelines.
  test_output_20260606_211025.log

[matulni]

Thanks, I'll look into it in the coming days. Don't worry about the typecheck failure in the CI, it's a known issue with mypy and qiskit. I thought I had fixed the dependency, but maybe qiskit 2.4 version is leaking from the mqt-bench plugin, I'll try to fix it.

[matulni]

Looking nice! I managed to run your bash script and tests are passing. I'll review your code shortly.

[matulni]

Thanks for your work, it looks very nice!

Here are some comments from my first pass

[ljcamargo]

hi @matulni I'm working on all your suggestions, I'll get back to you as soon as it is all done.

[ljcamargo]

hi @matulni except for the comment I left above the tests and lints are now passing for all the repo and the code suggestions applied, there is another new helper script run_validations that runs the previous mentioned commands but for all repo which ease the testing in local environments. Nonetheless, I will dig deeper in the Tuple v. tuple issue later today.

[matulni]

Hi @ljcamargo, we decided to award you the bounty for your contribution, congratulations!
In the remaining time before the hackathon ends, could you please address the remaining comments from my previous reviews ? I see many comments still unresolved. Thank you.
A second maintainer @thierry-martinez may briefly review your PR too. If that were the case, please address his comments too.

[ljcamargo]

Hi @matulni and @thierry-martinez thanks a lot for choosing me! I'm working on solving all the comments, they remain a few, nonetheless in the last batch of suggestion of yours I've applied, they raised many new lint and runtime errors, as they were many this time, I used ai analysis this time to asses which suggestion lead to which errors and got the following, I suggest reverting this to keep the tests running and tackle other fixes while you check this comments:

1. IndexError — remove_qubit crashes on all runs

  - branch = t[tuple(idx)].ravel()
  + branch = t[idx].ravel()

Cupy does not accept a plain list as an advanced index — it needs tuple(...). This caused every remove_qubit call
to raise IndexError, which cascaded into all 12 test failures (GPU remove, CPU remove pattern, and pattern
simulator).

2. ValueError: Unsupported dtype object — secondary crash in remove_qubit

  - branch_nrm2 = float(cp.sum(cp.abs(branch) ** 2))
  + branch_nrm2 = cp.sum(cp.abs(branch) ** 2)

cp.sum() returns a cupy scalar, not a Python float. math.isclose() downstream expects a Python float. The float()
cast is required.

3. mypy type error — with_capacity

  - backend = cls(**kwargs)
  - object.__setattr__(backend, "state", state_init)
  - return backend
  + return cls(state_init, **kwargs)

state has init=False on the frozen dataclass, so calling cls(state_init, ...) passes state_init (a Statevec) to the
first positional parameter node_index: NodeIndex. This is a type mismatch. The object.setattr pattern was the
correct fix for frozen dataclass fields with init=False.

(THIS ONE IS A CONCERN, WHICH AI ADVERTS OF POSSIBLE UNWANTED SIDE FX)
4. Semantic concern — fidelity now returns complex instead of squared modulus

  - ip = float(cp.dot(a.conj(), b))
  - return ip.real**2 + ip.imag**2
  + ip: float = cp.dot(a.conj(), b)    # cp.dot returns complex, not float
  + return ip                           # now returns complex, not |<ψ1|ψ2>|²

cp.dot() on complex arrays returns a complex value. The annotation ip: float is misleading — the value is complex.
More importantly, fidelity() was documented to return |⟨ψ1|ψ2⟩|² (a float in [0,1]), but now returns the raw
complex inner product. Calling code like isclose(), which does math.isclose(self.fidelity(other), 1, ...), will
break for states with non-real overlap (the result won't be comparable to 1). Needs float(...) and the
squared-modulus computation to match the docstring.

[thierry-martinez]

You've my approval! Thank for your contribution.

[matulni]

Hi @ljcamargo, thanks for your contribution and apologies for the misleading comments. We'll merge the PR!