Levenberg Marquardt Performance Regression on Master

[ChrisRackauckas]

Maybe someone can come up with a more minimal example. The problem I was investigating was parameter esitmation for ODEs. I setup an ODE and generated the data knowing that a=1.5:

using ParameterEstimation, OrdinaryDiffEq, ParameterizedFunctions, LsqFit,
      DiffEqBase
using Base.Test

f = @ode_def_nohes LotkaVolterraTest begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a=>1.5 b=1.0 c=3.0 d=1.0

u0 = [1.0;1.0]
tspan = [0;10.0]

prob = ODEProblem(f,u0)
sol = solve(prob,tspan)

using Plots
gr()
plot(sol)

t = collect(linspace(0,10,200))
randomized = [(sol(t[i]) + .01randn(2)) for i in 1:length(t)]

data = Matrix{Float64}(length(randomized),length(u0))
for i in 1:length(randomized)
  data[i,:] = randomized[i]
end
scatter!(t,data)

scatter(t,data)

f = LotkaVolterraTest(a=1.42)
prob = ODEProblem(f,u0)
sol = solve(prob,tspan)
plot!(sol)

From this we can build a model to optimize:

function model(t,p)
  for i in eachindex(f.params)
    setfield!(f,f.params[i],p[i])
  end
  sol = solve(prob,tspan)
  vecout = sol(t)
  y = Matrix{Float64}(length(t),length(u0))
  for i in 1:length(t)
    y[i,:] = vecout[i]
  end
  vec(y)
end

The Levenburg Marquardt implementation isn't very robust, so it can't stray very far from the true answer and actually get it right:

fit = curve_fit(model,t,vec(data),[1.45],show_trace=true,lambda=1000.0)

param = fit.param

That gets 1.50 on v0.6.1. However, on master it gets 0.04, i.e. it diverges. When the initial value is instead 1.49, it goes to 0.132. Thus on master it diverges even when starting right next to the true answer.

The only commit in this timeframe seems to be f0f5c18 . I think it should be reverted or investigated more deeply.

[pkofod]

I guess we will have to have a look at those changes once more. Reverting will reintroduce a bug where x could stray out of bounds, so that might not be optimal either. @matthieugomez @bjarthur

[matthieugomez]

I cannot reproduce your example because I cannot download the package you are using. Are you using bounds when calling levenberg_marquardt?

[ChrisRackauckas]

ParameterEstimation.jl is waiting to be registered (the PR is in METADATA), but can be cloned if you want to test this (there are a lot of dependencies though since it's pulling in a lot muscle from JuliaDiffEq). I didn't set bounds: I just increased the lambda since at the default lambda it really didn't work well (on both master and the latest release). How does one set the bounds for this? I can try it with different solver parameters and see what turns up.

Note I couldn't find much documentation on using this. It seems LM is left out of the Optim.jl docs (so I went to the code and found it had a docstring). It should really either be documented here, or moved out to LsqFit and documented there.

[pkofod]

Note I couldn't find much documentation on using this. It seems LM is left out of the Optim.jl docs (so I went to the code and found it had a docstring). It should really either be documented here, or moved out to LsqFit and documented there.

That is pretty much why I opened #207, but I guess that's still unresolved.

[pkofod]

@ChrisRackauckas maybe repost @ https://github.com/JuliaOpt/LsqFit.jl/issues ? You can argue that we might want to "backport" any solution you come up with to Optim while we're in the deprecation period. @matthieugomez @bjarthur