fix corner cases of binary arithmetic ops between sparse vectors and scalars (#21515)

[Sacha0]

This pull request fixes a few corner cases of binary arithmetic operations between sparse vectors and scalars. Ref. https://github.com/JuliaLang/julia/issues/21515#issuecomment-313868609 and downstream discussion. (Backport candidate? Tentatively added to the 0.6.x milestone.) Best!

[tkelman]

@nanosoldier runbenchmarks(ALL, vs=":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

[tkelman]

Hm, I thought these specializations were put in for performance reasons, but they do get the corner cases wrong. Not covered by benchmarks, or is generic broadcast just as good now?

[andreasnoack]

It probably isn't covered. I see a performance hit

julia> mulsp(x, a) = Base.Broadcast.broadcast_c(*, AbstractSparseArray, x, a)
mulsp (generic function with 1 method)

julia> x = sprand(10^6, 0.1);

julia> a = 2.3;

julia> @btime $x*$a;
  111.876 μs (5 allocations: 1.53 MiB)

julia> @btime mulsp($x,$a);
  204.978 μs (6 allocations: 1.53 MiB)

This is probably also within the scope of https://github.com/JuliaLang/julia/issues/22733. I guess the specialized broadcast methods should go in any case, though.

[Sacha0]

Hm, I thought these specializations were put in for performance reasons

IIRC you remember correctly, these specializations (and perhaps a few others like them?) were retained for performance.

Not covered by benchmarks, or is generic broadcast just as good now?

Echoing @andreasnoack, chances are the base benchmark suite does not cover these operations, and generic sparse broadcast is certainly somewhat slower.

This is probably also within the scope of #22733. I guess the specialized broadcast methods should go in any case, though.

👍 Hopefully we can continue to whittle away the remaining performance gap! :)

[Sacha0]

Other comments? Or in shape to merge? :)

[StefanKarpinski]

Bump – no comments so I assume this is good to go as soon as CI approves.

[andreasnoack]

Only half of this (removing the broadcast methods) is uncontroversial as mentioned in #22715 (comment). The other part depends on #22733

[Sacha0]

Only half of this (removing the broadcast methods) is uncontroversial as mentioned in #22715 (comment). The other part depends on #22733

Having the result of v * a differ from that of v .* a for SparseVector v and scalar a (and likewise with the order of v and a swapped, and for the analog with / and ./) seems problematic independent of the outcome of #22733. Why block this change on the outcome of #22733? Thanks!

[andreasnoack]

The reason is that v * a is a linear algebra operation and v .* a is a broadcast operation. They happen to be pretty similar but the current situation is that most/all of linear algebra follows the rule that 0*a = 0 for stored zeros whereas broadcasting doesn't follow this rule. The reason I opened #22733 was for us to make a decision about this. Before the issue has been resolved, I don't think it makes sense to make some linear algebra operation adopt the current broadcasting behavior. So to summarize, removing the specialized broadcast methods is fine and we can merge that right away but we should change the behavior for * and / before #22733 is closed.

[Sacha0]

Cheers, I will strip the changes to v * a and a * v from this pull request for now, the results of those operations deviating from those for v .* a and a .* v at most in stored entry structure where a is neither Inf nor NaN. The behavior of v / a where a is zero (the bug reported in https://github.com/JuliaLang/julia/issues/21515#issuecomment-313868609) concerns me more though. Would you be on board with at least fixing that behavior immediately? Best! :)

[andreasnoack]

I think we'll have to keep the current behavior for / for the same reason. For the symbolic computation to be useful, you'd have to assume finite α and you'd only be allowed to check after the symbolic computation. It would also be weird that the result of x/a was different from x*inv(α) (except for rounding).

[tkelman]

what symbolic computation? that's not how our api's work, with a handful of exceptions (not related to these specific functions) in internals that call cholmod or umfpack but aren't public. if we had split symbolic-numeric apis in more places, then the numeric step would be responsible for verifying that the input data satisfies the assumptions the symbolic step used. but that's not how these work.

[Sacha0]

The ability to perform structural/symbolic computations (here, for example, pattern intersection or union calculation) certainly has value and merits consideration. Such functionality presently exists implicitly in some operations, if somewhat ad hoc and inconsistently. Better (and explicitly) supporting such functionality would be lovely!

On the one hand, improving/extending the existing implicit support would constrain/compromise the semantics, implementation, and performance of both a variety of operations over sparse/structured objects (for example [most?] arithmetic operations, and associated higher order functions) as well as the embedded, implicit structural/symbolic functionality. This thread is a case in point.

On the other hand, providing separate, dedicated functionality for structural/symbolic computations would obviate the above constraints/compromises. Additionally, such separation should lead to greater clarity and consistency in semantics, documentation, and code, and potentially better reusability and composability. Hence in part the usual separation of such functionality I wager.

As such, providing separate, dedicated functionality for structural/symbolic computations strikes me as a better approach than continuing to yoke other operations with that responsibility.

Thoughts? Best! :)

[andreasnoack]

We can use a symbolic/numerical decoupling in implementations without exposing a symbolic API. So roughly speaking there are three solutions for each function: no decoupling (which is how most of the sparse broadcasting works), implicit decoupling (which is how most of the sparse linear algebra works) or explicit decoupling where we expose a symbolic API. #22733 is the place for discussing pros and cons. For now, the status is that broadcasting and linear algebra have different behavior. Hence, / shouldn't call ./.

[tkelman]

Doesn't correctness matter more than imposing an artificial discrepancy between linear algebra functions and broadcasting functions when the dimensions are such that they should behave identically? 0/0 == 0 isn't very defensible behavior.

Even if there were implicit numeric/symbolic calculation stages in the internal implementations of sparse linear algebra functions, it would be fairly straightforward to branch early and call the numeric-symbolic implementation when its assumptions are satisfied by the input data, and call more general consistent-with-dense fallbacks when the assumptions aren't satisfied.

For now, the status is that broadcasting and linear algebra have different behavior.

And part of the difference is a correctness bug that's being fixed here.

[Sacha0]

Even if there were implicit numeric/symbolic calculation stages in the internal implementations of sparse linear algebra functions, it would be fairly straightforward to branch early and call the numeric-symbolic implementation when its assumptions are satisfied by the input data, and call more general consistent-with-dense fallbacks when the assumptions aren't satisfied.

Sparse broadcast follows this design :).

[Sacha0]

I stripped the removal of the contentious non-broadcast arithmetic specializations. Thoughts? Thanks!

[Sacha0]

Thanks all! :)