[ViralBShah]

Sparse mapreduce folds the stored entries first and the unstored ones last, which gives a wrong result for an op that does not commute. This happens on main and on Julia 1.13, for a matrix, a column-subset view, a sparse vector, dims=1, dims=2 and sparse=true.

A = sparse([0 1; 2 0]); x = sparsevec([0, 1, 0, 2]); takelast = (a, b) -> b
mapreduce(x -> [x], vcat, A)                   # dense: [0, 2, 1, 0]
mapreduce(string, *, x)                        # dense: "0102"
mapreduce(identity, takelast, A; dims=1, init=7)   # dense: [2 0]
mapreduce(identity, takelast, A; dims=2, init=7)   # dense: [1; 0;;]

Before

[2, 1, 0, 0]
"0012"
[0 0]
[0; 0;;]

After

[0, 2, 1, 0]
"0102"
[2 0]
[1; 0;;]

An op that is not known to commute is now folded over the positions in element order. A run of unstored entries is folded by repeated squaring, in the type op gives it, and joined to the stored entry after it, so the cost stays O(nnz) calls of f plus O(log) calls of op per run; f(0) is evaluated once. The fast kernels remain for +, * over real or complex values, min, max, |, & and Base's internal reducers. _commutes becomes one definition, over Vararg types, shared by mapreduce(f, op, A) and mapreduce(f, op, A, B), so the two-array form now takes its fast kernel for | and & as well. With sparse=true the ordered case computes the dense reduction and stores the pattern the sparse kernels would produce.

An op the package does not recognise, such as an anonymous (a, b) -> a + b, went through Base's vectorized reduction of the stored values on main, and a strict fold in element order is 6 to 10 times slower than that. So the whole-array reduction first folds the stored values alone, pairwise and in their order, and tests in the same loop whether f(0) is a two-sided identity of op for each of them (op(z, x) === x and op(x, z) === x, in value and in type). Where it is, the unstored entries drop out wherever they lie, by associativity alone, and that fold is the result; f is called once per stored entry. Where it is not, the ordered fold runs: at once for vcat or a string product, which fail the test on the first entry, and after a wasted pass when the test fails late, as for a stored -0.0 under + (0.0 + -0.0 is 0.0).

Minimum of 30 runs, 10^6 stored entries in a 10^5 × 10^5 matrix A and in a length-10^7 vector x. The recognised reducers are unchanged:

	main	this PR
sum(A)	0.101 ms, 16 bytes	0.101 ms, 16 bytes
sum(A; dims=1)	0.952 ms, 803040 bytes	0.989 ms, 803040 bytes
sum(A; dims=2)	1.937 ms, 1605888 bytes	1.948 ms, 1605888 bytes
maximum(A; dims=2)	1.953 ms, 1607040 bytes	1.98 ms, 1607040 bytes
mapreduce(x -> x + 1, +, A; dims=2)	1.946 ms, 1605920 bytes	1.964 ms, 1605920 bytes
sum(x)	0.098 ms, 16 bytes	0.104 ms, 16 bytes

With add = (a, b) -> a + b, an unrecognised op is 1.0 to 1.7 times main where the identity test holds, no slower where it fails on the first entry (x -> x + 1), and 10 times where it fails late:

	main	this PR
mapreduce(abs, add, A)	0.155 ms, 64 bytes	0.198 ms, 64 bytes
mapreduce(identity, add, A)	0.101 ms, 64 bytes	0.162 ms, 64 bytes
mapreduce(abs, add, x)	0.154 ms, 48 bytes	0.197 ms, 48 bytes
mapreduce(identity, add, x)	0.098 ms, 48 bytes	0.162 ms, 48 bytes
mapreduce(identity, (a, b) -> max(a, b), x)	0.155 ms, 48 bytes	0.162 ms, 48 bytes
mapreduce(x -> x + 1, add, x)	7.992 ms, 48 bytes	6.57 ms, 48 bytes
mapreduce(identity, add, xm), xm being x with its last stored value -0.0	0.098 ms, 48 bytes	1.022 ms, 48 bytes

Along a dimension the kernels were not vectorized on main either, so they do not use the test. sparse=true allocates the dense reduction as well:

	main	this PR
mapreduce(abs, add, A; dims=1, init=0.0)	0.901 ms, 802944 bytes	0.982 ms, 802944 bytes
mapreduce(abs, add, A; dims=2, init=0.0)	1.939 ms, 1605792 bytes	1.793 ms, 1605792 bytes
mapreduce(abs, add, A; dims=1, init=0.0, sparse=true)	1.092 ms, 2408752 bytes	1.244 ms, 3326384 bytes
mapreduce(abs, add, A; dims=2, init=0.0, sparse=true)	2.312 ms, 3211648 bytes	2.287 ms, 3326432 bytes

A sweep of 36000 random cases (vectors, matrices, column views, twelve f/op pairs, stored zeros of either sign) agrees with foldl over the dense array in value, sign and type.

Left out: a non-associative op still differs from dense by grouping, which is outside mapreduce's contract.

Touches dispatch: _commutes (internal) goes from two type arguments to Vararg{Type,N}; no method of a Base function is added or widened.

Split from #1022.

Backport: this fixes a wrong result, but it is a large change to the reduction kernels and slows some unrecognised reducers, so it is a doubtful backport candidate.

Written by Claude Code.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 94.40559% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.78%. Comparing base (9705f28) to head (04d1460).

Files with missing lines	Patch %	Lines
src/reductions.jl	94.36%	8 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1063      +/-   ##
==========================================
- Coverage   94.79%   94.78%   -0.02%     
==========================================
  Files          16       16              
  Lines       10012    10146     +134     
==========================================
+ Hits         9491     9617     +126     
- Misses        521      529       +8     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ViralBShah]

Reviewed by Codex at 72463c7. One correctness regression found in the new ordered reduction path.

[ViralBShah]

[P2] Avoid probing op(z, z) when there is only one unstored entry

mapreduce(_ -> zeros(2, 3), *, spzeros(1)) returns the 2×3 mapped value on the base revision and for dense input, but now throws DimensionMismatch here. The reduction contains only one element, so no matrix multiplication is required; the idempotence probe introduces a multiplication that is not part of the reduction. The same problem affects isolated zero runs whose mapped value cannot multiply by itself. Defer the probe until a run actually needs at least two copies, rather than requiring every mapped zero to support op(z, z).

[ViralBShah]

Fixed in 0c491d1. _unstoredvalue now calls op(z, z) only when the counts guarantee two adjacent unstored entries (2·nstored < nslices·(len − 1)), which is exactly when the reduction calls it anyway; otherwise a run of two or more is folded by repeated squaring without the idempotence shortcut. mapreduce(_ -> zeros(2, 3), *, spzeros(1)) returns the 2×3 value again, and isolated zeros between stored entries match dense.

Written by Claude Code.