ForwardDiff - slowed down by Real vs. Float operations

General Usage
,
1

Hello, I am using ForwardDiff to solve a minimization problem.

Following the suggestions I read, I need to specify what were my Float64 arrays as Real arrays. I do so and my code works! However, it is really slow compared to before. Basically, all I do is:

M = zeros(Real,N, N,T);

instead of

M = zeros(N, N,T);

Doing operations with M is slowing my code down by a lot. Am I doing something wrong? I see 2 options to improve my performance:

  1. Change the way I define these matrices and do the operations
  2. Extract from the ForwardDiff Dual structure the Float I want and use that in my computations.
2

Instead of
M = zeros(Real, N, N, T)
do

function foo(x)
    # (create N and T somehow, they shouldn't be globals)
    M = zeros(eltype(x), N, N, T)
    # do stuff with M
end
ForwardDiff.gradient(foo, rand(31))

It sounds like this is what you meant by:

Change the way I define these matrices and do the operations
Extract from the ForwardDiff Dual structure the Float I want and use that in my computations.

The reason it will lead to much better performance than the Real arrays is that it lets your code still be type stable. Additionally, the dual numbers (or Float64s while not differentiating) can be stored contiguously inside the array.

3 Likes
3

Real is an abstract type and thus leads to dynamic lookup for every operation done instead of fast specialised code.

May I ask why the code only works with Real, but not Float64? What’s wrong when you use Float64 instead?

4

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-time which is slowing down your code.

5

Thank you. I see. But I am using this in conjunction with optimize! from JuMP, so though ForwardDiff is being used, I don’t really have control over it. Do you suggest I still define the function and then call it when creating the matrices of zeros?

6

Hi Sukera,

yes, the problem is with ForwardDiff, it seems to be a common problem (see here).

Best,

7

It will be a lot easier to give specific advice if the code you post is closer to a mwe.

1 Like
8

No, this is not a problem with ForwardDiff.

An abstract type is named so because it doesn’t really exist (or rather: no instances of it exist, it only exists in the abstract sense). An abstract type can have multiple different concrete (that is, instances can exist) subtypes. These subtypes may have different sizes and since you’ve restricted the compiler to pretend like everything is of unknown size, the compiler has to insert checks at runtime to select the correct method for the actual size of the element.

I suggest you consult the documentation about abstract and concrete types (here).

Moreover, wrapping your code in a function like suggested by @Elrod and @bashonubuntu above allows the compiler to specialise the code according to the concrete types (and thus the known sizes of the data). This is much faster than having to ask at every use of a variable of your type.

2 Likes
9

Yes, sorry. New to this.

Here is what is hopefully an easy version of my problem:

using JuMP, Ipopt

function myfun(x,y)
    X = zeros(10,10);
    X .+= x;
    res1 = x^2 - sqrt(y);
    res2 = x + log(y);
    return maximum(abs.([res1;res2]));
end
model = Model(Ipopt.Optimizer)
JuMP.register(model, :myfun, 2, myfun, autodiff=true)

@variable(model, x >=0, start=2)
@variable(model, y >= 0, start=0.7)
@NLobjective(model, Min, myfun(x,y))
optimize!(model)

and this gives the following error:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##107#109")){typeof(myfun)},Float64},Float64,2})
Closest candidates are:
  Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194....
10

If you just replace your definition of X with

 X = fill(x, 10,10);

that should fix the issue. The general goal is to write type stable code that will work no matter what type your input is.

11

I have tried a few of the suggestions. Here is a small functional example:

My original code:

using JuMP, Ipopt

function myfun(x,y)
    X = zeros(10,10);
    X .+= x;
    res1 = x^2 - sqrt(y);
    res2 = x + log(y);
    return maximum(abs.([res1;res2]));
end
model = Model(Ipopt.Optimizer)
JuMP.register(model, :myfun, 2, myfun, autodiff=true)

@variable(model, x >=0, start=2)
@variable(model, y >= 0, start=0.7)
@NLobjective(model, Min, myfun(x,y))
optimize!(model)

and the new code following the above suggestions:

using JuMP, Ipopt

function myfun(x,y)
    X = foo(rand(2));
    X .+= x;
    res1 = x^2 - sqrt(y);
    res2 = x + log(y);
    return maximum(abs.([res1;res2]));
end
model = Model(Ipopt.Optimizer)
JuMP.register(model, :myfun, 2, myfun, autodiff=true)

@variable(model, x >=0, start=2)
@variable(model, y >= 0, start=0.7)
@NLobjective(model, Min, myfun(x,y))
optimize!(model)


function foo(x::Array{T}) where {T<:Real}
    M = zeros(T, 10, 10)
end

Both give the same error, which is:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{getfield(JuMP, Symbol("##107#109")){typeof(myfun)},Float64},Float64,2})
Closest candidates are:
  Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
  Float64(::T<:Number) where T<:Number at boot.jl:718
  Float64(!Matched::Int8) at float.jl:60
12

This means that foo takes its type from rand(2), not from x or y. Perhaps it should be foo(x)?

Another way might be X = similar(x, 10, 10) .= 0 which will make an array of the same type and element type as x, and fill it with zeros.

Or maybe x is a scalar not an array? Then X = fill(x, 10,10) is a good answer, as above, or Z = zeros(typeof(x), 10, 10); X .+= x. (Presumably in your real example the output depends on X!)

1 Like
13

I can be completely lost as I do not use this kinda of optimization, but

using JuMP, Ipopt

function myfun(x :: T, y :: T) where {T <: Real}
    X = zeros(T, 10,10);
    X .+= x;
    res1 = x^2 - sqrt(y);
    res2 = x + log(y);
    return maximum(abs.([res1;res2]));
end
model = Model(Ipopt.Optimizer)
JuMP.register(model, :myfun, 2, myfun, autodiff=true)

@variable(model, x >= 0.0, start=2.0)
@variable(model, y >= 0.0, start=0.7)
@NLobjective(model, Min, myfun(x,y))
optimize!(model)

gives

ERROR: LoadError: DomainError with -9.995576033539066e-9:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x))

and my tentative of fixing it (that may not be correct)

    res1 = x^2 - sqrt(abs(y));
    res2 = x + log(abs(y));

Seems to work.

14

Thanks! This works. Probably not the nicest solution but this is what I end up doing, actually, it’s better.

Thanks!