0

I want to minimize an objective function with JuMP involving a bessel function. The complete function is rather complex, and I do not have a derivative of this function. It seems that JuMP during the auto-derivation calls this function with parameters equal to NaN which causes the function to fail.

Here is the full example:

using JuMP
using SpecialFunctions
using Ipopt

xdata = linspace(0.1,5,100);

data = 3 * besselk.(2,2 * xdata)

function fit4(p1,p2)
    @show p1,p2
    return sum((data - p1 * besselk.(2,p2 * xdata)).^2)
end

# is 0 as it should
@show fit4(3,2)

m = Model(solver=IpoptSolver())
@variable(m, 0.01 <= x[1:2] <= 5)
JuMP.register(m,:fit4,2,fit4; autodiff = true)

@NLobjective(m, Min, fit4(x[1],x[2]))
status = solve(m)

It produces an error:

 ERROR: AmosException with id 2: overflow.
 Stacktrace:
  [1] _besselk(::Float64, ::Complex{Float64}, ::Int32) at /home/abarth/.julia/v0.6/SpecialFunctions/src/bessel.jl:239
  [2] besselk(::Int64, ::Float64) at /home/abarth/.julia/v0.6/SpecialFunctions/src/bessel.jl:450
  [3] besselk at /home/abarth/.julia/v0.6/ForwardDiff/src/dual.jl:202 [inlined]
  [4] (::##8#10)(::ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}) at ./<missing>:0
  [5] macro expansion at ./broadcast.jl:155 [inlined]
  [6] macro expansion at ./simdloop.jl:73 [inlined]
  [7] macro expansion at ./broadcast.jl:149 [inlined]
  [8] _broadcast!(::##8#10, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2},1}, ::Tuple{Tuple{Bool}}, ::Tuple{Tuple{Int64}}, ::StepRangeLen{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}}}, ::Tuple{}, ::Type{Val{0}}, ::CartesianRange{CartesianIndex{1}}) at ./broadcast.jl:141
  [9] broadcast_t(::Function, ::Type{T} where T, ::Tuple{Base.OneTo{Int64}}, ::CartesianRange{CartesianIndex{1}}, ::StepRangeLen{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}},Base.TwicePrecision{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#fit4},Float64},Float64,2}}}) at ./broadcast.jl:270
[...]

It seems that besselk is a function from the Fortran Amos library. Does JuMP (and ForewardDiff et al.) require that all functions to be defined in pure Julia?

Alex338207
  • 1,825
  • 12
  • 16
  • 1
    In principle, `ForwardDiff` [does require](http://www.juliadiff.org/ForwardDiff.jl/stable/user/limitations.html) functions to be defined in pure Julia, yes, but there is a [rule](https://github.com/JuliaDiff/DiffRules.jl/blob/fbc4850e219ea83b1832f38d001fe0bf9210a50a/src/rules.jl#L121) in `DiffRules` for `SpecialFunctions.besselk`, so at least for the second parameter, automatic differentiation should work (and does so, when you try it in the console). – phipsgabler Mar 07 '18 at 11:05
  • @phg, this is very interesting. So in principle JuMP should be able to handle my objective function? – Alex338207 Mar 08 '18 at 12:40
  • 1
    _ForwardDiff_ should. I have actually no idea about JuMP, sorry... but I guess the error arises from something that gets passed into `besselk` by it, which is too small. Looking at the [implementation of it](https://github.com/JuliaLang/openspecfun/blob/master/amos/zbesk.f#L270), it's probably the case `CABS(Z) IS TOO SMALL`... eg., `besselk(2, 1e-200)` results in the same error. – phipsgabler Mar 08 '18 at 16:56

0 Answers0