Parallelize the misaligned direction of the core matmul path - #55
Merged
Merged
Conversation
`_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
fishidaho
force-pushed
the
feat/parallel-scatter-matmul
branch
from
September 18, 2026 03:50
da540b9 to
675272c
Compare
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
2 participants
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.
Problem
_ops._major_matvec/_major_matmat— the scatter direction, reached byVCSC @ BandB @ VCSR— ran single-threaded, losing to scipy by up to 4x:VCSC @ xVCSC @ B(k=8)x @ VCSRB @ VCSR(p=8)They cannot be parallelized over the major axis directly, since two major slices can collide on the same output index. #49 fixed this for the normalized view's
_vcs_matmulpath; this is the core array path in_ops, which is still serial.Fix
Thread-local accumulators — the shape
minor_sums/minor_extremaalready use — reduced across threads afterwards.The thread count is the whole design problem: the scatter is
nnz * widthwork split across threads while the reduction isnthreads * n_minor * width, so more threads is not better. On 6M nonzeros over 60k x 2k, the existing 64 MiB accumulator budget picked 48 threads for a width-8 product and ran it in 56.6 ms — slower than the serial kernel's 23.3 ms. 4 threads ran it in 11.6 ms.scatter_threads()takes the smaller of two bounds:sqrt(nnz / n_minor), where the derivative of scatter-plus-reduction is zero, and a 16 MiB byte cap that supplies the width-dependence the first bound lacks. Single-thread cases fall back to the serial kernels rather than pay an allocation and a reduction pass for one partial.Measurements
VCSC @ xVCSC @ B(k=8)x @ VCSRB @ VCSR(p=8)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.Dtype behaviour is unchanged: the wrappers promote before dispatching so the serial fallback cannot bypass #44's fix, and the kernels accumulate in
values.dtypeso the result dtype does not depend on the thread count the machine picked.uint16 @ uint16still wraps exactly as numpy and scipy do; mixeduint16 @ float64matches numpy to 3.7e-12 relative.Three benchmark cases now cover this direction; the existing
matvec_vs_scipy/matmat_vs_scipyonly ever ran VCSR aligned.Merge order
After #53. Conflicts with the pytest-memray PR in
benchmarks/— the cases added here carrypeak_alloc_mbkeys that the other PR removes; whichever lands second drops those keys.🤖 Generated with Claude Code
https://claude.ai/code/session_01Rvu3cf8ZL7F5EPX22eo6Je