Why is the product of Real and Float64 arrays of eltype Any?

General Usage
1

I just discovered that for instance

julia> eltype(LinearAlgebra.Diagonal(Real[1.0]) * [1.0;;])
Any

and

julia> Real[1.0;;] * [1.0;;]
1×1 Matrix{Any}:
 1.0

Why is the eltype here Any, rather than Real?

Element-wise products behave more cautiously in widening types:

julia> Real[BigFloat(1.0), 1.0] .* [1.0]
2-element Vector{AbstractFloat}:
 1.0
 1.0

Only relevant old thread I found: Sparse matrix product causes unnecessary casting to "Any" thereby causing issues? - #11 by cjdoris

2 Likes
2

Yes, the element-wise behavior is definitely better — it’s also what you get from the comprehension, which is just based on the types of the values that you end up with:

julia> [x*y for x in Real[1.0], y in [1.0]]
1×1 Matrix{Float64}:
 1.0

julia> [x*y for x in Real[BigFloat(1.0), 1.0], y in [1.0]]
2×1 Matrix{AbstractFloat}:
 1.0 # <- this is a big float!
 1.0 # <- Float64

The answer to why things are the way they are, it’s surely a combination of status quo, simplicity of implementation, and the fact that there’s not really a huge difference between the performance/behavior of one slightly narrower abstract supertype than Any if it’s still abstract. Trying to write code that dispatches on particular abstract eltypes is often fraught.

3 Likes
3

Thanks @mbauman, that’s really helpful.

I just thought about this a bit more, and it’s harder to fix than I realised. My original thought was that obviously the product of a Real and Float64 is going to be a Real, but since Real is an abstract type and users can make new subtypes of it, that might not actually hold. I can make a weird subtype of Real of my own, and define multiplication with Float64s for it to be whatever I want. So if we preallocate the return array, the eltype has to be Any.

What e.g. comprehensions and element-wise operations do is they go through all the elements and narrow the result type based on not the eltype of the inputs but the types of the elements. This is type unstable, but results in a narrower eltype for the return value. Presumably we want to avoid doing that sort of scanning the elements for matrix products in general. You could maybe put in methods that do things differently for abstract and concrete element types, but that sounds unaesthetical.

4 Likes
4

The infrastructure to do the incremental widening of comprehensions is a little bit of magic to do something fundamentally type-unstable in a fast and efficient way (and, importantly, handle empty arrays without messing up the non-empty paths). Its implementation is indeed a little ugly… but you can use it yourself without any ugliness. Just collect a generator!

instead of:

A = zeros(size(iteration_space(#= ... =#)))
for (i,j,k) in iteration_space(#= ... =#)
    # loop body
    A[i,j,k] = #something
end

You write (backwards-ly):

collect((begin
    # loop body
    # something
end) for (i,j,k) in iteration_space(#= ... =#))
5

Yeah, that’s more or less correct. Because of how insanely polymorphic julia is, it’s often the case that knowing a value is of some abstract type tells you almost nothing about what that value can or cannot do.

This means that abstract types are often useless for inference, but still useful for code organization and the writing of generic methods.

2 Likes