Fix wrong results in sparse `prod`, `mapreduce` and `findmin`, and make `map!`/`broadcast!` safe with `missing`, fixed destinations and a throwing `f` by ViralBShah · Pull Request #1022 · JuliaSparse/SparseArrays.jl

ViralBShah · 2026-10-09T03:54:50Z

Several reductions return a value that differs from the dense one, the non-zero-preserving map/broadcast kernels throw on missing, and map!/broadcast! can leave the destination structurally invalid. Found in a soundness review of the whole package; this is the broadcast and reductions batch.

Each result is shown on Julia 1.13.1 (with the SparseArrays it ships), on main before this PR, and after this PR.

julia> prod(sparse([0.0 Inf; 1.0 2.0]))
0.0                         # 1.13
0.0                         # before this PR
NaN                         # after this PR, as dense

julia> mapreduce(x -> x == 0 ? Inf : x, +, sparse([1.0 2.0; 0.0 3.0]); dims=2)
[NaN; Inf;;]                # 1.13
[NaN; Inf;;]                # before this PR
[3.0; Inf;;]                # after this PR

julia> mapreduce(x -> x > 5, |, sparse([1.0 2.0; 0.0 3.0]); dims=1, init=true)
Bool[0 0]                   # 1.13
Bool[0 0]                   # before this PR
Bool[1 1]                   # after this PR

julia> mapreduce(x -> [x], vcat, sparse([0 1; 2 0]))
[2, 1, 0, 0]                # 1.13
[2, 1, 0, 0]                # before this PR
[0, 2, 1, 0]                # after this PR

julia> findmin(sparse([1 missing; 0 2]))
ERROR: TypeError: non-boolean (Missing) used in boolean context   # 1.13
ERROR: TypeError: non-boolean (Missing) used in boolean context   # before this PR
(missing, CartesianIndex(1, 2))                                   # after this PR

julia> A = sparse([1.0 0; 0 2.0]); map(x -> x == 0 ? 1.0 : missing, A)
ERROR: TypeError: non-boolean (Missing) used in boolean context   # 1.13
ERROR: TypeError: non-boolean (Missing) used in boolean context   # before this PR
[missing 1.0; 1.0 missing]                                        # after this PR

julia> Array(sparsevec([Inf, 0.0]) .* sparsevec([0.0, 1.0])')
[0.0 Inf; 0.0 0.0]          # 1.13
[0.0 Inf; 0.0 0.0]          # before this PR
[NaN Inf; 0.0 0.0]          # after this PR

julia> B = sparse([1.0 0 0 4; 0 2 0 0; 0 0 3 0; 1 1 0 0]);
julia> B .= ((x, y) -> x == 3 ? error("boom") : x + y).(B, sparse(ones(4, 4)))   # throws
julia> Array(B)
ERROR: AssertionError: _goodbuffers(S)   # 1.13
ERROR: AssertionError: _goodbuffers(S)   # before this PR
4×4 zeros                                # after this PR

julia> C = spzeros(Float64, Int8, 12, 12); C .= 1.0   # InexactError
julia> SparseArrays.getcolptr(C)
Int8[1, 13, 25, …, 121, 1, 1]            # 1.13
Int8[1, 13, 25, …, 121, 1, 1]            # before this PR
Int8[1, 1, …, 1]                         # after this PR, untouched

julia> Ci = sparse([1 0; 0 2]); Ci .= Ci .+ 0.5   # InexactError
julia> nnz(Ci), Array(Ci) == [1 0; 0 2]
(4, false)                               # before this PR: every entry stored, values uninitialized
(2, true)                                # after this PR, untouched

Reductions: prod keeps its zero short-circuit only when f is integer-valued, since a stored Inf or NaN times zero is NaN. Otherwise it multiplies the stored entries, and when that product is not finite it folds them into the zero instead, so a product that merely overflows (prod(sparse([0.0 1e200; 1e200 1.0]))) is still 0.0. mapreduce(f, +; dims=2) no longer adds f(0) * 0 to a row that stores every entry. The |/& column kernel folds into init instead of overwriting it. An op not known to commute is folded in element order, for a matrix, a column-subset view, a sparse vector, dims=1, dims=2 and sparse=true; runs of unstored entries are folded by repeated squaring, in the type op gives them, and f(0) is evaluated once. The fast paths remain for +, * over real or complex values, min, max, |, & and Base's internal reducers; _commutes is now one definition for mapreduce(f, op, A) and mapreduce(f, op, A, B), so the two-array form takes the fast kernel for | and & too. findmin/findmax use Base's replacement test, so missing works and a stored -0.0 ranks as in dense. any/all with a non-Boolean predicate throw TypeError, as dense does.

Map and broadcast: the non-zero-preserving kernels compare a result with the fill value by identity (!==), which fixes missing and a fill value equal to, but of another type than, the result. A result type without zero (string.(A), A .=> 1) works. x .* y' leaves the outer-product shortcut when a stored entry times zero is not zero. cos.(S) goes through dense broadcast when the index type cannot address every entry. _densestructure! checks the index type before writing, and a kernel that fills the zeros converts the fill value before it densifies the destination. A Diagonal argument keeps -0.0.

Behaviour to review: when f (or a conversion) throws inside map!/broadcast!, a destination that is not fixed is left in one of three states, never structurally invalid.

Untouched when the failure can be found first: the index type cannot address every entry, or the fill value f(0) does not convert.
Partly written, as a dense destination is, for one argument whose pattern the destination already has (C .= f.(C)). This is what main does.
Empty otherwise, since the pattern, or the filled values, are then no longer the destination's own. It is not restored. Truncating at the failing column instead needs loop state in the catch, which costs the hot loop.

Minimum of 7 runs, about 10⁶ stored entries (A is 10⁵×10⁵, x has length 10⁷); allocations are identical on every row:

Call	`main` (ms)	branch (ms)
`sum(A)`	0.109	0.115
`prod(A)`	0.000	0.136
`sum(A; dims=1)` / `dims=2`	1.011 / 2.030	0.973 / 2.092
`maximum(A; dims=2)`	2.039	2.002
`mapreduce(x -> x + 1, +, A; dims=2)`	2.016	1.991
`mapreduce(abs, (a, b) -> a + b, A)`	0.162	0.961
`mapreduce(abs, (a, b) -> a + b, x)`	0.161	0.894
`findmin(A)` / `findmax(A; dims=1)`	2.143 / 3.341	2.147 / 3.370
`map(sin, A)`	6.207	6.081
`map(+, A, B)` / `A .* B`	10.676 / 8.835	10.139 / 8.891
`map!(+, C, A, B)` / `C .= A .+ B`	9.797 / 9.769	9.981 / 9.983

Two slowdowns are the price of the correct answer: prod of a float matrix with unstored entries was O(1) and is O(nnz), and a whole-array mapreduce with an op the package does not recognise (an anonymous function, say) is about 6× slower because it folds in order.

Left out: prod of a sparse vector, and prod(...; dims), still give NaN when the stored product overflows next to a zero, as on main; x .* y' for an eltype without zero (a vector of matrices) still throws, as on main; a broadcast of two or more arguments where one stores every entry and f returns a type without zero (string.(A, sparse(ones(2,2)))) still throws, since the per-element zero test would need a tolerant check that costs every non-Number eltype; a non-associative op still differs from dense by grouping, which is outside mapreduce's contract.

The wrong-result fixes are backport candidates. The sparse-vector half of map(*, x, y) is in #1020, and unaliasing of fixed inputs in #1018. Broadcasting another shape into a fixed destination uses the helper #1033 added on main.

Beyond the routine suites: the randomized broadcast sweeps against dense from the review (about 107,000 cases) showed no failures, and the trimmed test app built and ran, both before the latest commit (the merge of main and the review fixes); neither was rerun on it.

Written by Claude Code.

🤖 Generated with Claude Code

…-finite values, order, `missing` and failed writes Reductions: - `prod` no longer returns zero when a stored entry is `Inf` or `NaN`. - `mapreduce(f, +, A; dims=2)` does not add `f(0)` to a row that stores every entry. - `mapreduce(pred, |, A; dims=1, init)` and `&` fold into `init`. - `mapreduce` with an `op` not known to commute folds the elements in their order, for a matrix, a vector, a view of columns, along a dimension and with `sparse = true`. - `findmin` and `findmax` of a matrix accept `missing`. - `any` and `all` reject a predicate that is not Boolean. `map` and `broadcast`: - The kernels for a function that is not zero at zero compare with the fill value by identity, so `missing` and a value equal to it in another type are stored. - A function whose value has no zero (`String`, `Tuple`, `nothing`, `Pair`) is accepted. - `broadcast!` of another shape into a fixed destination computes the result aside and checks it against the pattern before writing. - `x .* y'` falls back to the general kernel when a stored entry is not finite. - A dense result is not built through an index type that cannot address every entry. - `map!` and `broadcast!` leave the destination empty, and valid, when `f` or a conversion throws part-way; densifying checks the index type before writing. - A `Diagonal` argument keeps the sign of a zero on its diagonal. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

codecov · 2026-10-09T04:03:19Z

Codecov Report

❌ Patch coverage is 97.67442% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.78%. Comparing base (c0f75f0) to head (0d23e20).

Files with missing lines	Patch %	Lines
src/reductions.jl	97.81%	3 Missing ⚠️
src/higherorderfns.jl	97.10%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1022      +/-   ##
==========================================
+ Coverage   94.75%   94.78%   +0.02%     
==========================================
  Files          16       16              
  Lines        9885    10037     +152     
==========================================
+ Hits         9367     9514     +147     
- Misses        518      523       +5

☔ 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 with Codex. Reproducible findings are attached inline.

ViralBShah · 2026-10-09T13:57:46Z

+# The assertion rejects a predicate that is not Boolean, as the dense reductions do.
 Base._any(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) =
-    iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)
+    iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)::Union{Bool,Missing}


[P2] Validate predicate values before combining them in any/all

The result assertion does not catch a non-Boolean predicate value that bitwise reduction absorbs into missing:

using SparseArrays A = sparse(reshape(Union{Missing,Int}[missing, 1], 2, 1)) any(A) # missing all(A) # missing any(Matrix(A)) # TypeError all(Matrix(A)) # TypeError

This leaves the non-Boolean-input fix incomplete: missing | 1 and missing & 1 produce a final value that passes the new assertion. Please validate each predicate result that the reduction consumes, rather than only the aggregate, and apply the same change to the sparse-vector methods.

Fixed in ab82a14. The predicate is now wrapped so that each value it returns is checked to be Bool or missing before | or & sees it, in the matrix methods and the sparse-vector ones; the assertion on the aggregate is gone. any/all on the Union{Missing,Int} matrix and vector throw TypeError, as dense does, and these are tests. The wrapper costs nothing measurable: any(isnan, A) at 10⁶ stored entries is 0.10 ms on main and on the branch, and both still infer Bool.

One difference from dense remains by construction: dense any short-circuits, so any(Any[true, 1]) returns true, while the sparse reduction visits every stored value and throws.

Written by Claude Code.

Asserting the type of the result let `missing | 1` through as `missing`. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…ns-soundness # Conflicts: # src/higherorderfns.jl

ViralBShah · 2026-10-13T02:38:47Z

Closing in favour of one PR per change; links to follow here as they are opened.

ViralBShah · 2026-10-13T03:21:04Z

Split into one PR per change:

Let map and broadcast return values of a type without zero #1055 map/broadcast returning a type without zero
Throw from any and all of a sparse array on a non-Boolean value, as dense does #1056 any/all on a non-Boolean value
Fold mapreduce in element order unless op is known to commute #1063 mapreduce folds in element order unless op is known to commute
Return NaN from prod of a sparse matrix that stores an Inf or NaN next to a zero #1065 prod with a stored Inf or NaN next to a zero
Fix mapreduce(f, +, A; dims=2) adding f(0) to a row stored in full #1068 mapreduce(f, +, A; dims=2) on a row stored in full
Fold init into mapreduce with | or & along the rows #1069 mapreduce with | or & along the rows folds into init
Validate before densifying the destination of map! and broadcast! #1071 map!/broadcast!: validate before densifying the destination
Rank missing in findmin and findmax of a sparse matrix as dense does #1072 findmin/findmax with missing
Keep a non-finite value times an unstored zero in the sparse outer product x .* y' #1077 x .* y' with a non-finite value times an unstored zero
Leave a valid destination when f throws inside map! or broadcast! #1082 map!/broadcast!: valid destination when f throws (stacked on Validate before densifying the destination of map! and broadcast! #1071)
Keep the sign of a zero on a Diagonal argument of a sparse broadcast #1085 Diagonal broadcast argument keeps -0.0
Return the dense result of a broadcast whose index type cannot address every entry #1086 dense broadcast result when the index type cannot address every entry
Compare with the fill value by identity in the non-zero-preserving map and broadcast kernels #1087 non-zero-preserving kernels compare with the fill value by identity

The tests for broadcast! of another shape into a fixed destination are not carried over: test/fixed.jl already has them since #1033.

#1068) `mapreduce(f, +, A; dims=2)` returns `NaN` for a row that stores every entry when `f(0)` is not finite, where dense input gives the sum of the row. The same kernel serves `dims=1` of an adjoint or transpose and a column view. ```julia using SparseArrays f = x -> x == 0 ? Inf : x D = [1.0 2.0; 0.0 3.0] mapreduce(f, +, sparse(D); dims=2), mapreduce(f, +, D; dims=2) ``` **Before** (`main`, and Julia 1.13) ``` ([NaN; Inf;;], [3.0; Inf;;]) ``` **After** ``` ([3.0; Inf;;], [3.0; Inf;;]) ``` The `+` method of `_mapreducecols!` counts the entries each row does not store and then adds `f(0) * count` to every row. For a row with a count of zero that is `Inf * 0`, which is `NaN`. Such a row has no `f(0)` to add, so it is now skipped. The path is only reached when the matrix has some unstored entry and `f(0)` is not zero; a matrix stored in full already avoided evaluating `f(0)`. The `dims=1` kernel and the whole-array reduction go through `_mapreducezeros`, which returns early on a count of zero, and a sweep of `mapreduce(f, +, ...)` and `sum(f, ...)` over `dims` in `1`, `2`, `(1, 2)` and `:`, with and without `sparse=true`, on a matrix, its adjoint and a column view agrees with dense on this branch. Split from #1022. A candidate for backport: a wrong result, with no change to any type or signature. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…as dense does (#1056) `any` and `all` of a sparse matrix or vector reduce the predicate's values with `|` and `&` and never check that they are Boolean, so a non-Boolean value is returned as the result where dense throws a `TypeError`. Next to a `missing` it is absorbed, and the result is `missing`. **Before** ```julia julia> any([1 0; 0 0]) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> any(sparse([1 0; 0 0])), all(sparse([1 0; 0 0])), any(sparsevec([1, 0])) (1, 0, 1) julia> M = reshape(Union{Missing,Int}[missing, 1], 2, 1); julia> any(M) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> any(sparse(M)), all(sparse(M)), any(sparse(vec(M))) (missing, missing, missing) ``` **After** ```julia julia> any(sparse([1 0; 0 0])) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> all(sparse([1 0; 0 0])) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> any(sparsevec([1, 0])) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> any(sparse(M)) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> all(sparse(M)) ERROR: TypeError: non-boolean (Int64) used in boolean context julia> any(sparse(vec(M))) ERROR: TypeError: non-boolean (Int64) used in boolean context ``` `Base._any` and `Base._all` for `SparseMatrixCSCOrColumnSubset` and `SparseVectorOrView` pass the predicate straight to the `|`/`&` `_mapreduce` kernels, which accept any value those operators do. They now pass `_boolean ∘ f`, where `_boolean` lets a `Bool` or `missing` through and asserts `::Bool` on anything else, the check the dense reductions make on each value. The composed predicate still specializes: `any(iszero, A)` and `all(!iszero, A)` at 10^6 stored entries take the same time and allocations as on `main`. The reductions along a dimension (`dims=1`) are not touched here. Split from #1022. This is a bug fix and could be backported, with the caveat that calls which returned a non-Boolean value now throw. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…#1071) In-place `map!` and `broadcast!` with a function that is not zero at the zeros store every entry of the destination. They densified it first and only then hit failures that could have been found beforehand, so a throw left the destination corrupt: column pointers half rewritten, or a full pattern with uninitialized values. ```julia using SparseArrays C = spzeros(Float64, Int8, 12, 12) try; C .= 1.0; catch e; println(typeof(e)); end println(SparseArrays.getcolptr(C)) Ci = sparse([1 0; 0 2]) try; Ci .= Ci .+ 0.5; catch e; println(typeof(e)); end println(nnz(Ci), " ", Ci == [1 0; 0 2]) ``` **Before** ``` InexactError Int8[1, 13, 25, 37, 49, 61, 73, 85, 97, 109, 121, 1, 1] InexactError 4 false ``` **After** ``` InexactError Int8[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] InexactError 2 true ``` `_densestructure!` resized the buffers and wrote the column pointers before the last one, `m * n + 1`, turned out not to fit the index type, and the fill kernels called it before `fill!` converted the fill value to the element type. `_checkdensestructure` now holds what `_densestructure!` checks before it writes anything (the fixed-pattern check, and for a matrix that the last column pointer converts), and `_guardfill` runs it and the conversion of the fill value ahead of `_map_notzeropres!`, `_broadcast_notzeropres!` and the zero-argument `C .= f.()` path. The conversion is skipped for an empty destination, which stores nothing. Fixed destinations throw the same errors as before. `map!(cos, C, A)` and `C .= A .+ 1.0` at 10^6 entries take the same time and allocations as on `main`. Left out: an `f` that throws part-way through the kernel, after the destination is densified, still leaves it with values that are not its own. That is a follow-up stacked on this branch. Split from #1022. A backport candidate: the destination is corrupt after the throw. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

… does (#1072) `findmin` and `findmax` of a sparse matrix that holds a `missing` throw, or return a `NaN` that dense ranks below the `missing`. **Before** ```julia julia> findmin(sparse([1 missing; 0 2])) ERROR: TypeError: non-boolean (Missing) used in boolean context julia> findmax(sparse([NaN missing; 0.0 2.0])) (NaN, CartesianIndex(1, 1)) ``` **After** ```julia julia> findmin(sparse([1 missing; 0 2])) # as findmin([1 missing; 0 2]) (missing, CartesianIndex(1, 2)) julia> findmax(sparse([NaN missing; 0.0 2.0])) (missing, CartesianIndex(1, 2)) ``` The two comparison helpers of `_findr` spelled out the `NaN` rule with `==` and `!=`, which return `missing` next to a `missing`. They now use the replacement test of Base's dense `findmin` and `findmax`, `Base.isgreater(b, a)` and `isless(b, a)`, which are total: `missing` ranks above `NaN`, `NaN` above every number, and `-0.0` below `0.0`. Results for `NaN` and a stored `-0.0` are unchanged. A sweep of 240,000 comparisons against dense over random matrices up to 4×4 drawn from `0.0`, `-0.0`, `±1`, `±Inf`, `NaN` and `missing`, with a random subset of the zeros stored, finds no mismatch. | 10^4 × 10^4, 10^6 stored entries | `main` | branch | |---|---|---| | `findmin(A)` | 1.127 ms | 1.126 ms | | `findmax(A; dims=1)` | 1.795 ms | 1.795 ms | Split from #1022. A bug fix, so a backport candidate. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

`map` and `broadcast` over a sparse array throw a `MethodError` when the function returns a value of a type that has no `zero`, such as a `String`, a `Tuple`, a `Pair` or `nothing`. ```julia using SparseArrays A = sparse([1.0 0.0; 0.0 2.0]) string.(A) A .=> 1 map(x -> (x,), A) map(x -> nothing, A) ``` **Before** ``` string.(A) MethodError: no method matching zero(::String) A .=> 1 MethodError: no method matching zero(::Pair{Float64, Int64}) map(x -> (x,), A) MethodError: no method matching zero(::Tuple{Float64}) map(x -> nothing, A) MethodError: no method matching zero(::Nothing) ``` **After** ``` string.(A) ["1.0" "0.0"; "0.0" "2.0"] A .=> 1 sparse([1, 2, 1, 2], [1, 1, 2, 2], [1.0 => 1, 0.0 => 1, 0.0 => 1, 2.0 => 1], 2, 2) map(x -> (x,), A) sparse([1, 2, 1, 2], [1, 1, 2, 2], [(1.0,), (0.0,), (0.0,), (2.0,)], 2, 2) map(x -> nothing, A) sparse([1, 2, 1, 2], [1, 1, 2, 2], [nothing, nothing, nothing, nothing], 2, 2) ``` To choose between the zero-preserving kernel and the one that fills the result, the dispatchers evaluate `f` at the zeros and call `_iszero` on the value, which calls `zero` for anything but a `Number`. The new `_iszeroresult` is `_iszero` for a `Number`, an `AbstractArray` or `Missing`, and for any other value returns `false` when `_iszero` throws a `MethodError`: a value without a zero is not one, and fills the result like any other nonzero. It replaces `_iszero` wherever `f` at the zeros is tested: the `map!` and `broadcast!` dispatchers, `_noshapecheck_map`, the zero-argument `broadcast!`, `_diffshape_broadcast` and `_preserveszeros`, where only the call of `f` stays inside the `try`. The `try` is reached only for a non-`Number` result, so the result types of `map(sin, A)` and `A .+ B` are inferred as before. Left out: a broadcast of two or more arguments where one of them stores every entry, such as `string.(A, sparse(ones(2, 2)))`, still throws `MethodError: no method matching zero(::String)`. That path tests each computed entry for zero, and a tolerant check there would cost every eltype that is not a `Number`. Split from #1022. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…NaN` next to a zero (#1065) `prod` of a sparse matrix with an unstored entry returns zero without looking at the stored ones, but a stored `Inf` or `NaN` times zero is `NaN`, which is what the dense `prod` gives. **Before** (`main`, and Julia 1.13) ```julia julia> prod(sparse([0.0 Inf; 1.0 2.0])), prod([0.0 Inf; 1.0 2.0]) (0.0, NaN) julia> prod(sparse([0.0 NaN; 1.0 2.0]); dims=(1, 2)) 1×1 Matrix{Float64}: 0.0 ``` **After** ```julia julia> prod(sparse([0.0 Inf; 1.0 2.0])), prod([0.0 Inf; 1.0 2.0]) (NaN, NaN) julia> prod(sparse([0.0 NaN; 1.0 2.0]); dims=(1, 2)) 1×1 Matrix{Float64}: NaN ``` The `*` kernel of `Base._mapreduce` for a `SparseMatrixCSC` or a column view of one short-circuits when the unstored entries multiply to zero. That is kept only where `f` returns an integer (`Base.promote_op(f, T) <: Integer`), which has no value that survives a zero. Otherwise the zero is multiplied by the product of the stored entries. That product can overflow to `Inf` on its own and turn the zero into a `NaN` that the dense fold never sees, so when the result is nonzero the stored entries are folded into the zero one at a time instead: `prod(sparse([0.0 1e200; 1e200 1.0]))` is still `0.0`. The cost is that `prod` of a float matrix with unstored entries now reads the stored entries, O(nnz) where it was O(1). For a 2000×2000 `Float64` matrix with 10^6 stored entries: | | time | allocations | |---|---|---| | `main` | 0.00006 ms | 0 bytes | | this PR | 0.096 ms | 0 bytes | Left out: `prod` of a sparse vector and `prod(A; dims=1)` or `dims=2` go through other kernels and still give `NaN` where the stored entries of the vector or slice overflow next to a zero, as on `main` (`prod(sparsevec([0.0, 1e200, 1e200]))` is `NaN`, dense gives `0.0`). Split from #1022. This is a wrong result on released Julia, so it is a backport candidate. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…!` (#1082) When `f`, or the conversion of its result to the destination's element type, throws part-way through an in-place `map!` or `broadcast!`, the sparse destination is left with a pattern that is half-written: column pointers of the new result next to those of the old one. The next use of it fails an internal assertion or reads the wrong entries. **Before** (on #1071; the same on `main` and Julia 1.13) ```julia julia> B = sparse([1.0 0 0 4; 0 2 0 0; 0 0 3 0; 1 1 0 0]); julia> B .= ((x, y) -> x == 3 ? error("boom") : x + y).(B, sparse(ones(4, 4))) ERROR: boom julia> Array(B) ERROR: AssertionError: _goodbuffers(S) ``` **After** ```julia julia> B .= ((x, y) -> x == 3 ? error("boom") : x + y).(B, sparse(ones(4, 4))) ERROR: boom julia> Array(B) 4×4 Matrix{Float64}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 julia> nnz(B) 0 ``` The kernels rebuild the pattern of the destination, or densify it, as they go, so a throw from `f` leaves it between two patterns. The dispatchers of `map!` and `broadcast!` now call the kernel inside a `try` (`_guardpattern` for the zero-preserving kernels, the existing `_guardfill` for those that fill, and the zero-argument `C .= f.()` path), and on a throw empty the destination before rethrowing. The kernels themselves are unchanged. Behaviour to review: after a throw, a destination that is not fixed is in one of three states, never structurally invalid. - **Untouched** when the failure can be found first: the index type cannot address every entry, or the fill value `f(0)` does not convert. That is #1071. - **Partly written**, as a dense destination is, for one argument whose pattern the destination already has (`C .= f.(C)`). This is what `main` does. - **Empty** otherwise, since the pattern, or the filled values, are then no longer the destination's own. It is not restored. A fixed destination keeps its pattern throughout, as before. Truncating at the failing column instead of emptying needs loop state in the `catch`, which costs the hot loop. The `try` is outside the kernels. Minimum of 15 runs, time in ms and allocated bytes, `n = 10^5` square with about 10^6 stored entries per argument (`map!(cos, C, A)`: a 1000×1000 destination, all 10^6 entries stored): | | #1071 | this PR | |---|---|---| | `map(sin, A)` | 10.548, 16815672 | 10.458, 16815672 | | `map(+, A, B)` | 15.268, 32839224 | 15.249, 32839224 | | `A .* B` | 12.570, 32839320 | 12.599, 32839320 | | `map!(+, C, A, B)` | 12.577, 48 | 12.800, 48 | | `C .= A .+ B` | 12.563, 144 | 12.689, 144 | | `map!(cos, C, A)` | 0.207, 48 | 0.211, 48 | The out-of-place rows do not go through the new code. The inferred result types of all six are unchanged. Split from #1022. Stacked on #1071 (`vs/densify-validate-first`); only the last two commits are new. Backport: this fixes a corrupted destination, but it also changes what a destination holds after a throw (empty, where it used to be invalid but with some of its old entries readable), so it is a behaviour change as well as a fix. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…oduct `x .* y'` (#1077) `x .* y'` for sparse vectors drops a non-finite product with an unstored zero: `Inf * 0` and `NaN * 0` are `NaN` in the dense result and an unstored zero in the sparse one. `kron(x, y')` goes through the same broadcast and has the same defect. **Before** ```julia julia> x, y = sparsevec([Inf, 0.0]), sparsevec([0.0, 1.0]); julia> Array(x .* y') 2×2 Matrix{Float64}: 0.0 Inf 0.0 0.0 julia> Array(x) .* Array(y)' 2×2 Matrix{Float64}: NaN Inf 0.0 0.0 ``` **After** ```julia julia> Array(x .* y') 2×2 Matrix{Float64}: NaN Inf 0.0 0.0 ``` `Array(kron(x, y'))` changes in the same way and now equals `kron(Array(x), Array(y)')`. The outer-product shortcut `_outer` multiplies stored entries with stored entries alone, which is the product only when every stored entry times a zero is a zero. `_copy(*, x, y')` now checks that with `_outerkeepszeros`, an O(nnz) pass over both vectors, and otherwise takes the general broadcast kernel. Both branches return the same `SparseMatrixCSC` type, so the result stays inferable. `x .* y'` with 10^3 stored `Float64` entries in each vector of length 10^4 (10^6 stored results), minimum of 20 runs on a shared machine: | | time | allocations | |---|---|---| | `main` | 2.87 ms | 9 (16087432 bytes) | | this PR | 2.77 ms | 9 (16087432 bytes) | Left out: `x .* y'` for an eltype without `zero`, such as a vector of matrices, throws a `MethodError` as it does on `main`. Split from #1022. Backport candidate: a wrong result, also on Julia 1.13. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…ap` and `broadcast` kernels (#1087) `map` and `broadcast!` of a function that is not zero at zero compare each result with the fill value `f(0)` using `!=`. That throws when the comparison is `missing`, and it leaves the fill value in place of a result that equals it but has another type. **Before** ```julia julia> A = sparse([1.0 0; 0 2.0]); julia> map(x -> x == 0 ? 1.0 : missing, A) ERROR: TypeError: non-boolean (Missing) used in boolean context julia> typeof.(map(x -> x == 0 ? 1 : 1.0, A)) 2×2 Matrix{DataType}: Int64 Int64 Int64 Int64 ``` **After** ```julia julia> map(x -> x == 0 ? 1.0 : missing, A) 2×2 SparseMatrixCSC{Union{Missing, Float64}, Int64} with 4 stored entries: missing 1.0 1.0 missing julia> typeof.(map(x -> x == 0 ? 1 : 1.0, A)) # as for Array(A) 2×2 Matrix{DataType}: Float64 Int64 Int64 Float64 ``` The `_map_notzeropres!` and `_broadcast_notzeropres!` kernels (one, two and N arguments) fill the destination with `f(0)` and then skip the store of a result that equals it. They now skip only a result that is identical to it (`!==`/`===`), which is always a `Bool` and is the only case in which the store is redundant. Time at 10^6 stored entries (4000×4000, `x -> x + 1.0` and `(x, y) -> x + y + 1.0`, `Float64`, minimum of 15 runs): the one-argument kernel is unchanged (`map!` 13.0 ms on both), the two-argument kernel goes from 18.6 ms to 21.1 ms in `map!` and from 36–39 ms to 38–40 ms in `map`. Allocations are unchanged. Split from #1022. Backport candidate: a bug fix, with no new method and no change of a result type that was computed correctly before. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

#1068) `mapreduce(f, +, A; dims=2)` returns `NaN` for a row that stores every entry when `f(0)` is not finite, where dense input gives the sum of the row. The same kernel serves `dims=1` of an adjoint or transpose and a column view. ```julia using SparseArrays f = x -> x == 0 ? Inf : x D = [1.0 2.0; 0.0 3.0] mapreduce(f, +, sparse(D); dims=2), mapreduce(f, +, D; dims=2) ``` **Before** (`main`, and Julia 1.13) ``` ([NaN; Inf;;], [3.0; Inf;;]) ``` **After** ``` ([3.0; Inf;;], [3.0; Inf;;]) ``` The `+` method of `_mapreducecols!` counts the entries each row does not store and then adds `f(0) * count` to every row. For a row with a count of zero that is `Inf * 0`, which is `NaN`. Such a row has no `f(0)` to add, so it is now skipped. The path is only reached when the matrix has some unstored entry and `f(0)` is not zero; a matrix stored in full already avoided evaluating `f(0)`. The `dims=1` kernel and the whole-array reduction go through `_mapreducezeros`, which returns early on a count of zero, and a sweep of `mapreduce(f, +, ...)` and `sum(f, ...)` over `dims` in `1`, `2`, `(1, 2)` and `:`, with and without `sparse=true`, on a matrix, its adjoint and a column view agrees with dense on this branch. Split from #1022. A candidate for backport: a wrong result, with no change to any type or signature. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 763cc76) Backport adjustment: applied to src/sparsematrix.jl and test/sparsematrix_ops.jl (for src/reductions.jl and test/reductions.jl).

…#1071) In-place `map!` and `broadcast!` with a function that is not zero at the zeros store every entry of the destination. They densified it first and only then hit failures that could have been found beforehand, so a throw left the destination corrupt: column pointers half rewritten, or a full pattern with uninitialized values. ```julia using SparseArrays C = spzeros(Float64, Int8, 12, 12) try; C .= 1.0; catch e; println(typeof(e)); end println(SparseArrays.getcolptr(C)) Ci = sparse([1 0; 0 2]) try; Ci .= Ci .+ 0.5; catch e; println(typeof(e)); end println(nnz(Ci), " ", Ci == [1 0; 0 2]) ``` **Before** ``` InexactError Int8[1, 13, 25, 37, 49, 61, 73, 85, 97, 109, 121, 1, 1] InexactError 4 false ``` **After** ``` InexactError Int8[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] InexactError 2 true ``` `_densestructure!` resized the buffers and wrote the column pointers before the last one, `m * n + 1`, turned out not to fit the index type, and the fill kernels called it before `fill!` converted the fill value to the element type. `_checkdensestructure` now holds what `_densestructure!` checks before it writes anything (the fixed-pattern check, and for a matrix that the last column pointer converts), and `_guardfill` runs it and the conversion of the fill value ahead of `_map_notzeropres!`, `_broadcast_notzeropres!` and the zero-argument `C .= f.()` path. The conversion is skipped for an empty destination, which stores nothing. Fixed destinations throw the same errors as before. `map!(cos, C, A)` and `C .= A .+ 1.0` at 10^6 entries take the same time and allocations as on `main`. Left out: an `f` that throws part-way through the kernel, after the destination is densified, still leaves it with values that are not its own. That is a follow-up stacked on this branch. Split from #1022. A backport candidate: the destination is corrupt after the throw. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 4cbae1e) Backport adjustment: release-1.13 densifies with its own `_densestructure!` in src/higherorderfns.jl (without the fixed-pattern checks of main), so `_checkdensestructure` and `_guardfill` are defined there and check only the index type and the fill value. `S .= x` with a scalar goes through `fill!` and `_fillnonzero!` on release-1.13, not through `_densestructure!` as on main, so `_fillnonzero!` gets the same checks before it resizes anything. The test is a new testset in test/higherorderfns.jl that compares with `==`.

… does (#1072) `findmin` and `findmax` of a sparse matrix that holds a `missing` throw, or return a `NaN` that dense ranks below the `missing`. **Before** ```julia julia> findmin(sparse([1 missing; 0 2])) ERROR: TypeError: non-boolean (Missing) used in boolean context julia> findmax(sparse([NaN missing; 0.0 2.0])) (NaN, CartesianIndex(1, 1)) ``` **After** ```julia julia> findmin(sparse([1 missing; 0 2])) # as findmin([1 missing; 0 2]) (missing, CartesianIndex(1, 2)) julia> findmax(sparse([NaN missing; 0.0 2.0])) (missing, CartesianIndex(1, 2)) ``` The two comparison helpers of `_findr` spelled out the `NaN` rule with `==` and `!=`, which return `missing` next to a `missing`. They now use the replacement test of Base's dense `findmin` and `findmax`, `Base.isgreater(b, a)` and `isless(b, a)`, which are total: `missing` ranks above `NaN`, `NaN` above every number, and `-0.0` below `0.0`. Results for `NaN` and a stored `-0.0` are unchanged. A sweep of 240,000 comparisons against dense over random matrices up to 4×4 drawn from `0.0`, `-0.0`, `±1`, `±Inf`, `NaN` and `missing`, with a random subset of the zeros stored, finds no mismatch. | 10^4 × 10^4, 10^6 stored entries | `main` | branch | |---|---|---| | `findmin(A)` | 1.127 ms | 1.126 ms | | `findmax(A; dims=1)` | 1.795 ms | 1.795 ms | Split from #1022. A bug fix, so a backport candidate. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 0671810) Backport adjustment: applied to src/sparsematrix.jl and test/sparsematrix_ops.jl (for src/reductions.jl and test/reductions.jl).

…NaN` next to a zero (#1065) `prod` of a sparse matrix with an unstored entry returns zero without looking at the stored ones, but a stored `Inf` or `NaN` times zero is `NaN`, which is what the dense `prod` gives. **Before** (`main`, and Julia 1.13) ```julia julia> prod(sparse([0.0 Inf; 1.0 2.0])), prod([0.0 Inf; 1.0 2.0]) (0.0, NaN) julia> prod(sparse([0.0 NaN; 1.0 2.0]); dims=(1, 2)) 1×1 Matrix{Float64}: 0.0 ``` **After** ```julia julia> prod(sparse([0.0 Inf; 1.0 2.0])), prod([0.0 Inf; 1.0 2.0]) (NaN, NaN) julia> prod(sparse([0.0 NaN; 1.0 2.0]); dims=(1, 2)) 1×1 Matrix{Float64}: NaN ``` The `*` kernel of `Base._mapreduce` for a `SparseMatrixCSC` or a column view of one short-circuits when the unstored entries multiply to zero. That is kept only where `f` returns an integer (`Base.promote_op(f, T) <: Integer`), which has no value that survives a zero. Otherwise the zero is multiplied by the product of the stored entries. That product can overflow to `Inf` on its own and turn the zero into a `NaN` that the dense fold never sees, so when the result is nonzero the stored entries are folded into the zero one at a time instead: `prod(sparse([0.0 1e200; 1e200 1.0]))` is still `0.0`. The cost is that `prod` of a float matrix with unstored entries now reads the stored entries, O(nnz) where it was O(1). For a 2000×2000 `Float64` matrix with 10^6 stored entries: | | time | allocations | |---|---|---| | `main` | 0.00006 ms | 0 bytes | | this PR | 0.096 ms | 0 bytes | Left out: `prod` of a sparse vector and `prod(A; dims=1)` or `dims=2` go through other kernels and still give `NaN` where the stored entries of the vector or slice overflow next to a zero, as on `main` (`prod(sparsevec([0.0, 1e200, 1e200]))` is `NaN`, dense gives `0.0`). Split from #1022. This is a wrong result on released Julia, so it is a backport candidate. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 775677b) Backport adjustment: applied to src/sparsematrix.jl and test/sparsematrix_ops.jl (for src/reductions.jl and test/reductions.jl).

…oduct `x .* y'` (#1077) `x .* y'` for sparse vectors drops a non-finite product with an unstored zero: `Inf * 0` and `NaN * 0` are `NaN` in the dense result and an unstored zero in the sparse one. `kron(x, y')` goes through the same broadcast and has the same defect. **Before** ```julia julia> x, y = sparsevec([Inf, 0.0]), sparsevec([0.0, 1.0]); julia> Array(x .* y') 2×2 Matrix{Float64}: 0.0 Inf 0.0 0.0 julia> Array(x) .* Array(y)' 2×2 Matrix{Float64}: NaN Inf 0.0 0.0 ``` **After** ```julia julia> Array(x .* y') 2×2 Matrix{Float64}: NaN Inf 0.0 0.0 ``` `Array(kron(x, y'))` changes in the same way and now equals `kron(Array(x), Array(y)')`. The outer-product shortcut `_outer` multiplies stored entries with stored entries alone, which is the product only when every stored entry times a zero is a zero. `_copy(*, x, y')` now checks that with `_outerkeepszeros`, an O(nnz) pass over both vectors, and otherwise takes the general broadcast kernel. Both branches return the same `SparseMatrixCSC` type, so the result stays inferable. `x .* y'` with 10^3 stored `Float64` entries in each vector of length 10^4 (10^6 stored results), minimum of 20 runs on a shared machine: | | time | allocations | |---|---|---| | `main` | 2.87 ms | 9 (16087432 bytes) | | this PR | 2.77 ms | 9 (16087432 bytes) | Left out: `x .* y'` for an eltype without `zero`, such as a vector of matrices, throws a `MethodError` as it does on `main`. Split from #1022. Backport candidate: a wrong result, also on Julia 1.13. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 3d2e6df) Backport adjustment: applied by hand to src/higherorderfns.jl, where the outer product dispatches on `SparseVectorUnion` on release-1.13. release-1.13 has no `_densezero`, so the zeros are `zero(eltype(x))` and `zero(eltype(y))`, and `wrapperop` is called qualified. The test is a new testset in test/higherorderfns.jl, unguarded.

…ap` and `broadcast` kernels (#1087) `map` and `broadcast!` of a function that is not zero at zero compare each result with the fill value `f(0)` using `!=`. That throws when the comparison is `missing`, and it leaves the fill value in place of a result that equals it but has another type. **Before** ```julia julia> A = sparse([1.0 0; 0 2.0]); julia> map(x -> x == 0 ? 1.0 : missing, A) ERROR: TypeError: non-boolean (Missing) used in boolean context julia> typeof.(map(x -> x == 0 ? 1 : 1.0, A)) 2×2 Matrix{DataType}: Int64 Int64 Int64 Int64 ``` **After** ```julia julia> map(x -> x == 0 ? 1.0 : missing, A) 2×2 SparseMatrixCSC{Union{Missing, Float64}, Int64} with 4 stored entries: missing 1.0 1.0 missing julia> typeof.(map(x -> x == 0 ? 1 : 1.0, A)) # as for Array(A) 2×2 Matrix{DataType}: Float64 Int64 Int64 Float64 ``` The `_map_notzeropres!` and `_broadcast_notzeropres!` kernels (one, two and N arguments) fill the destination with `f(0)` and then skip the store of a result that equals it. They now skip only a result that is identical to it (`!==`/`===`), which is always a `Bool` and is the only case in which the store is redundant. Time at 10^6 stored entries (4000×4000, `x -> x + 1.0` and `(x, y) -> x + y + 1.0`, `Float64`, minimum of 15 runs): the one-argument kernel is unchanged (`map!` 13.0 ms on both), the two-argument kernel goes from 18.6 ms to 21.1 ms in `map!` and from 36–39 ms to 38–40 ms in `map`. Allocations are unchanged. Split from #1022. Backport candidate: a bug fix, with no new method and no change of a result type that was computed correctly before. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> (cherry picked from commit 079bcc1) Backport adjustment: the test is a new testset in test/higherorderfns.jl, with a literal matrix in place of the `fixturepair` fixture.

#1085) A `Diagonal` argument of a sparse broadcast loses the sign of a `-0.0` on its diagonal, so the result differs from the dense one. **Before** ```julia julia> sparse(ones(2, 2)) ./ Diagonal([-0.0, 1.0]) 2×2 Matrix{Float64}: Inf Inf Inf 1.0 julia> ones(2, 2) ./ Diagonal([-0.0, 1.0]) 2×2 Matrix{Float64}: -Inf Inf Inf 1.0 ``` **After** ```julia julia> sparse(ones(2, 2)) ./ Diagonal([-0.0, 1.0]) 2×2 Matrix{Float64}: -Inf Inf Inf 1.0 ``` Broadcast and `map` convert a structured argument with `_sparsifystructured`. For a `Diagonal` that went through `SparseMatrixCSC(D)`, which drops the zeros of the diagonal, and the kernel then read the dropped `-0.0` as a structural `+0.0`. `Bidiagonal`, `Tridiagonal` and `SymTridiagonal` arguments already convert with their whole band stored, zeros included, so they were not affected. `_sparsifystructured(::Diagonal)` now stores the whole diagonal as well. `SparseMatrixCSC(D)` and `sparse(D)` are unchanged. The stored pattern of a zero-preserving result does not change: the kernels store only the nonzero values they compute, so `S .* D`, `S .+ D`, `S .- D`, `map(+, S, D)` and `broadcast!(*, Z, S, D)` with zeros (of either sign) on the diagonal have the same `nnz` and pattern as on `main`. Touches dispatch: adds the internal method `_sparsifystructured(::Diagonal)`. Split from #1022. A bug fix with no change of result type, so a backport candidate. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…s every entry (#1086) A broadcast whose function is not zero at the zeros returns a dense array, but it is computed by the sparse kernels, which store every entry in a sparse array with the index type of the arguments. When that index type cannot address every entry, the call throws although the dense result is representable. ```julia S = sparse(Int8[1], Int8[1], [1.0], 12, 12) # 144 entries, `Int8` indices R = cos.(S) ``` **Before** ``` ERROR: InexactError: trunc(Int8, 133) ``` **After** ``` julia> typeof(R), R == cos.(Array(S)) (Matrix{Float64}, true) ``` `_densenonzero` now asks, from the promoted index type of the arguments and the shape of the result, whether that type addresses every entry, and otherwise hands the broadcast to `_densebroadcast`, as it already does for a `Bool` or non-concrete element type. The check is on a type and folds away, so `cos.(A)` for an `Int`-indexed matrix stays inferred as `Matrix{Float64}`; at 10^6 entries its time (2.0 ms) and allocations (16.0 MB) are the same as on `main`. Left out: the in-place forms on the same matrix (`S .= 1.0`, `map!`), which fail for the same reason in a sparse destination and are a separate change. Split from #1022. A bug fix with no change of result type, so a candidate for backporting. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

`mapreduce(pred, |, A; dims=1, init)` and `mapreduce(pred, &, A; dims=1, init)` on a sparse matrix ignore `init`: the result is the reduction of each column alone, where dense folds it into the initial value. **Before** ```julia julia> A = sparse([1.0 2.0; 0.0 3.0]); julia> mapreduce(x -> x > 5, |, A; dims=1, init=true) 1×2 Matrix{Bool}: 0 0 julia> mapreduce(!iszero, &, A; dims=1, init=false) 1×2 Matrix{Bool}: 0 1 ``` **After** ```julia julia> mapreduce(x -> x > 5, |, A; dims=1, init=true) 1×2 Matrix{Bool}: 1 1 julia> mapreduce(!iszero, &, A; dims=1, init=false) 1×2 Matrix{Bool}: 0 0 ``` The short-circuiting `_mapreducerows!` kernel for `|` and `&` assigned each column's reduction to `R[1, ii]`, which on entry holds the initial value. It now reads that value first, skips a column the initial value already decides, and stores the column's reduction otherwise. The `&` method reuses the `|` kernel with the negated predicate as before, but passes a `negate::Bool` flag in place of the `post` function, so that the kernel can also interpret the initial value. `mapreduce(>(2.0), |, A; dims=1)` on a 10^4 × 10^3 matrix with 10^6 stored entries, none of which satisfies the predicate, takes 273 μs and 1280 bytes before and 274 μs and 1280 bytes after. Split from #1022. This is a wrong result on released versions too (Julia 1.13 gives the same output), so it is a backport candidate. Written by Claude Code. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

) 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`. ```julia 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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> Co-authored-by: Viral B. Shah <ViralBShah@users.noreply.github.com>

ViralBShah commented Oct 9, 2026

View reviewed changes

ViralBShah and others added 2 commits October 9, 2026 10:12

Check every predicate value in sparse any and all

ab82a14

Asserting the type of the result let `missing | 1` through as `missing`. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Merge remote-tracking branch 'origin/main' into vs/broadcast-reductio…

0d23e20

…ns-soundness # Conflicts: # src/higherorderfns.jl

ViralBShah closed this Oct 13, 2026

Navigation Menu

Repository navigation

Fix wrong results in sparse `prod`, `mapreduce` and `findmin`, and make `map!`/`broadcast!` safe with `missing`, fixed destinations and a throwing `f` - #1022

Fix wrong results in sparse `prod`, `mapreduce` and `findmin`, and make `map!`/`broadcast!` safe with `missing`, fixed destinations and a throwing `f`#1022
ViralBShah wants to merge 3 commits into
mainfrom
vs/broadcast-reductions-soundness

ViralBShah commented Oct 9, 2026 •

edited

Loading

Uh oh!

codecov Bot commented Oct 9, 2026 •

edited

Loading

Uh oh!

ViralBShah left a comment

Uh oh!

ViralBShah Oct 9, 2026

Uh oh!

ViralBShah Oct 9, 2026

Uh oh!

ViralBShah commented Oct 13, 2026

Uh oh!

ViralBShah commented Oct 13, 2026

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

ViralBShah commented Oct 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov Bot commented Oct 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

ViralBShah left a comment

Choose a reason for hiding this comment

Uh oh!

ViralBShah Oct 9, 2026

Choose a reason for hiding this comment

Uh oh!

ViralBShah Oct 9, 2026

Choose a reason for hiding this comment

Uh oh!

ViralBShah commented Oct 13, 2026

Uh oh!

ViralBShah commented Oct 13, 2026

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

ViralBShah commented Oct 9, 2026 •

edited

Loading

codecov Bot commented Oct 9, 2026 •

edited

Loading