[BUG] Wrong answer in `sparse_mul!`

[lkdvos]

Currently, the code for sparse_mul! contains a small error. Whenever there are multiple entries contributing to a single output, the beta parameter will in fact be used that many times.

SparseArraysBase.jl/src/abstractsparsearrayinterface.jl

Line 352 in aceca8c

function sparse_mul!(

MWE

using SparseArraysBase: SparseArrayDOK
a = SparseArrayDOK{Float64}(2, 2)
for I in eachindex(a)
    a[I] = 1
end

julia> a
2×2 SparseMatrixDOK{Float64, typeof(SparseArraysBase.default_getunstoredindex)}:
 1.0  1.0
 1.0  1.0

julia> a * a
2×2 SparseMatrixDOK{Float64, typeof(SparseArraysBase.default_getunstoredindex)}:
 1.0  1.0
 1.0  1.0

julia> Array(a) * Array(a)
2×2 Matrix{Float64}:
 2.0  2.0
 2.0  2.0

[mtfishman]

Yikes, good catch.

I see in SparseArrayKit.jl you just scale the destination by β at the beginning and then don't use β after: https://github.com/Jutho/SparseArrayKit.jl/blob/v0.4.0/src/tensoroperations.jl#L113, though here we also want to handle sparse arrays with non-scalar elements so maybe that isn't the best strategy for sparse_mul!.

[mtfishman]

I guess we could store a set keeping track if I_dest has been written to, and if it has we can set β to one (obviously we would have to be careful about race conditions if we multithread over elements/blocks but that's something we would have to worry about anyway once we get to that).

[lkdvos]

It's always a bit of a trade-off, I guess. In principle, if the data is really stored in a vector, scaling by a number is orders of magnitude faster than the actual multiplication, even if the entries are not scalar. Especially since the current implementation uses an allocating inner multiplication (if I'm not mistaken?), I would not really worry too much about the performance implications of separating out the scalar multiplication, unless it starts actually showing up in a profiler? Of course, if the blocks are really big and the matrices are really sparse, this could end up becoming an issue.

Keep in mind that it is not solved by simply keeping track of which I_dest has been written to by the multiplication, also the untouched elements have to be scaled correctly: the following should also work

julia> c
2×2 SparseMatrixDOK{Float64, typeof(SparseArraysBase.default_getunstoredindex)}:
 1.0  1.0
 1.0  1.0

julia> a
2×2 SparseMatrixDOK{Float64, typeof(SparseArraysBase.default_getunstoredindex)}:
 ⋅  ⋅
 ⋅  ⋅

julia> b
2×2 SparseMatrixDOK{Float64, typeof(SparseArraysBase.default_getunstoredindex)}:
 ⋅  ⋅
 ⋅  ⋅

julia> mul!(c, a, b) # beta = 0 so result should be zero!
2×2 SparseMatrixDOK{Float64, typeof(SparseArraysBase.default_getunstoredindex)}:
 1.0  1.0
 1.0  1.0

julia> mul!(Array(c), Array(a), Array(b))
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

At the end of the day, I'd say that probably the overhead most likely is larger than just looping through the memory twice.

[mtfishman]

I agree it is a subleading term in the overall scaling (I think you're right that the only limit it would matter is the very sparse element-wise limit but that limit is hard to optimize for anyway). Good point that in the most general case we need to scale the destination anyway. We can just keep it simple for now and scale the destination as a first step, and then profile and optimize as needed.