Parameter estimation of an ODE in Julia is slower than in R

General Usage Performance
, , ,
21

Could you explain what is happening here? Why is the last line necessary? What does the save_start=true option do?

  • save_start: Denotes whether the initial condition should be included in the solution type as the first timepoint. Defaults to true.

The rest I can glean from the source:

Regarding my earlier solution, I think the minimum change involves:

function DEsol_l(parms)
    u0            = @view(parms[1:n])
    integrator.p[1] = parms[n+1]
    integrator.p[2] = parms[n+2]
    reinit!(integrator, u0)
    sol           = solve!(integrator)
    sol_v         = vec(transpose(@view(sol[:,:])))   
    return sol_v
end

That seems to still return the correct answer.

22

Your solution is fine. Mine is just doing some extra tricks to keep the sol.u and sol.t caches and decrease the memory usage a bit more. It seems that memory pressure isn’t the core left in the compute though, so it doesn’t matter much.

2 Likes
23

Something additional has come up with this example. I’m trying to solve it with other methods (BFGS and Newton) just to check if I get the same values. It is not working with the BFGS method.

If I do this:

result_BFGS = optimize(b -> loglik1(AGE, HDOM, b), inits1, BFGS(), Optim.Options(iterations=1000000))

It gets solved but the minimum value is 1024.27 instead of 35.08.

If I use AD, I get an error related to the dimensions.

func_once = OnceDifferentiable(b -> loglik1(AGE, HDOM, b), inits1; autodiff=:forward);
result_once = optimize(func_once, inits1, BFGS())

I get this warning:

Warning: First function call produced NaNs. Exiting. Double check that none of the initial conditions, parameters, or timespan values are NaN.

And then:
DimensionMismatch(“arrays could not be broadcast to a common size; got a dimension with lengths 50 and 450”)

Remember 50 was the size of the vector of initial values, and 450 is what we get if we multiply the ages for which this was solved (9) by 50.

I don’t have this problem with Newton.I get the same value as with the Nelder-Mead method:

func_twice = TwiceDifferentiable(b -> loglik1(AGE, HDOM, b), inits1; autodiff=:forward);
@time result_twice = optimize(func_twice, inits1, Newton())

Why can’t I get it to work with BFGS?

24

Did you try ForwardDiff NaNSafe mode? I mention in the debugging talk this exact case:

1 Like
25

Yes. I tried that before posting. I tried with this but it didn’t solve it.

set_preferences!(ForwardDiff, "nansafe_mode" => true)

I also tried different ODE solvers to see if the problem was consistent, and it was.

26

Can you isolate it? What parameters cause the NaN? And Why? Most likely it’s parameters that are just resulting in an unstable model and doing something like reducing the initial stepnorm on BFGS will stabilize it.