Semantics of map! and broadcast! over sparse arrays

[Sacha0]

#19239 concluded with generic map/broadcast over sparse arrays not storing numerical zeros in the result. But whether map!/broadcast! should do the same, or instead store some numerical zeros dependent on / reflecting the inputs' structures, wasn't settled. Hence this issue.

Originally (c. #19239) I thought map!/broadcast! should retain the inputs' structure. But working on #19371 changed my thinking: map!/broadcast! should either (1) store the same result as map/broadcast or (2) retain solely the output argument's structure, perhaps throwing an error where broadcast(f, inputargs...) does not fit within the output argument's structure. Primary reason:

What "retaining the inputs' structure" means is not well defined: Should broadcast!(f, C, A, B) retain A's structure, B's structure, the union of A's and B's structures, the union of A's, B's, and C's structures, or the intersection of C's structure with the union of A's and B's structures? What if either of A or B has a singleton dimension or is scalar?

Only the output argument seems logically privileged, hence proposal (2) above. (2) has a few possible variations: When broadcast(f, inargs...) contains a nonzero value outside the equivalent broadcast!'s output argument's structure, broadcast! could either (2a) error, (2b) ignore that nonzero, or (2c) grow the output argument's structure.

Some thoughts on proposals (1) and (2):

(1) has the advantages of ease of use (it always works) and strict consistency with map/broadcast.

(2a) behaves like broadcast!(f, S, args...) where S is a structured (e.g. Diagonal, Tridiagonal, etc) matrix: strictly non-allocating, but not as flexible / easy to use as (1).

(2b) is interesting: On the one hand, it requires iteration only over the parts of the input arguments' structure that intersect the output argument's structure, potentially enabling greater efficiency where the result's structure is known/predictable. In other words, this approach provides a form of escape hatch for exploiting known detailed zero-preserving behavior of the map/broadcast function (ref. discussion in #18590), and that escape hatch could be extended to map!/broadcast! for structured matrices. On the other hand, map/broadcast and map!/broadcast! could produce inconsistent results.

(2c) seems like a half measure with little if any advantage over either (1) or (2a).

Thoughts? Thanks and best!

[martinholters]

I think silently ignoring a nonzero result is dangerous, especially in generic code. If f(x) = broadcast!(cos, similar(x), x) errors when called with a sparse x, that might be ok (unavoidable even, if x is so huge that its densified version doesn't fit into memory), but silently leaving zeros as zeros would be a nasty surprise.

[kmsquire]

Since the semantic of broadcasting over only nonzero/structural values is different than broadcasting over all elements, why not define a new set of methods for doing this, e.g., mapnz, broadcastnz?

[stevengj]

My first inclination would be to use the output argument's nonzero structure, and throw an error if this is not possible (i.e. if f(input) is nonzero at locations not stored in the output array).

I don't see the point of using broadcast! here if you don't want to pre-allocate the nonzero storage in the output. If you want to resize the output to fit f(input), then you might as well use broadcast.

[Sacha0]

Cheers, seems we can easily eliminate (2b) and (2c).

My first inclination would be to use the output argument's nonzero structure, and throw an error if this is not possible (i.e. if f(input) is nonzero at locations not stored in the output array).

My thinking as well. A slight generalization of (2a): Rather than require that the result's nonzero pattern fit within the output argument's existing structure (pattern), we could require the weaker condition that the result fit within the output argument's existing storage. Then, instead of needing to know the result's nonzero pattern, you need only know how much space it requires. Regarding implementation, this approach should also be less complex and more efficient.

I don't see the point of using broadcast! here if you don't want to pre-allocate the nonzero storage in the output. If you want to resize the output to fit f(input), then you might as well use broadcast.

Agreed. The only counterargument I see is generic programming. And on top of the variation of (2a) described above, the implementation and performance costs of allowing reallocation in broadcast! when storage runs out mid-stream may be negligible. If that bears out, providing the additional flexibility / generality / fault tolerance seems appealing?

Proposal for action: If the variation of (2a) described above (i.e. without reallocation) seems agreeable, I'll sketch an implementation of map! to get a better sense of the implementation considerations. Then revisit the reallocation question if necessary?

Thanks and best!

[tkelman]

I'm more in favor of 1. The number of nonzero entries for this kind of operation might be a bit complicated to determine ahead of time.

[Sacha0]

I'm more in favor of 1. The number of nonzero entries for this kind of operation might be a bit complicated to determine ahead of time.

Is that to say you would also be on board with the variation of (2a) described in #19372 (comment) with reallocation (which is basically (1) but only allocates if necessary)? Thanks!