CMAE ES available in Julia

[finmod]

Following the recent exchange https://github.com/JuliaDiffEq/DiffEqParamEstim.jl/issues/60#issuecomment-356646400, I have uploaded here https://github.com/finmod/DiffEqParamEstim.jl/blob/master/test/lorenz%20true%20test.ipynb my experience of estimating from data the parameters of the Lorenz system with a range of optimizers including BBO and NLopt families. In this version, the issue of accuracy of the solver has now been solved by adopting Vern9. This has shifted the spotlight to the efficiency of the optimizers and BBO remains the main tool to recover the parameters when you dont know anything about the starting values (at least positive). BBO is very slow to converge after finding rather quickly a narrow range for the parameters. A practical solution would be BBO + Local optimizer.

From your benchmarking exercise, CMAE-ES appears to be an interesting challenger to BBO in the search for a global optimum. Which version of CMAE-ES in Julia do you advise trying?

[jonathanBieler]

Hey, I don't think there's a robust CMA-ES available in Julia right now, I'm working on one, so you can give it a try, but it's very rough.

https://gist.github.com/jonathanBieler/456b37c4cbae70d269ae6899058092ba

For 3 parameters you probably want a population size (λ) a bit smaller (I would try 12 maybe).

Otherwise you can call the python version using PyCall, see here:

https://github.com/jonathanBieler/BlackBoxOptimizationBenchmarking.jl/blob/master/scripts/optimizers_interface.jl#L32

You also need to have the package installed in python:

https://pypi.python.org/pypi/cma

[finmod]

@ChrisRackauckas Please advise on the sequence:

• PR to put https://github.com/finmod/DiffEqParamEstim.jl/blob/master/test/lorenz%20true%20test.ipynb on DiffEqBenchmarks and then insert CMAE-ES in the family of available optimizers or
• sort out first the best CMAE-ES in Julia between cmaes.jl, PyCMA and other sources compatible with JuliaDiffEq ecosystem. Then insert it in the benchmarks notebook?

[ChrisRackauckas]

Either way is fine. I'd prefer the first since we might as well get those benchmarks out, and then add to them. A lot of our benchmarks get regularly updated so this won't be out of the ordinary.

[finmod]

@ChrisRackauckas I notice that CMAE-ES is also available at Evolutionary.jl. It is in one of your runtests somewhere?

[ChrisRackauckas]

Yes but it needs master since Evolutionary needs a tag.

[ChrisRackauckas]

Master on Evolutionary.jl

[finmod]

@ChrisRackauckas The benchmark for parameter estimation of the lorenz system from a dataset is available here: https://github.com/finmod/DiffEqParamEstim.jl/blob/master/test/parameter%20estimation%20benchmark%20for%20lorenz%20system.ipynb

Can you please check that the set up of the data and model are correct because I got lost in the latest updates with saveat, callback etc. In short, the model should be the same as @Datseris to avoid further confusions with prototype models and datasets.

On the optimizers, the default BBO adaptive_de_rand_1_bin_radiuslimited() remains best for now. Julia versions of CMAES from Evolutionary and others fail to impress. pyCMA should be tried here but I have to solve how to install cma in Conda.jl with Conda.add(). Hence, I would not expect markedly different results from @jonathanBieler until an efficient Julia version of CMAES can match pyCMA.

The estimation problem should move from least square (minimum distance) to likelihood (minimum variance) to increase the speed of convergence. A simple interface to EM, MAP and MLE estimators is needed. Suggestions and solutions can be sought from @genkuroki using Distributions, or from @tpapp with DynamicHMC associating an auxiliary Bayesian model to the dataset or from @Datseris on whether restructure() may allow a recovery of the parameters from the Lyapunov's exponents.

[Datseris]

Hello. I need some explanation on why I was tagged here, as I don't really know what to do :P

[ChrisRackauckas]

You should make the first tspan equal to 30. That's why it plots a little weird, since it continues to solve until 50 and then stops (and saves the last point). I don't think it effects the optimization though. Other than that, the setup looks right. Would you mind PRing that to DiffEqBenchmarks.jl?

On the optimizers, the default BBO adaptive_de_rand_1_bin_radiuslimited() remains best for now. Julia versions of CMAES from Evolutionary and others fail to impress. pyCMA should be tried here but I have to solve how to install cma in Conda.jl with Conda.add(). Hence, I would not expect markedly different results from @jonathanBieler until an efficient Julia version of CMAES can match pyCMA.

I can't wait to see it. From @jonathanBieler 's benchmarks it seems it may be way better than what we've been using so far.

The estimation problem should move from least square (minimum distance) to likelihood (minimum variance) to increase the speed of convergence.

This is the same exact thing. I've mentioned this. Least square minimization in this case is equivalent to maximum likelihood estimation for the problem where you assume that the error is homoscedastic normally-distributed measurement error. The heteroscedastic case is equivalent to using a weight vector in the cost function. Then you can also add regularization to the cost function. You haven't played with the weight vectors or regularizations yet though. Weight vectors don't make sense in that test case because we generated it with homoscedastic noise though. Regularization would probably help and should be tested.

A simple interface to EM, MAP and MLE estimators is needed.

Suggestions and solutions can be sought from @genkuroki using Distributions

Distributions are being tracked here (https://github.com/JuliaDiffEq/DiffEqParamEstim.jl/issues/39), but in the optimization sense they are just equivalent to using different metrics for the cost function. In the Bayesian sense they are more interesting. Using Distributions.jl in the optimization sense then only makes sense for stochastic problems because that's the only case where you'd have a meaningful sample of the likelihood, which is what's discussed there.

A simple interface to EM, MAP and MLE estimators is needed.

EM isn't possible for this because of the nonlinear transformation. SAEM (a stochastic approximation of the likelihoods) is, and should be considered (good paper). MAP is a good idea too with a nice paper. And as above, MLE is the optimization approach. I don't think there's good systematic benchmarking of what does well in the ODE case though, so we should implement them all and find out. I started an issue for it: (https://github.com/JuliaDiffEq/DiffEqParamEstim.jl/issues/62) and will write this up as a possible Google Summer of Code project (both have nice papers, should be a straightforward student project to implement and benchmark).

What you're missing though are the ODE-specific methods. Did you try the two-stage method to see how close it gets? Also, we need multiple shooting methods.

or from @tpapp with DynamicHMC associating an auxiliary Bayesian model to the dataset

See DiffEqBayes.jl (https://github.com/JuliaDiffEq/DiffEqBayes.jl/issues/16). We will likely have a GSoC for this. That sets it up with DynamicHMC.jl, Stan.jl, and Turing.jl. Don't look at this as a magical better method though since it's generally more expensive so the optimization approach is much much faster.

[ChrisRackauckas]

Hello. I need some explanation on why I was tagged here, as I don't really know what to do :P

Do you have any parameter estimation for chaotic systems stuff? Just curious.

[Datseris]

I do not have anything like that in the DynamicalSystems.jl ecosystem. Even though parameter estimation seems to be really closely related to dynamical systems, it is really far of from what I know as "nonlinear dynamics and chaos" which is where my packages are aimed.

I also have not used it in (personal) research. As far as I know, parameter estimation is not really a big thing in nonlinear dynamics and chaos communities. It is a big thing in modelling for sure, but I have not seen or heard a talk or paper on chaotic systems that references it directly... Of course, I am a very young scientist so my word doesn't really have any weight behind it.

Never-the-less, regardless of my familiarity with the overall "vibe" of parameter estimation, it is clear to me that it can help anybody interested in studying dynamical systems. Hopefully with the latest update that allows a ContinuousDS to be directly related with its ODEProblem people can take full advantage of all parameter-estimation methods offered by DifferentialEquations. I think the only thing left is to allow/be able to handle ParameterizedFunctions for the parameter handling functions of my packages (like e.g. orbit diagrams).

Or make the parameter estimation methods work with normal functors :P

[finmod]

EM isn't possible for this because of the nonlinear transformation. SAEM (a stochastic approximation of the likelihoods) is, and should be considered (good paper). MAP is a good idea too with a nice paper. And as above, MLE is the optimization approach. I don't think there's good systematic benchmarking of what does well in the ODE case though, so we should implement them all and find out. I started an issue for it: (JuliaDiffEq/DiffEqParamEstim.jl#62) and will write this up as a possible Google Summer of Code project (both have nice papers, should be a straightforward student project to implement and benchmark).

Let me add this paper https://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjsje-owNzYAhUIyKQKHbRjDhsQFggxMAA&url=https%3A%2F%2Farxiv.org%2Fabs%2F1609.03508&usg=AOvVaw0zy9BIJnaJ_yy2vrbCZzqr on synthetic likelihood SAEM. There is a corresponding MATLAB code https://github.com/umbertopicchini/SAEM-SL and an application to LV ODE mentioned in the paper. @umbertopicchini

[jonathanBieler]

@finmod The CMA-ES I linked above should have very similar performance as the python one, the issues are more with the interface, options and error handling, but I have to wait for some changes in Optim to improve that.

For python I found easier to point PyCall to my normal python distribution and install packages there with pip.

[finmod]

In fact, I am trying to run cmaes.jl on the bench. For this input:

optimizers = [ OptimRestart(NelderMead()), CMAESoptim, NelderMead(), BlackBoxOptimMethod(:adaptive_de_rand_1_bin_radiuslimited), BlackBoxOptimMethod(:adaptive_de_rand_1_bin),

NLoptOptimMethod(:LN_BOBYQA), # fails too

#NLoptOptimMethod(:LN_PRAXIS), #this guy often fails
#NLoptOptimMethod(:LN_SBPLX), # this guy segfault
BlackBoxOptimMethod(:xnes),
BlackBoxOptimMethod(:generating_set_search),
BlackBoxOptimMethod(:de_rand_2_bin),
#BlackBoxOptimMethod(:resampling_memetic_search),#poor performance 
#PyMinimize("Nelder-Mead"),
#PyCMA(),

]

opt_strings = map(string,optimizers)

I get he following error:

WARNING: using BlackBoxOptim.Parameters in module Main conflicts with an existing identifier. WARNING: Method definition minimum( MethodError: no method matching string(::TypeName) You may have intended to import Base.string Closest candidates are: string(::BlackBoxOptimMethod) at C:\Users\Denis.julia\v0.6\BlackBoxOptimizationBenchmarking\scripts\optimizers_interface.jl:70 string(::NLoptOptimMethod) at C:\Users\Denis.julia\v0.6\BlackBoxOptimizationBenchmarking\scripts\optimizers_interface.jl:12 string(::Optim.Optimizer) at C:\Users\Denis.julia\v0.6\BlackBoxOptimizationBenchmarking\scripts\optimizers_interface.jl:48 ...

Stacktrace: [1] string(::OptimRestart{Optim.NelderMead{Optim.AffineSimplexer,Optim.AdaptiveParameters}}) at C:\Users\Denis.julia\v0.6\BlackBoxOptimizationBenchmarking\scripts\optimizers_interface.jl:63 [2] _collect(::Array{Any,1}, ::Base.Generator{Array{Any,1},#string}, ::Base.EltypeUnknown, ::Base.HasShape) at .\array.jl:488 [3] map(::Function, ::Array{Any,1}) at .\abstractarray.jl:1868 [4] include_string(::String, ::String) at .\loading.jl:522

Optim.OptimizationResults) in module Optim at C:\Users\Denis.julia\v0.6\Optim\src\api.jl:3 overwritten in module Main at C:\Users\Denis.julia\v0.6\BlackBoxOptimizationBenchmarking\scripts\optimizers_interface.jl:45. WARNING: Method definition minimizer(Optim.OptimizationResults) in module Optim at C:\Users\Denis.julia\v0.6\Optim\src\api.jl:2 overwritten in module Main at C:\Users\Denis.julia\v0.6\BlackBoxOptimizationBenchmarking\scripts\optimizers_interface.jl:46.