Inconsistant precision?

General Usage
1

Hello

I am a bit confuse as to how precision is assigned or indicated in Julia.
Consider:

julia> a = Complex(1.0, 0.0)
1.0 + 0.0im

julia> typeof(a)
ComplexF64

So by default a Complex number will be built with Float64 elements.
However depending on how I create a matrix of complex numbers the precision is different.

Case1: starting with an undef matrix

julia> b = Matrix{Complex}(undef, 2,2)
2×2 Matrix{Complex}:
 #undef  #undef
 #undef  #undef

julia> b[1,1] = 1.0+0.0im
1.0 + 0.0im

julia> b[1,2] = 0.0+0.0im
0.0 + 0.0im

julia> b[2,1] = 0.0+0.0im
0.0 + 0.0im

julia> b[2,2] = 1.0+0.0im
1.0 + 0.0im

julia> b
2x2 Matrix{Complex}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

julia> typeof(b)
Matrix{Complex} (alias for Array{Complex, 2})

Case2 : I create the matrix with numbers

julia> a = [Complex(1.0,0.0) Complex(0.0,0.0) ; complex(0.0,0.0) Complex(1.0,0.0) ]
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

julia> typeof(a)
Matrix{ComplexF64} (alias for Array{Complex{Float64}, 2})

So, with a scalar Julia will give me a 64 bits version whereas a matrix will be different. Is that not incoherent? Am I missing something?
If I see a function having as argument a Matrix{Complex} I have to be careful in the way I create the matrix I will be using since it may end up being Matrix{ComplexF64} and vise versa.

The behavior seems to be different with integers. Int64 is used all the time.

julia> c = Matrix{Int}(undef, 2,2)
2×2 Matrix{Int64}:
 1473840714320  1473840714384
 1473840714352  1473840714416

BTW, why do I not get #undef in this case as I got with Matrix{Complex}(undef, 2,2)?

I am using Julia 1.8.4 on Windows10

Thanks for your help clarifying how types are assigned in Julia.

2

This hasn’t got anything to do with the precision. In case 1 you create a container (matrix) parameterized by an abstract type (Complex). This makes the matrix able to contain objects of any subtype of the abstract type as elements.

Using Complex is equivalent to using something like Complex{F} where {F <: Any}.

You probably want to use Complex{Float64} instead of Complex. An alias for Complex{Float64} is ComplexF64.

1 Like
3

yet.

If I write the following function:

julia> function fun(m::Matrix{Complex})
       return nothing
       end
fun (generic function with 1 method)

And I use

julia> a = [Complex(1.0,0.0) Complex(0.0,0.0) ; complex(0.0,0.0) Complex(1.0,0.0) ]
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

julia> fun(a)
ERROR: MethodError: no method matching fun(::Matrix{ComplexF64})
Closest candidates are:
  fun(::Matrix{Complex}) at REPL[30]:1
Stacktrace:
 [1] top-level scope
   @ REPL[31]:1

Don’t you think this is confusing?
Of course I don’t mind going ComplexF64 all the way. The problem is that when I see a function I did not create asking for Matrix{Complex} I have to start converting my matrices because it will reject my Matrix{ComplexF64}.

4

This is because type parameters in Julia are invariant. I.e. Matrix{ComplexF64} is not a subtype of Matrix{Complex} even if ComplexF64 is a subtype of Complex. If anyone writes a function to operate on general complex matrices, it should not have the signature fun(m::Matrix{Complex}), but either just fun(m), or fun(m::Matrix), or fun(m::Matrix{<:Complex}) signifying a subtype of Complex.

For performance you should probably not have a matrix of type Matrix{Complex} at all. The entries can be a mixture of Complex{Int}, Complex{Float16}, Complex{Float64}, or any other Real type that somebody has made, so the compiler can’t optimize properly.

4 Likes
5

Well, that’s why one should read the Manual :grin:

Julia is more powerful than, e.g., Python, but the cost attached to that power is having to read the manual.

6

More specifically, in this case, Types · The Julia Language.

2 Likes
7

That explains my confusion.

Performance wise, if I use fun(m::Matrix{<:Complex}) , will the compiler create a function for each possible case of complex numbers? I have no problem to stay with ComplexF64, but if someone want to use it and has a huge number of matrices, he may want to use ComplexF32 instead to save memory space. So I’ll have to write a function for each case. Is that correct? Never mind matrices with mixture of different types of complex numbers, that too weird, I don’t care.

8

Yes. A function is not compiled before it’s called for the first time, with some actual parameters. It will use the types of these parameters for compiling the function. If you later call the same function with parameters of another type, it will compile an instance of the function for these new parameters too. There is nothing performance-wise to lose by writing the function with Matrix{<:Complex}. If some types of Complex numbers can benefit from another algorithm, you can even have things like fun(m::Matrix{Complex{<:Integer}}) in addition to the more general signature.

I have no problem to stay with ComplexF64, but if someone want to use it and has a huge number of matrices, he may want to use ComplexF32 instead to save memory space. So I’ll have to write a function for each case. Is that correct?

Yes. If you specify ComplexF64 the function can only be called with ComplexF64. It would be a waste to copy the same code for ComplexF32. If the function takes a Matrix{<:Complex}, it will cover both cases (in addition to Complex{Float16} and whatnot), with no loss of performance.

1 Like
9

Thanks. I now recall having read that. I forgot it.

10

And if you think the function might be used in many different contexts, like somebody calls it with fun(mat'), the adjoint of mat, which is of type Adjoint, you should consider using the signature fun(m::AbstractMatrix{<:Complex}). This will allow various subtypes of matrices as well, like Hermitian, Diagonal, Adjoint, Bidiagonal etc. Then it will work as now, and you can later add specializations for these with e.g. function fun(m::Diagonal{<:Complex}) ..., if needed.

Though with AbstractMatrix comes a Pandora’s box of indexing. Somebody might have a matrix with indices starting at -23 or some such stuff. There exists stuff for handling such things, it’s described in the manual.

1 Like
11

Do you expect your algorithm to fail for real-valued matrices? Wouldn’t that be inconvenient? Perhaps you should use something like

fun(m::AbstractMatrix)

or

fun(m::AbstractMatrix{<:Number})

?

12

Thanks for the hint @DNF.
However, in my specific case, the formulas for real numbers are much simpler and address a very specific situation. So, I will write a separate function.

13

Thanks so much @sgaure this is very relevant.

Obviously, I need to read again the manual.

14

That sounds reasonable. Then you would have two methods

fun(m::AbstractMatrix{<:Complex})
# and
fun(m::AbstractMatrix{<:Real})