# Inconsistant precision?

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.

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

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}`.

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

Well, that’s why one should read the Manual

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

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

2 Likes

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.

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

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

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

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})
``````

?

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.

Thanks so much @sgaure this is very relevant.

Obviously, I need to read again the manual.

That sounds reasonable. Then you would have two methods

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