GBSKernels

A GPU-native, batched library of the #P-hard matrix functions behind photonic quantum sampling — the permanent, hafnian, loop hafnian, and torontonian — with an explicit floating-point accuracy model that ranges from native double precision to rigorous a-posteriori error bounds.

Gaussian boson sampling and its displaced- and threshold-detector variants, and standard boson sampling, reduce computationally to evaluating large batches of independent, medium-sized instances of these four functions. Each is an alternating signed sum whose value can be many orders of magnitude smaller than its terms, so it is cancellation-prone in floating point: both throughput and known accuracy matter. GBSKernels targets exactly this workload — batched-first kernels on CPU and CUDA, every function validated against independent combinatorial ground truth, and a precision model in which the accuracy of every result is either measured or provably bounded.

GBSKernels is designed to complement The Walrus, the canonical CPU/C++ reference implementation, rather than to replace it: a GPU-native, batched companion with an explicit accuracy characterization. The Walrus is used here as one of several independent test oracles.

Status. The CPU library is complete and pip-installable (pure Python, numpy + mpmath only). Historical release validation covers NVIDIA hardware (RTX 4090, A100); the v0.2 changes pass the CPU host-shim gates and still require a fresh on-device validation session. This is pre-1.0 software and the public API may change.

Version 0.2.0 adds the fail-closed Jiuzhang confirmatory workflow and corrects the real threshold-torontonian construction used by the Jiuzhang examples. See the changelog for the release scope and compatibility notes.

Installation

Not yet on PyPI; install from source. The CPU backend needs only numpy + mpmath; the CUDA extension is a separate, optional build (see bindings/README.md):

uv build                        # -> dist/gbskernels-*.whl, then pip install it
# pip install gbskernels        # planned once published to PyPI

Quick start

import numpy as np
import gbskernels

A = np.array([[1., 2.], [3., 4.]])
gbskernels.perm(A)                              # (10+0j)

# The API is batched-first: one call evaluates a stack of matrices.
rng = np.random.default_rng(0)
stack = rng.standard_normal((4096, 8, 8))
stack = stack + stack.transpose(0, 2, 1)        # symmetric, as the hafnian needs

gbskernels.haf_batched(stack)                   # double precision (default)
gbskernels.haf_batched(stack, precision="auto") # double precision + cancellation guard
gbskernels.haf_batched(stack, backend="gpu")    # CUDA, if the extension is built
gbskernels.haf(stack[0], precision="ref")       # arbitrary-precision reference

A short guided tour — closed-form checks, the auto tier catching real double-precision cancellation, a GBS distribution, and a sampling run — runs with wheel dependencies only:

python examples/gbs_demo.py

Precision model

Precision is always an explicit tier, so the accuracy of a result is never left implicit:

Features

• All four functions, batched, on CPU and CUDA. Field-standard algorithms (Glynn/BB–FG for the permanent, the power-trace form for the hafnian and loop hafnian, subset-determinant inclusion–exclusion for the torontonian), each with an independent naive reference used for verification.
• A measured accuracy boundary. The double-precision–vs–double-double crossover is characterized for every function against an arbitrary-precision reference, on physical and adversarial inputs.
• Certified evaluation. Every value can be returned with a rigorous error bound whose enclosure of the true value is a hard test invariant — verified on closed forms, random ensembles, and adversarial cancellation families.
• A device-resident workspace. gbskernels.Workspace evaluates ragged batches with a persistent stream, pinned staging, and a zero-copy DLPack output that CuPy / PyTorch / JAX can consume without a device-to-host copy (docs/device_resident_contract.md).
• A conditional GBS sampler (sampling/) that draws photon-number samples by the chain rule, evaluating each mode's batch of hafnians on the GPU, validated distributionally (total-variation distance below 0.03 against both the exact distribution and The Walrus, and a chi-square goodness-of-fit p > 10⁻³ against The Walrus).
• Structure-aware kernels for the sampling workload: a repeated-row finite-difference sieve for the loop hafnian, and a recursive prefix-Cholesky torontonian — including a single-large mode that splits one evaluation across the GPU to dimension 64 (32 modes).
• A fail-closed confirmatory workflow for the public Jiuzhang example, covering exposure audit, frozen design, future-beacon selection, content-addressed evaluation, refusal recovery, simultaneous coherence-grid inference, predictive checks, and release verification (protocol, tools). The repository ships templates and validation contracts, not a fabricated registration or completed v2 outcome.

Verification

Verification is treated as the deliverable (docs/DESIGN.md, §8). Five layers are exercised by the test suite on every commit:

1. Independent combinatorial ground truth — perfect-matching counts and closed forms (perm(J_n) = n!, haf(J_{2n}) = (2n−1)!!, …), sharing no code with any existing library.
2. Differential testing against The Walrus across sizes, seeds, and real/complex inputs.
3. Property-based invariants (Hypothesis): scaling, permutation, transpose, the loop-hafnian → hafnian reduction, and others.
4. Statistical / end-to-end validation — the kernels compute GBS distributions that sum to one and match The Walrus per pattern in all three detector regimes.
5. Numerical-accuracy characterization against the arbitrary-precision reference, including adversarial cancellation families that defeat double precision and must survive double-double.

The Jiuzhang v2 tooling adds workflow-contract tests for canonical hashing, registration timing, exclusion-ledger completeness, selection, immutable run reduction, refusal recovery, inference, reconstruction, and release assembly.

The development discipline is CPU-first: host-shim-compatible CUDA kernels are compiled and run through the CPU pre-flight (part of the test suite), then the full source set is checked on-device against independent CPU references before any throughput result is considered publishable. Device-only variants such as the warp-specialized permanent are covered by the on-device gates; profiler diagnostics collected before those gates remain provisional until they pass.

Performance

The library is batched-first: throughput comes from evaluating thousands of independent instances per launch, and the benchmark harness reports where the CPU wins (small single evaluations) as plainly as where the GPU wins (large batches), always alongside the accuracy each point was measured at. The structure-aware kernels give substantial speedups in their regimes — the recursive torontonian over the LU form, and the finite-difference sieve over the expanded power-trace at repeated-row patterns. Raw benchmark data is written, append-only, to results/; the protocol is described in docs/benchmark_protocol.md.

python -m bench.accuracy                          # the accuracy boundary, all four functions
python -m bench.throughput --func perm --sizes 4,6,8,10
bash scripts/gpu_session.sh                       # full on-device session (on a CUDA host)

Scope and limitations

These are correctness-first kernels with a deliberately narrow scope:

• One evaluation per GPU thread is the default mapping. Warp/block-cooperative variants exist and are dispatched only where they measurably win (the cooperative permanent); for the hafnian, loop hafnian, and torontonian the per-thread linear algebra is memory-bound, and the cooperative strategy does not help them.
• The four functions do not share a single subset-enumeration engine. Only the permanent uses the Gray-code delta walk; the others enumerate masks independently, because their per-subset work is memory-bound.
• Fixed size limits (per-thread buffers): permanent ≤ 28, hafnian ≤ 20 (DD ≤ 16), loop hafnian ≤ 20 (DD ≤ 14), torontonian 2n ≤ 24. Larger inputs raise on backend="gpu"; use the CPU backend. The single-large recursive torontonian extends to dimension 64.
• The double-double tier is internal: kernels carry double-double through the cancelling summation and collapse to complex128 on output — they recover a correct double-precision answer under cancellation rather than exposing extended-precision results. The double-double torontonian is real-input only.
• The sampler is pure-state and hybrid host-orchestrated (the host drives the chain rule; the GPU evaluates each step's batch of hafnians). A fully on-device variant is available but limited to double precision within a submatrix-size cap.
• Single GPU, CUDA only — no multi-GPU, ROCm, or automatic differentiation.

Repository layout

gbskernels/     batched-first Python API (the public surface)
cpu_ref/        independent CPU reference implementations (double + double-double)
highprec_ref/   arbitrary-precision reference (mpmath)
core/           CUDA C++: subset-enumeration utilities, the kernels, on-device gates
bindings/       nanobind extension exposing the GPU backend to Python
sampling/       boson-sampling orchestration and the conditional GBS sampler
examples/       runnable CPU demo plus Jiuzhang reproduction and v2 workflow tools
tests/          the five-layer verification suite
bench/          accuracy and throughput harnesses
scripts/ envs/  scripted GPU-session runner and container definitions
docs/           design, benchmark/device contracts, and v2 protocol/templates
results/        append-only validation, gate, benchmark, and profiler artifacts

Development

uv sync                      # Python 3.12 environment, dependencies, editable install
uv run pytest                # the CPU verification suite
uv run pytest -m layer1      # the independent combinatorial ground truth
uv run pytest -m "not slow"  # the fast tier

Building the CUDA extension requires nvcc and is a separate, optional step (bindings/README.md); the pure-Python wheel never requires it. The full design rationale — algorithms, precision strategy, and the verification contract — is in docs/DESIGN.md.

Relationship to The Walrus

GBSKernels would not exist without The Walrus, which remains the canonical reference implementation and this library's primary differential oracle. The scope here is deliberately narrower: the four batched kernels and their precision characterization.

License

Apache-2.0.