Fix vector/matrix product (@) silently overflowing for narrow value dtypes - #44
Merged
Merged
Conversation
…types _major_matvec/_minor_matvec/_major_matmat/_minor_matmat allocated their output accumulator as values.dtype, the compact per-element storage dtype (e.g. uint16 for single-cell counts), instead of a dtype sized for the accumulated total. Contributions landing in the same output slot wrapped modulo 2**bits instead of promoting, unlike sum() which already accumulates in float64. Promote values and the other operand to their common numpy dtype (np.result_type) before running the kernels, matching the dtype promotion an equivalent dense-array product would get. Fixes #43. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
fishidaho
added a commit
that referenced
this pull request
Sep 18, 2026
`_ops._major_matvec`/`_major_matmat` -- the scatter direction, reached by
`VCSC @ B` and `B @ VCSR` -- ran single-threaded with no strategy at all,
losing to scipy by up to 4x:
VCSC @ x 8.7 ms vs scipy 6.7 ms 1.30x
VCSC @ B (k=8) 20.1 ms vs scipy 21.0 ms 0.96x
x @ VCSR 20.5 ms vs scipy 5.1 ms 3.99x
B @ VCSR (p=8) 28.5 ms vs scipy 14.6 ms 1.96x
They cannot be parallelized over the major axis directly, since two major
slices can collide on the same output index.
Regrouping a chunk into the opposite format and running an aligned kernel --
the strategy `_vcs_matmul` uses for the normalized view -- was tried first and
is much worse here. That path caches its transposed copy and amortizes it
across an iterative algorithm's many products; a bare `A @ B` pays the regroup
once per call, which turned the 8.7 ms product into 700 ms. Reverted.
Thread-local accumulators instead -- the shape `minor_sums`/`minor_extrema`
already use -- reduced across threads afterwards. The thread count is the
whole design problem: the scatter is `nnz * width` work split across threads
while the reduction is `nthreads * n_minor * width`, so more threads is not
better. Measured on 6M nonzeros over 60k x 2k, letting the existing 64 MiB
accumulator budget pick 48 threads for a width-8 product ran it in 56.6 ms,
*slower than the serial kernel's* 23.3 ms; 4 threads and 15 MB ran it in
11.6 ms.
`scatter_threads()` takes the smaller of two bounds. Setting the derivative of
scatter-plus-reduction to zero puts the optimum at `sqrt(nnz / n_minor)`,
which is 10 for that array against a measured best of 16/8/4 at widths 1/4/8 --
right in magnitude, but blind to width. A 16 MiB byte cap supplies the missing
width-dependence, admitting fewer threads exactly as the accumulator grows.
Together they track the measured optimum across widths:
VCSC @ x 2.4 ms vs scipy 6.7 ms 0.35x (was 1.30x)
VCSC @ B (k=8) 12.0 ms vs scipy 21.6 ms 0.56x (was 0.96x)
x @ VCSR 1.8 ms vs scipy 5.2 ms 0.36x (was 3.99x)
B @ VCSR (p=8) 5.2 ms vs scipy 14.3 ms 0.36x (was 1.96x)
Peak allocation stays bounded by the cap rather than by `nnz`: 5.3 MB at
width 1, 19.2 MB at width 8, 23.0 MB at width 16. Single-thread cases fall
back to the serial kernels rather than pay an allocation and a reduction pass
for one partial.
Three benchmark cases now cover this direction; the existing
`matvec_vs_scipy`/`matmat_vs_scipy` only ever ran VCSR aligned.
Rebased from the pre-0.4.0 branch onto main. Two adaptations were needed,
both from #44's dtype-promotion fix, which landed after this work was
written:
- The wrappers now call `_promote()` before dispatching, so the serial
fallback (nthreads <= 1) does not bypass it.
- The kernels accumulate in `values.dtype` rather than a hardcoded
float64, and reduce the thread partials into an explicitly-typed
output. `partial.sum(axis=0)` widens a narrow dtype (uint16 -> uint64),
which would have made the result dtype depend on the thread count the
machine happened to pick.
Verified that `uint16 @ uint16` still wraps exactly as numpy and scipy do,
and that the mixed `uint16 @ float64` case matches numpy to 3.7e-12
relative -- tighter than scipy's own agreement with numpy on the same input.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_019ZwuJMcTRUacxqARwjMngE
aarmey
pushed a commit
that referenced
this pull request
Sep 18, 2026
`_ops._major_matvec`/`_major_matmat` -- the scatter direction, reached by
`VCSC @ B` and `B @ VCSR` -- ran single-threaded with no strategy at all,
losing to scipy by up to 4x:
VCSC @ x 8.7 ms vs scipy 6.7 ms 1.30x
VCSC @ B (k=8) 20.1 ms vs scipy 21.0 ms 0.96x
x @ VCSR 20.5 ms vs scipy 5.1 ms 3.99x
B @ VCSR (p=8) 28.5 ms vs scipy 14.6 ms 1.96x
They cannot be parallelized over the major axis directly, since two major
slices can collide on the same output index.
Regrouping a chunk into the opposite format and running an aligned kernel --
the strategy `_vcs_matmul` uses for the normalized view -- was tried first and
is much worse here. That path caches its transposed copy and amortizes it
across an iterative algorithm's many products; a bare `A @ B` pays the regroup
once per call, which turned the 8.7 ms product into 700 ms. Reverted.
Thread-local accumulators instead -- the shape `minor_sums`/`minor_extrema`
already use -- reduced across threads afterwards. The thread count is the
whole design problem: the scatter is `nnz * width` work split across threads
while the reduction is `nthreads * n_minor * width`, so more threads is not
better. Measured on 6M nonzeros over 60k x 2k, letting the existing 64 MiB
accumulator budget pick 48 threads for a width-8 product ran it in 56.6 ms,
*slower than the serial kernel's* 23.3 ms; 4 threads and 15 MB ran it in
11.6 ms.
`scatter_threads()` takes the smaller of two bounds. Setting the derivative of
scatter-plus-reduction to zero puts the optimum at `sqrt(nnz / n_minor)`,
which is 10 for that array against a measured best of 16/8/4 at widths 1/4/8 --
right in magnitude, but blind to width. A 16 MiB byte cap supplies the missing
width-dependence, admitting fewer threads exactly as the accumulator grows.
Together they track the measured optimum across widths:
VCSC @ x 2.4 ms vs scipy 6.7 ms 0.35x (was 1.30x)
VCSC @ B (k=8) 12.0 ms vs scipy 21.6 ms 0.56x (was 0.96x)
x @ VCSR 1.8 ms vs scipy 5.2 ms 0.36x (was 3.99x)
B @ VCSR (p=8) 5.2 ms vs scipy 14.3 ms 0.36x (was 1.96x)
Peak allocation stays bounded by the cap rather than by `nnz`: 5.3 MB at
width 1, 19.2 MB at width 8, 23.0 MB at width 16. Single-thread cases fall
back to the serial kernels rather than pay an allocation and a reduction pass
for one partial.
Three benchmark cases now cover this direction; the existing
`matvec_vs_scipy`/`matmat_vs_scipy` only ever ran VCSR aligned.
Rebased from the pre-0.4.0 branch onto main. Two adaptations were needed,
both from #44's dtype-promotion fix, which landed after this work was
written:
- The wrappers now call `_promote()` before dispatching, so the serial
fallback (nthreads <= 1) does not bypass it.
- The kernels accumulate in `values.dtype` rather than a hardcoded
float64, and reduce the thread partials into an explicitly-typed
output. `partial.sum(axis=0)` widens a narrow dtype (uint16 -> uint64),
which would have made the result dtype depend on the thread count the
machine happened to pick.
Verified that `uint16 @ uint16` still wraps exactly as numpy and scipy do,
and that the mixed `uint16 @ float64` case matches numpy to 3.7e-12
relative -- tighter than scipy's own agreement with numpy on the same input.
Claude-Session: https://claude.ai/code/session_019ZwuJMcTRUacxqARwjMngE
Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
vector @ VCSRArray/VCSRArray @ vector(and the 2-D@ matrixvariants) allocated their output accumulator usingvalues.dtype— the compact stored-element dtype (e.g.uint16for single-cell counts) — instead of a dtype wide enough for the accumulated total, so results silently wrapped modulo2**bitswhilesum()(which already accumulates infloat64) stayed correct.valuesand the other operand to their commonnp.result_typebefore running the_ops.pykernels, matching the dtype promotion an equivalent dense-array product would get.Test plan
test_matmul_does_not_overflow_narrow_value_dtypeintests/test_ops.py, covering both@directions and both vector/matrix operands, for bothVCSRArrayandVCSCArray, reproducing the issue's uint16 repro.uv run pytest -q— 1662 passed, 109 skipped.uv run ruff check— clean.Fixes #43.
🤖 Generated with Claude Code