broadcast makes sparse matrix full

[mauro3]

julia> sp = sprand(3,3,0.1)
3x3 sparse matrix with 1 Float64 entries:
        [2, 2]  =  0.885289

julia> broadcast(*, [1,1,1], sp)
3x3 Array{Float64,2}:
 0.0  0.0       0.0
 0.0  0.885289  0.0
 0.0  0.0       0.0

[pao]

@mauro3 Which version of Julia?

[mauro3]

Sorry, both 0.3 and 0.4. (edit: oddly broadcasting .* works in 0.4 but not in 0.3: #11475)

[tkelman]

This is because of the Array on this line:

julia/base/broadcast.jl

Line 236 in 6b35206

broadcast(f::Function, As...) = broadcast!(f, Array(promote_eltype(As...), broadcast_shape(As...)), As...)

You should be able to do either broadcast!(*, copy(sp), [1,1,1], sp), or the more-specialized but non-exported SparseMatrix.broadcast_zpreserving(*, [1,1,1], sp), ref #10106

[mauro3]

So, should it be a similar there instead?

And the broadcast! trick works on 0.3 as well, thanks for the work-around. (Which is good enough for me, so no need to fix this on 0.3 from my point of view).

[tkelman]

Well, similar to which of the inputs? We'd need something to the effect of promote_containertype which we don't really have right now, and would seem to be quite difficult to do fully generically. We should try to move in that direction of being smarter and more generic about container types instead of falling back to dense Array all the time, but it's going to take some time and effort to get there.

[ViralBShah]

I am pretty sure the dense output is due to a generic fallback.

[ViralBShah]

Perhaps if one of the inputs is sparse, we should dispatch to the sparse version.

julia> @which broadcast(*, [1,1,1], sp)
broadcast(f::Function, As...) at broadcast.jl:236

julia> @which broadcast(*, sparse([1,1,1]), sp)
broadcast{Tv1,Ti1,Tv2,Ti2}(f::Function, A_1::Base.SparseMatrix.SparseMatrixCSC{Tv1,Ti1}, A_2::Base.SparseMatrix.SparseMatrixCSC{Tv2,Ti2}) at sparse/sparsematrix.jl:854

[tkelman]

Not unless we know a-priori that the function is zero-preserving?

On second thought, maybe that would make sense, but the behavior of the function on zeros would have to be checked for that case.

[toivoh]

With only one sparse input, I think it doesn't make sense to produce a
sparse output unless the function products a zero if either input is zero,
like * does (but not so many other functions).

[ViralBShah]

This seems like it should be an easy thing to test. The unfortunate thing is that we will end up with a type unstable solution. I doubt that will be a problem, but still worth keeping in mind.

[tkelman]

I don't think we can add a whole lot of cleverness to broadcast on sparse matrices in time for 0.5. There are sparse-specialized internal methods but you have to look for them and your operation needs to satisfy the invariants.

[tkelman]

@Sacha0 is there zero-preserving-detection machinery in sparse broadcast that should be kicking in on this case now? This still gives me the same behavior at 1362268

[Sacha0]

is there zero-preserving-detection machinery in sparse broadcast that should be kicking in on this case now?

The existing machinery should only kick in if all arguments are SparseMatrixCSCs. In the example, one argument is an Array. Whether broadcast over combinations of SparseMatrixCSCs and Arrays should yield an Array or a SparseMatrixCSC has not been decided yet if memory serves? Best!

[Sacha0]

Consensus seems in favor of returning an Array from broadcast when the arguments include at least one Array (alongside sparse/structured array types). As such and with sparse broadcast now functional, anything left to address in this issue? Best!

[tkelman]

where's that? for concatenation we promote combinations of dense and sparse to sparse. the potential benefit here for operations like multiplication to promote results to sparse is large. dense in sparse format is not a terrible fallback when the operation acts more like addition, and you could use broadcast! if you really wanted a dense format output

[nalimilan]

We really need a common and generic mechanism to choose the return type for cat and promote. Cf. #2326.

[Sacha0]

where's that?

Vague impression from reading comments in the week prior to that post, but that vague impression could well have been the product of my then sleep-deprivation-addled mind :).

for concatenation we promote combinations of dense and sparse to sparse.

Good point!

the potential benefit here for operations like multiplication to promote results to sparse is large.

Yes, given still more sophisticated generic sparse broadcast, though not presently. Best!

[Sacha0]

We really need a common and generic mechanism to choose the return type for cat and promote. Cf. #2326.

Agreed. Base.Broadcast.broadcast_containertype seems to be evolving into that mechanism? Best!

[nalimilan]

Agreed. Base.Broadcast.broadcast_containertype seems to be evolving into that mechanism? Best!

Yes, maybe, in which case it could be renamed to something more general. You understand these issues better than I do.

[Sacha0]

You understand these issues better than I do.

Doubtful :).

Near term (i.e. for 0.6), when broadcast'ing over combinations including both dense arrays and sparse/structured arrays, we could treat dense arrays as we do(/will, #19926) structured arrays, i.e. promote them to sparse. That change should be straightforward and would close this issue. Thoughts? Thanks!

[nalimilan]

Shouldn't the type of the result depend on whether the operation is zero-preserving? It wouldn't be great to have sp .+ a return a sparse matrix with no zero entries; better return a dense array in that case.

[tkelman]

That wouldn't be type stable though, since we can't dispatch on zero-preserving-ness from the type of a function alone.

[KristofferC]

Memory of a dense sparse matrix is only ~2x a dense matrix. User can always convert it to a full matrix.

[tkelman]

or use broadcast! or .= into a dense output

[nalimilan]

That wouldn't be type stable though, since we can't dispatch on zero-preserving-ness from the type of a function alone.

OK, sorry, I haven't followed the discussions regarding this closely enough. But Sacha didn't believe me...

[Sacha0]

The example from the original post now works (#20009):

julia> sp = sprand(3,3,0.1)
3×3 sparse matrix with 2 Float64 stored entries:
  [1, 1]  =  0.930125
  [3, 2]  =  0.781933

julia> broadcast(*, [1,1,1], sp)
3×3 sparse matrix with 2 Float64 stored entries:
  [1, 1]  =  0.930125
  [3, 2]  =  0.781933

Best!