Reconsider uses of promote_op

[martinholters]

At the moment, the result's element type of broadcast is determined in a straight-forward way by typejoining the type of all the result's elements. (Exception: The empty case, where inference is invoked.) However, for certain operations, the more elaborate mechanism of promote_op is used, which especially aims at doing magical stuff if the input arrays have non-leaf element types. E.g.

julia> x = ones(Real, 1)
1-element Array{Real,1}:
 1

julia> x .+ x
1-element Array{Int64,1}:
 2

julia> x + x
1-element Array{Real,1}:
 2

While the behavior of + here does look useful, I wonder whether it is worth the inconsistency and the effort. I tend to say no, but before starting any real work on the issue, I'd like some feedback.

[martinholters]

As a quick test, I did

--- a/base/arraymath.jl
+++ b/base/arraymath.jl
@@ -68,7 +68,7 @@ promote_array_type{S<:Integer}(F, ::Type{S}, ::Type{Bool}, T::Type) = T
 
 for f in (:+, :-, :div, :mod, :&, :|, :xor)
     @eval ($f)(A::AbstractArray, B::AbstractArray) =
-        _elementwise($f, promote_eltype_op($f, A, B), A, B)
+        _elementwise($f, Broadcast._broadcast_type($f, A, B), A, B)
 end
 function _elementwise(op, ::Type{Any}, A::AbstractArray, B::AbstractArray)
     promote_shape(A, B) # check size compatibility

which at least gives

julia> x + x
1-element Array{Int64,1}:
 2

for x as above and it breaks no tests. So in case we decide to keep the current behavior, we should add tests...

[martinholters]

Poking around a bit, the code lines for mixed scalar/array cases pointed out by @andreasnoack in JuliaLang/LinearAlgebra.jl#388 are more interesting. They are responsible e.g. for

julia> 1.0 + Float32[1.0]
1-element Array{Float32,1}:
 2.0

julia> 1.0 + [1]
1-element Array{Float64,1}:
 2.0

Preserving the array element type is a no-go for the second example. Blindly using broadcast would give a Vector{Float64} in the first case. The question is whether this doubling of memory consumption would be too surprising for the casual user (also see #14252). Consistency with .+ still seems preferable to me, but here, it is maybe not as clear anymore.

[pabloferz]

Regarding the preservation of the array element type, I believe we should consider the following for a scalar s and an array A

1. Doing s ⨳ A creates a new array applying the ⨳ operation elementwise. The result of each such operation has a type which, in case of being the same for all array elements, should be the eltype of the new array, or a type wide enough to hold all the values otherwise. For the empty case, use inference. This is the behavior of broadcast.

2. If you need to preserve the type of your container you now have many options:
   a. Modify your array in place. E.g. Float32[1.0] .+= 1.0
   b. In some cases you could use an scalar such that the return type is the same as the array type. Float32[1.0] + 1.f0
   c. Use a comprehension. T=eltype(A) and T[a ⨳ s for a in A] or [T(a ⨳ s) for a in A]
   d. Take advantage of dot fusion. T.(s .⨳ A)

This and the fact that the code would be more easily maintainable, as @martinholters has mentioned before, makes me think that we should follow the route proposed in this issue (point 1 above).

[andreasnoack]

I don't think that Array+Number can meaningfully be anything but broadcast (or an error, see #5810). * is different but I'd guess that most operators are similar to + here. Consequently, the behavior should be exactly that of broadcast.

[martinholters]

I have collected some of the low-hanging fruit in this branch. I probably won't spend time on this during the holiday season, so if anyone wants to pick up from here, feel free. (Note the remarks in the commit regarding Diagonal!)

Edit: I have delete the branch as most of the work there has made it into master, and the work on Diagonal wasn't really usable as it was.

[stevengj]

I think that the performance issues with using .* etc in Diagonal should be fixed by #19639?

[martinholters]

@stevengj, is that in response to the changes in the mentioned branch? Then no, the problem is that transposition cannot be fused into broadcast, leading to additional temporaries. (broadcast_At_B, anyone? 😄)

[Sacha0]

Another voice of support for this thread's proposal and the overall direction in the branch linked above. (Re. the first commit on that branch / @stevengj's comment, ref. #19672.) Best!

[nalimilan]

This is quite convincing. I think we should just list the typical cases where preserving abstract element types could have been useful (if any), just to be sure we won't regret losing this feature too much.

[stevengj]

@martinholters, ah, I see. Once we have some kind of Transpose type (#6837 or #19670), fusion of transposition won't be a problem.

[andyferris]

I don't think that Array+Number can meaningfully be anything but broadcast (or an error, see #5810). * is different but I'd guess that most operators are similar to + here. Consequently, the behavior should be exactly that of broadcast.

I would prefer +, -, *, etc take their linear algebra meanings only (which are scalar*vector, vector+vector, and the understanding that arrays of any dimension are still living in a vector space (they are like a tensor), plus we have a variety of matrix/vector multiplications).

So Array + Number should definitely be an error, IMO (with a nice error message telling people to use .+).

[stevengj]

We already tried deprecating array + scalar and it was too inconvenient, nor is there a strong mathematical need to restrict the meaning of + to linear algebra. On the other hand, we are now deprecating most vectorized functions like sin(array), so there is greater consistency in deprecating array + scalar at this point (and potential performance advantages, but this is limited by the fact that we won't deprecate array + array or array * scalar).