broadcast!(f, C) for sparse vector/matrix C

[Sacha0]

This pull request provides broadcast!(f, C) for sparse vector/matrix C, with semantics close to those for generic two-argument broadcast!. Specifically, this implementation checks at runtime whether f() yields zero. If f() does yield zero, this implementation empties C and returns without additional f() calls. If f() does not yield zero, this implementation densifies and fills C via independent f() calls. (@stevengj, thoughts on these semantics?) Best!

[tkelman]

are you wanting to test the end values of z, fz?

[Sacha0]

Only the contents of V/fV and C/fC (which implicitly test the end values of z and fz)?

[Sacha0]

(AV x86_64 failure unrelated.)

[Sacha0]

Any thoughts on these semantics? If not, and absent objections or requests for time, I plan to merge this Monday morning PST. Best!

[martinholters]

Oh, I don't know why I missed this one, but actually, I do find the semantics a bit dubious: If f is assumed pure, independent calls seem wasteful, otherwise, just checking whether it yields zero once is insufficient. Wouldn't e.g. f() = Rand(Bool) behave funnily?

[Sacha0]

I do find the semantics a bit dubious: If f is assumed pure, independent calls seem wasteful, otherwise, just checking whether it yields zero once is insufficient. Wouldn't e.g. f() = Rand(Bool) behave funnily?

Agreed on all points in principle :). Thoughts on a (presently realizable) better compromise between the generic AbstractArray code's behavior for broadcast!(f, C) and sparse broadcast[!]'s behavior elsewhere? Thanks!

[martinholters]

I see the following options:

1. Clearly document that broadcasting a non-pure function is UB and change broadcast!(f, C) for AbstractArrays to only evaluate f once to keep people from relying on any other behavior. This would leave room for future optimizations, e.g. freedom for auto-parallelization or when broadcasting from a vector v into a matrix m, just evaluate f once for every element of v and then make copies to fill the matrix m. However, this will likely surprise people trying to fill an array with independent instances by broadcasting a constructor.
2. Always evaluate f exactly once for each element of the destination array. Consistent, least surprising, extremely wasteful for sparse arrays and pure, zero-preserving f.
3. Guess the programmers intent to do the right thing at the cost of inconsistency. If well done, will often give good performance while doing what's desired, but may probably lead to very subtle bugs in generic code.
4. Allow call-site annotation to decide between 1 and 2, e.g. by a macro that rewrites broadcast[!] and . expressions to broadcast_chosen_option[!]. Would nicely cover all cases, but requires the most effort to realize and makes the otherwise highly usable dot-notation more obscure.

None of these are completely satisfying, my gut feeling is to prefer 4.

[martinholters]

The inlining of constants makes the present state really hairy:

foo(p) = rand() < p ? rand() : 0.0

x = spzeros(100)
x .= foo.(0.1)

y = spzeros(100)
p = 0.1
y .= foo.(p)

After this, y is a sparse vector with about 10% stored non-zero values as one would expect. OTOH, x has with 90% probability no stored values, and with 10% probability 100 stored values, of which about 90% are zero.

[tkelman]

The inlining of constants makes the present state really hairy:

The fact that it's not that hard to write code where that optimization ends up changing semantics like this makes me think it's not an unequivocally good thing to be doing in all cases. Why are we doing it in lowering instead of letting usual LLVM constant propagation optimizations do their thing more generally and safely?

[Sacha0]

I see the following options:

Background reading re. broadcast!(f, A)'s semantics in general: #12277.

The inlining of constants makes the present state really hairy:

Wonderful example :).

The fact that it's not that hard to write code where that optimization ends up changing semantics like this makes me think it's not an unequivocally good thing to be doing in all cases. Why are we doing it in lowering instead of letting usual LLVM constant propagation optimizations do their thing more generally and safely?

cc @stevengj. Best!

[StefanKarpinski]

For now I would lean towards option 1: making this edge case explicitly undefined. That gives us some leeway to figure out the right thing in the next release (which might still be letting it be ub).

[stevengj]

Option 5 would be to document that broadcast assumes purity for sparse arrays only?

In algorithms with dense arrays, it's pretty useful to be able to do e.g. x .+= rand.() in stochastic algorithms, so that you don't have to allocate an auxiliary array of random numbers.

In algorithms with sparse arrays, using stochastic functions like this doesn't really make sense because they destroy sparsity, and using non-pure functions in general will defeat the ability to detect whether broadcast preserves sparsity (which is essential for performance).

[martinholters]

Option 5 would be option 3 with documentation. Getting that documentation right/exactly specifying the behavior is tricky, though. Consider

x .= f.(a, b)

For what combinations of sparse/dense a, b, and x would f be guaranteed to be evaluated exactly length(x) times? Is only the type of x of importance here? And for f.(a, b) also the return type? Is Diagonal sparse? Is UpperTriangular?

But even with this sorted out in a reasonable way, I still think that if the semantics differ by input type even though the same semantics could be applied for all types, this is very likely a severe gotcha when writing generic code.

[StefanKarpinski]

Uncertainty about this is why I'm in favor of leaving this undefined until 1.0 – hopefully we'll figure out the right choice in the next several months and then we can make that the defined behavior.

[martinholters]

Of course, we could combine options 3 (or 5) and 4: Do something reasonable by default, but allow opting in to assume-pure or assume-non-pure behavior.

x .= f.(a, b) # semantics type dependent
@assume_pure x .= f.(a,b) # may do all kinds of optimizations to reduce number of times f is called
@assume_non_pure x .= f.(a,b) # exactly length(x) calls to f in linear indexing order within a single thread

These annotations would be macros (with better names, hopefully) that do the same transformations that lowering does, but replacing broadcast! with broadcast_pure! and broadcast_non_pure!. (Similar for the non-! case, of course).