Julia matrix multiplication type behavior

Viewed 424

I'd like to ask about some bizarre behavior I noticed in Julia. Consider the following snippet:

function wut()
   mat::Array{Number,2}=[1 2 3; 5 6 7]
   mat2::Array{Number,2}=[1 1; 2 3]
   return typeof(mat2*mat)
end
wut()

Enormously proud of itself, this outputs Array{Any,2}. What is going on here? Why isn't the output Array{Number,2}? If one replaces Number with Float64 the output is Array{Float64,2} as it should be, but why does julia think that the multiplication of two matrices made of abstract "numbers" should come out as a matrix made up of "anything under the sun"?

3 Answers

Here's an alternative answer: this is a quirk and implementation detail that should ideally not exist. It's not a result of Julia's type design or dispatch, nor is it really a great design pattern. It's just the status quo.

Matrix multiplication is based upon an in-place API. That is, A*B becomes mul!(output_array, A, B). Thus, we need to pre-allocate the result before actually knowing what will happen. The computation of this output element type is done by a quirky and ill-specified function: promote_op. It's something that really should be removed but would require a huge and difficult refactor... and thus we have strange cases like this.

For more details, see the docstring of Base.promote_op (note that it's unexported and doesn't even appear in the online manual):

help?> Base.promote_op
  promote_op(f, argtypes...)

  Guess what an appropriate container eltype would be for storing results of
  f(::argtypes...). The guess is in part based on type inference, so can
  change any time.

  │ Warning
  │
  │  Due to its fragility, use of promote_op should be avoided. It is
  │  preferable to base the container eltype on the type of the actual
  │  elements. Only in the absence of any elements (for an empty result
  │  container), it may be unavoidable to call promote_op.

For more internal details on promote_op see issue #19669, comments in #25689, and their discussion.

The problem is that for this to happen, the compiler must be able to prove that any subtype of number multiplied by any other subtype of number is still a subtype of a number. While this is probably true, there are probably thousands of number subtypes (complex and rational both blow up the number of types), so proving this automatically is basically impossible.

First of all, declaring a Matrix with elements of type Number, does usually not make much sense. You wanna specify the Element to some Concrete number type, such as Int64 in your case, since this takes care of type stability and speeds up your code by a lot. The cases where some values of your Array need to be a Float64 and some others an Int8 are rare, and mostly you don't want to have such a scenario. Julia itself never uses on its own Abstract types, they are only needed for method definitions. I.e. you define a function for a Number type and one for a Integer type, but you will only call them with Float64 or Int64 (or any other concrete type). If you wanna know more about type hierarchy you can read the Documentation.

Anyways, what you are doing is multiplying a Matrix with elements of types Number.

A = Array{Number}(undef,1,1)

A[1] = 1

B = Array{Number}(undef,1,1)

B[1] = 1

typeof(A) 
# Array{Number,2}

typeof(A) <: AbstractMatrix
# true

So you are calling the method declared for AbstractMatrices, and you can find that in the package LinearAlgebra.jl

matprod(x, y) = x*y + x*y

function (*)(A::AbstractMatrix, B::AbstractMatrix)
    TS = promote_op(matprod, eltype(A), eltype(B))
    mul!(similar(B, TS, (size(A,1), size(B,2))), A, B)
end

So we can shrink down the to the Base.promote_op(matprod, Number, Number).

help?> Base.promote_op
  promote_op(f, argtypes...)

  Guess what an appropriate container eltype would be for storing results of    
  f(::argtypes...). The guess is in part based on type inference, so can
  change any time.
│ Warning
│
│  Due to its fragility, use of promote_op should be avoided. 
|  It is preferable to base the container eltype on the type of the actual elements. 
|  Only in the absence of any elements (for an empty result container), it may be
│  unavoidable to call promote_op.

This function then uses the Core.Compiler.return_type function, which is basically a ccall, that tries to infer to which type to promote the result. Doing this with concrete types is easy, as you can try yourself with Base.promote_op(+, Float64, Int8) which gives you a Float64. The trouble starts with Abstract types, where it does not always make much sense. As an example, this is a bit misleading:

julia> Base.promote_op(+,  Complex, Real)
Complex

julia> Base.promote_op(+,  Complex, Complex)
Complex

julia> Base.promote_op(+,  Complex, Number)
Any

julia> Base.promote_op(+,  Real, Real)
Any

So, you see, instantiating Abstract types is nothing you want to do in Maths.

Related