Julia: dense .+ sparse -> sparse?

[abraunst]

(reported first on this thread). I find this counterintuitive:

julia> using SparseArrays

julia> A=rand(100,100);

julia> B=sprand(100,100,0.01);

julia> typeof(A + B)
Array{Float64,2}

julia> typeof(A .+ B)
SparseMatrixCSC{Float64,Int64}

I would have expected the result to be dense. Is there a rationale for this?

[ViralBShah]

It should be dense.

[mbauman]

Coming up with a general rule here is hard. We want A .* B to be sparse, yah? Both that and A .+ B and f.(A, B) all use the same code path currently (IIRC).

[ViralBShah]

Yes, for the .*, sparse result makes sense, but if one of the inputs is dense, it is quite likely that the dense output won't blow up the memory.

[abraunst]

I have not thought very thoroughly about this, but could it make sense to make it depend on the order of the operands? I.e.

dense .op sparse -> dense
sparse .op dense -> sparse

it would be easy to remember because it would follow the behavior of .op=. The only problem I see is that sometimes you can't choose the order (i.e. when the operand is non-commutative).

EDIT: OTOH, this restriction is already present in op= or .op=, i.e. you can't really choose the ordering.

[StefanKarpinski]

Basing it on ordering seems very unJulian.

I think it makes sense to just special case .* since it’s one of the few operations where sparse/dense will be sparse.

[ViralBShah]

With ./ you get a bunch of Infs. You also have all the functions like sin, where sparse -> sparse makes sense, and cos where sparse -> dense makes sense. The key test is whether an operation on with a zero yields a zero or not.

[abraunst]

Basing it on ordering seems very unJulian.

Yeah, I guessed something like this based on the rate of thumbs down to thumbs up (infinity). :smile:

I think it makes sense to just special case .* since it’s one of the few operations where sparse/dense will be sparse.

I think it would be better to have a consistent rule instead of special casing. The point is that if we don't give the choice to the user, then it should be predictable.

[mbauman]

It gets hard with combinations of dense and sparse arrays. We can check if f(0, 0) yields a zero quite easily, but what you really want in the case of f.(dense, sparse) is to check f(nonzero, 0). But of course we can't check all nonzero values nor the general behavior and just examining one nonzero value might find the one nonzero value that ends up yielding a zero.

[ViralBShah]

Yes, it cannot be value dependent. It should be based on f(nonzero), and ideally, a commutative f(nonzero,0).