Initialise array without specific types

[baptnz]

I wrote a program that works fine but has several hard-coded types (Float64, etc.) which I’m now trying to get rid of, because they get in the way of automatic differentiation.

I’ve hit a snag with the following: I initialise a matrix before a for loop, and then fill it one row at a time. My original initialisation step was

cext = Array{Float64}(undef,(3,4))

for (...)

        cext[ii,:] = ...

end

How would I got about this, without specifying the type Float64 in the constructor? I don’t mind if the results are temporarily stored in a different container (vector of rows), but again I’m not sure of the syntax to use (Vector(undef,3) works, but I’m not sure if it’s a good idea as the size isn’t meaningfully allocated). The matrix isn’t too big, rows and columns in the hundreds typically.

[rdeits]

If you know that your matrix will have a type which matches some function input, then you can use that:

julia> function build_matrix(x)
         out = Matrix{typeof(x)}(undef, 2, 2)
         for i in 1:2
           out[:, i] .= x
         end
         out
       end
build_matrix (generic function with 1 method)

julia> build_matrix(1)
2×2 Matrix{Int64}:
 1  1
 1  1

julia> build_matrix(1.0)
2×2 Matrix{Float64}:
 1.0  1.0
 1.0  1.0

or, equivalently:

julia> function build_matrix(x::T) where {T}
         out = Matrix{T}(undef, 2, 2)
         ...

Alternatively, if you have easy access to the individual rows, you can do:

cext = permutedims(reduce(hcat, rows))

where rows is a vector of vectors.

By the way, you will get better performance if you transpose your data so that you are frequently accessing individual columns of the matrix rather than rows. This is because Julia arrays are column-major, just like fortran and matlab but unlike C and numpy.

[bashonubuntu]

@baptnz

Array{Float64}(undef,(3,4))

The array needs to be a subtype of Real (but no smaller) for automatic differentiation to work.

I think you need something like

julia> function fill_mat_new(x::T) where {T<:Real} 
           Array{T}(undef,(3,4)) 
end

fill_mat_new (generic function with 1 method)

julia> fill_mat_new(2)
3×4 Array{Int64,2}:
 0  0  0  0
 0  0  0  0
 0  0  0  0

julia> fill_mat_new(2.3)
3×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0  
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

Here is a quick example which may be helpful

function f_new(meta::Array{T}, init) where {T<:Real}
    params = Array{T}(undef, size(init)[1], 1)
    for i=1:5
        params[i] = init[i]
    end
    for i=2:5
        params[i] = meta[1] * params[i-1]
    end
    return sum(params)
end

s = [1.,2.,3.,4.,5.]
ForwardDiff.gradient(x -> f_new(x, s), [2.0])

julia> ForwardDiff.gradient(x -> f_new(x, s), [2.0])
1-element Array{Float64,1}:
 49.0

You may also find these threads helpful

ForwardDiff cannot assign ForwardDiff.Dual to input array

Try using eltype() to debug. This often how I catch errors when I’m using ForwardDiff in my code. It’s likely that one of the many parameters might be having its type forced into something like an array of floats. Same code using eltype() below. Notice that special ForwardDiff.Dual type being created. This is why you want to avoid forcing the type of a parameter to be Float64 because Julia can’t convert ForwardDiff.Dual type to a Float64, hence the "“cannot convert ForwardDiff.Dual to Float64 e…

ForwardDiff - slowed down by Real vs. Float operations

Building off on @Elrod’s suggestion, note that you can also do julia> function foo(x::Array{T}) where {T<:Real} M = zeros(T, 1, 1, 1) # do stuff with M end foo (generic function with 1 methods) julia> ForwardDiff.gradient(foo, rand(31)) Note that parametrizing the type of the input function will also lead to the kind of specialized code which I think @Sukera is hinting at in his reply. The current issue seems to be that the conversion to Real is happening at run-t…

[Tamas_Papp]

@baptnz

they get in the way of automatic differentiation

The first thing I try in these situations is work around this by using a functional style (mapreduce etc) that figures out the type itself. But without an MWE it is hard to suggest something specific.

Please read: make it easier to help you

Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

[baptnz]

@bashonubuntu

function fill_mat_new(x::T) where {T<:Real} 
           Array{T}(undef,(3,4)) 
end

This sounds like a great solution but I wonder how I can apply it to a variable that is not in the function arguments? Must I define a constructor like fill_mat_new, or is there a way to signal types in the case below:

function do_something(x::Int) 

    # ... 
    
        # intermediate_thing = Array{T}(undef,(3,4)) where {T<:Real} # fails
        intermediate_thing = Array{Real}(undef,(3,4)) 

    # ...
    
    intermediate_thing[1,1] = 3.0
    
    # ...

   return(intermediate_thing[1,1] * x)
end

do_something(1)

[baptnz]

@Tamas_Papp

The first thing I try in these situations is work around this by using a functional style ( mapreduce etc) that figures out the type itself. But without an MWE it is hard to suggest something specific.

Thanks – I usually go for functional style but this particular code is very matrix-centric (it’s all about solving a linear system) and it’s easier to think about it by filling the matrix in a nested for loop. The matrix is also quite large, and gets reused multiple times by replacing values in-place. For these two reasons I chose to pre-allocate the matrix at the start. If I can’t work the other suggestion(s) I’ll try to make the problem into a minimal example more representative of the actual problem.

[bashonubuntu]

You can do something like and pass an instance of Array{T}(undef, 0, 0), like Array{Real}(undef, 0, 0) in the y argument.

function do_something(x::Int, y::Array{T}) where {T<:Real} 

  
    intermediate_thing = Array{T}(undef,(3,4))

    intermediate_thing[1,1] = 3.0
    
   return(intermediate_thing[1,1] * x)
end

julia> do_something(1, Array{Real}(undef, 0, 0))
3.0

julia> do_something(1, Array{Int}(undef, 0, 0))
3

julia> do_something(1, Array{Float64}(undef, 0, 0))
3.0

Notice how the output changes when I change the instance of y that I pass into do_something(x, y)

[rdeits]

@bashonubuntu

You can do something like and pass an instance of Array{T}(undef, 0, 0) , like Array{Real}(undef, 0, 0) in the y argument.

Why wouldn’t you just pass the type T instead of an unused array?

[bashonubuntu]

Yes, you can do that as well. Just picked Array{T}(undef, 0, 0) arbitrarily.

[rdeits]

@baptnz

Must I define a constructor like fill_mat_new , or is there a way to signal types in the case below:

It’s hard to give an answer that will fit your use case without knowing more about what you actually want. The answer in the example you’ve given is easy because it’s just Float64, but presumably that’s not the case in your real code. But the data you’re filling the matrix with must come from somewhere, and there may be an easy way to pre-allocate the matrix if you can tell us something about where that somewhere is.

[baptnz]

@rdeits

If you know that your matrix will have a type which matches some function input, then you can use that:

That was my thought – thanks for providing the example; in my case things get trickier (maybe?) because some of the intermediate quantities (not the input or the output), such as the elements of this matrix, are complex numbers, and they tend to cause problems with AD. So I’m not entirely sure how to go about making

 out = Matrix{ Complex{typeof(x)} }(undef, 2, 2)

would that make sense for, say, dual numbers? It seems to work with Float input, just wanted to check I’m doing something sensible.

Actually one of those matrices should be a unit matrix, so presumably I could even do

Matrix{Complex{typeof(x)}}(I, 2, 2)

[baptnz]

OK thanks, but do I understand correctly that one way or another the where {T<:Real} can only be applied to one of the arguments of the function? If my only argument was x, and not of type T, I would need to add another just to specify the type T?

[bashonubuntu]

Here’s a more complete example with ForwardDiff working. This also takes care of the comment by @rdeits. I had actually assumed that @baptnz wants to eventually pass this into ForwardDiff.gradient, so we would need to pass an Array as an argument.

julia> function do_something(x::Array{T}) where {T<:Real}

          intermediate_thing = Array{T}(undef,(3,4))
              intermediate_thing[1,1] = 3.0

         return(intermediate_thing[1,1] * x)[1]
      end
do_something (generic function with 1 method)

julia> ForwardDiff.gradient(x -> do_something(x), [2])
1-element Array{Int64,1}:
3

julia> ForwardDiff.gradient(x -> do_something(x), [4])
1-element Array{Int64,1}:
3

Now notice when I change from T<:Real to T<:Float64 the ForwardDiff call will fail

julia> function do_something_2(x::Array{T}) where {T<:Float64}

           intermediate_thing = Array{T}(undef,(3,4))
               intermediate_thing[1,1] = 3.0

          return(intermediate_thing[1,1] * x)[1]
       end
do_something_2 (generic function with 1 method)

julia> ForwardDiff.gradient(x -> do_something_2(x), [4.0])
ERROR: MethodError: no method matching do_something_2(::Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#509#510",Float64},Float64,1},1})

julia> ForwardDiff.gradient(x -> do_something(x), [4.0])
1-element Array{Float64,1}:
 3.0

And, this failure occurs because

julia> ForwardDiff.Dual{ForwardDiff.Tag{var"#509#510",Float64},Float64,1} <: Real
true

julia> ForwardDiff.Dual{ForwardDiff.Tag{var"#509#510",Float64},Float64,1} <: Float64
false

So, you need the array type to be no smaller (in terms of subsets) than Real

[baptnz]

The code’s logic is to have a high-level function where this array is pre-allocated, and then will be re-used in place many times by some internal functions, a bit like the example below,

function update_matrix!(m)

  x = [1.2, 3.4]
  m[1,1] = exp(1im*x[1]) + sin(x[2])
  m[2,1] = exp(1im*x[2]) + sin(x[2])
  m[1,2] = exp(1im*x[1]) + sin(x[1])
  m[2,2] = exp(1im*x[2]) + sin(x[1])

end

function high_level(a::Int, x::Real)

  m = Matrix{Complex{Real}}(undef, 2,2)

  for i in 1:a
    update_matrix!(m)
  end

  return(real(m[1,1])*x)
end

high_level(2,1.5)

So in a sense the high-level function has little idea of what the type of m should be, as it’s left to the update_matrix function to deal with it.

I could obviously refactor the code and have everything in a monolithic function, but it’s much nicer with the current structure (and I’ve seen a few viable options already in this thread to pre-allocate m – in fact I think this example is just fine now, with m = Matrix{Complex{typeof(x)}}(undef, 2,2) and I probably don’t need something more general than that).

[baptnz]

I think I get it more or less, but how would I pass the type in practice? (the code below fails)

# function do_something(x::Int, T) where {T<:Real} # fails
function do_something(x::Int, y::T) where {T<:Real} 

  
    intermediate_thing = Array{T}(undef,(3,4))

    intermediate_thing[1,1] = 3.0
    
   return(intermediate_thing[1,1] * x)
end

do_something(1.2, Float16) # also fails

[bashonubuntu]

You need to pass an instance (or, example) of a type, so in this case any real number would work like so

function do_something_3(x::Int, y::T) where {T<:Real} 

  
    intermediate_thing = Array{T}(undef,(3,4))

    intermediate_thing[1,1] = 3.0
    
   return(intermediate_thing[1,1] * x)
end

julia> do_something_3(1, 1.2)
3.0

Also, note you’re making another mistake when you write do_something(1.2, Float16). The argument x in do_something(x, y) needs to be an integer and not a float. For example, this fails (notice the stacktrace)

julia> do_something_3(1.0, 1.2)

ERROR: MethodError: no method matching do_something_3(::Float64, ::Float64)

Closest candidates are:
  do_something_3(::Int64, ::T) where T<:Real at REPL[293]:1

but this doesn’t

julia> do_something_3(1, 1.2)
3.0

[baptnz]

Ah OK, thanks (and the 1.2 was a typo from copy-paste). I see how this can work but I’ll probably try to stick to the other solution (inferring the type from other variables) because I prefer not to have dummy variables in the function signature, especially for a high-level function (that would confuse other users of my code who shouldn’t need to worry about implementation details).

[bashonubuntu]

You don’t need a function signature in any of these cases, btw. You can leave the function completely unsigned. You can check by removing all the function signatures in the examples I gave.

One of the beauties of Julia is that you can write completely generic code, and signing functions like we did doesn’t necessarily help with speed (but can sometimes help in figuring out errors like we did or can help catch mistakes during compile time). I personally find it helpful to remember what the type of each argument is, but it’s a style thing.

@baptnz But I don’t understand why you call it a “dummy” type. You will be passing x anyway, so you don’t have to create an extra argument just to pass an instance of a type (in case that’s your worry). I did give an example of this above too.

julia> function do_something(x::Array{T}) where {T<:Real}

          intermediate_thing = Array{T}(undef,(3,4))
              intermediate_thing[1,1] = 3.0

         return(intermediate_thing[1,1] * x)[1]
      end
do_something (generic function with 1 method)

julia> ForwardDiff.gradient(x -> do_something(x), [2])
1-element Array{Int64,1}:
3

[rdeits]

@bashonubuntu

You need to pass an instance (or, example) of a type, so in this case any real number would work like so

This isn’t correct. You don’t need to pass an “example” of a type if all you want is the type itself:

julia> function foo(::Type{T}) where {T}
        return Matrix{T}(undef, 2, 2)
      end
foo (generic function with 1 method)

julia> foo(Int)
2×2 Matrix{Int64}:
0  0
0  0

julia> foo(Float64)
2×2 Matrix{Float64}:
0.0  0.0
0.0  0.0

[rdeits]

@baptnz

in fact I think this example is just fine now, with m = Matrix{Complex{typeof(x)}}(undef, 2,2) and I probably don’t need something more general than that).

That actually looks great. The compiler should do a good job of figuring out what Complex{typeof(x)} means, and you should end up with code that runs efficiently.

This is the point I was trying to get across about where your data comes from. In your case, you know that the element type will be the complex version of the type of x, so Complex{typeof(x)} solves the problem perfectly. In other cases you might know something like “the element type is some number which can hold the type of x or y”, in which case you could do promote_type(typeof(x), typeof(y)) or even something like 'the element type is the result of dividing x by y, in which case you can do Base.promote_op(/, typeof(x), typeof(y)).