Robust RK45 ODE integrator

[Sebastian Weber]

Hi!

Currently, the use of the ODE integrator is very limited as it may freeze easily Stan whenever a unlucky parameter is tried out which renders the ODE system too stiff to attain an accuracy as requested.

That is, odeint does not fail in these situations which is what would make the system robust. I figured that we were unlucky with the choice of the integration scheme in odeint. That is, the dense_output scheme does not have any means to protect against this situation. Now, the controlled_stepper scheme does protect against this! That is, after 1000 tries, the controlled_stepper will throw an exception which gives Stan the required robustness.

Issue: https://github.com/stan-dev/math/issues/179
Branch with fix for it: https://github.com/stan-dev/math/tree/feature/issue-179-robust-integrator

I tested this on a 2-cmt oral dosing system which did kill Stan in iteration 3 so far with the old scheme. Now, the integrator simply rejects a few draws at the beginning, but then Stan continues sampling.

I hope this can get in quickly... as now people do not freeze Stan anymore.

Sebastian

[Bob Carpenter]

Yes, we very much want to get this in for the next release.
It's a huge feature just to get Stan not to hang. It
needs to raise the right kind of exception so that it
gets recorded as a divergence, but otherwise, not many restrictions.

Did you test it for time? Mine and Daniel's understanding was
that the interpolation was important if you had multiple solutions.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Sebastian Weber]

Hi Bob!

Great that this goes in soon. I did only some brief performance comparisons which showed performance in the same ballpark for the tests. However, my biggest interest was in a model which caused Stan to hang and with the patch things run with a couple of exceptions through.

Not sure what you refer to with multiple solutions. My understanding is that the dense output has some advantages to use if the step_size dt can be chosen very large, i.e. large here with respect to the density of points. The controlled stepper increases step sizes between time-points when you request the solution, but when it approaches a point where the solution is requested, the controlled steppers choses dt such that the stepper "steps" on the time-point when the solution is needed. Hence for functions which can have large dt (straight lines) and very dense sampling scheme t_{j+1} - t_j is at the order of a typical dt, then the dense output is better. Otherwise, there should´t be a noticeable difference which I think is true for many applications for which Stan is used.

I just followed an example of odeint which combines the dense output explicitly with the controlled range kuttta - maybe this gives us both (robust & dense) properties, I have to test it tomorrow on the example which causes Stan to hang.

BTW, have a look here:

https://github.com/headmyshoulder/odeint-v2/blob/master/include/boost/numeric/odeint/integrate/detail/integrate_times.hpp

There you see that the controlled stepper has an explicit statement about too many tries which triggers an exception whenever 1000 tries are not sufficient to get the solution with the desired accuracy.

Best,
Sebastian

[Sebastian Weber]

Hi Bob!

I think I understood the problem here. This is possibly a bug of odeint. I filed an issue

https://github.com/headmyshoulder/odeint-v2/issues/173

Let's see what happens, but I am afraid the only way to get a robust integrator with odeint right now is with the controlled error stepper. The issue above explains my conclusion from what is happening; but maybe we are using the library just in a wrong way. Let's see.

Best,
Sebastian

[Sebastian Weber]

Hi!

After some long discussions with the odeint author (many thanks to him!!!), I can now propose a solution which is straightforward and does not rely on side-effects which the current solutions relies on.

In essence: One wants to limit the total number of calls to the system functor. So if the system functor is called excessively, then an exception must be raised. This means to count the number of evaluations of the coupled ode system functor and throw an exception whenever this exceeds a very big number. This is according to the odeint author the appropriate way of doing things and it will still allow us to use the dense output scheme.

Sebastian

[Daniel Lee]

That sounds reasonable!

Can you propose a number of iterations as a starting point? We can move forward from there.



Daniel

[Sebastian Weber]

See what I just pushed, i.e. a total of 1e5 max. However, there is still a problem in that the total number of steps is limited now - but you want to limit the number of steps between two time-points where the user requests the solution. This is what CVODE does - they limit by default the number of steps between two consecutive steps to 500 by default. To make this very conservative a value of 1e3 would be good, I think. Now that means the current proposal needs to be changed by either

1) giving the observer object a reference to the coupled ode system such that the observer can reset the number of counted evaluations (clean, but creates dependencies)

2) or, alternatively, you can make the number of iterations when the exception is thrown depend on the number of time-points the user wants the solution (less dependencies, but one more constructor argument for coupled_ode_system).

Sebastian

BTW, the erase operation after the integrate_times could easily be avoided (it is of order O(size of vector ts) ) by instructing the coupled ode observer to always discard the very first call to it.

On Tuesday, October 6, 2015 at 1:11:00 PM UTC+2, Daniel Lee wrote:
> That sounds reasonable!
>
> Can you propose a number of iterations as a starting point? We can move forward from there.
>
>
>
> Daniel
>
>
> > On Oct 6, 2015, at 5:26 AM, Sebastian Weber
> >

[Sebastian Weber]

Hi!

Thinking about it a bit more, we should consider going with a maximal density of steps. That would mean to set the maximal number of steps allowed to

(t_max - t_min)/step_size * 1E3

This would say we allow for 1E3 steps per step_size time-units. Makes sense?

Sebastian

[Bob Carpenter]

It's just a heuristic. Nothing's going to make
sense everywhere because the scale of the function's
smoothness relative to the time scale isn't something
we control. I'd think maybe 1E3 calls per solution
time would be a more reasonable default. But it sounds
like something we need to control, because that could be
too low for some cases.

- Bob

[Sebastian Weber]

Of course, it is a heuristic and it should be controlled. However, since abs+rel tolerance are fixed this determines to some degree what smoothes you expect, i.e. setting step_size to sqrt(abs_tol) maybe is a very safe choice? To my understanding you designed the ODE function such that 1 is a unit which is reasonable (a good modeling principle).

Sebastian