Definitely agree there. That said at this point, seeing more descriptions of RK --> Rosenbrock and Adams --> BDF, so I’m wondering whether those flows canonically work more easily somehow and whether the two-way stiffness detection works more easily for those families.
Dumb question – would you happen to know how to estimate Jacobian eigenvalues at each step for the Boost algorithms? Are eigenvalue estimates already optionally output in each step as a part of Boost (guessing not)? For systems with < 30 variables, I think we can use stiffness detection to look for trouble regions in RK, and as you said, then one-way jump to CVODE BDF as long as the norm bounding process is straight forward.
Haha, yeah, felt nervous asking just in case it was an obvious thing, but glad it was not. Yeah, I think that’s the next step to figure out. A lot more to be done, but it looks like things are much more straightforward if there is a way of extracting the eigenvalues from Boost.
I don’t want to stop anyone from putting work into this… but I really think that adjoint sensitivity analysis is the much more important thing to tackle if we ever want to solve larger systems. Speed problems in smaller ODE systems can be dealt with whenever we talk about a hierarchical problem (using brute-force reduce_sum parallelization)…but large ode systems are a killer for Stan right now.
For sure, this is definitely very secondary to the adjoint sensitivity priority, a side project if you will. I can iteratively plug at the switching algorithm with minimal bothering of the devs so that y’all can focus on adjoint sensitivity.
@bbbales2 Having done a little bit more digging, it looks like the Boost-to-CVODE plan won’t be possible because stiffness detection is not included within Boost’s odeint. Also, going between libraries will also be tricky if there are a few of arguments being passed between (for example, in the case of a forcing vector). For instance, the remaining portion of a forcing vector would have to be passed through, which could be a nuisance.
So, there’s a couple of ways to go from here regarding a switching algorithm:
Wrapping a C-based pre-CVODE LSODA solver such as liblsoda (https://github.com/sdwfrost/liblsoda) that uses the aging norm-based dominant eigenvalue stiffness detection method for an intermediate solution to switch between Adams and BDF.
Implementing a separate RK - Rosenbrock switching method natively in Stan Math. There are some newer RK methods available (like Tsit5 and Vern 6+). We would have to go this route if we want something switching from RK. According to some in numerical analysis, the RK paradigm is still overall the best default family of non-stiff methods. I’m willing to do some more research into this as an exercise for myself. This route affords us more flexibility with updated means of stiffness and non-stiff detection. (Would one of you guys happen to know what the workflow would be for building a solver in Stan Math, like which folders the solver files would go into, etc.?)
Any thoughts on which route to go for there?
Again, I don’t want to step on any toes, and clearly, the adjoint sensitivity detection is the biggest ODE priority. This would be something I’d be up for working together with people on down the line. For the time being, before I get to a phase in my project where speed really is an issue, I’m going to stick with old trusty CVODE BDF for testing and trialing purposes.
Great, thanks for the pointer! I’ll read through those examples. We’ll go ahead and sketch in a tentative side project plan for an RK --> Rosenbrock then (again, apologies in advance for the slow speed on getting work done on this). Linking to the Shampine’s 1982 paper that describes the RK --> Rosenbrock algorithm: http://www.hysafe.org/science/eAcademy/docs/ACMTransactionsOnMathematicalSoftware_v8_p93to113.pdf
@xiehw, @wds15 I am a new Stan user. I was curious if number 6 and 7 from the original priority list were ever implemented in stan and if so is can I please be pointed towards the relevant documentation. Thanks.
Addressing another thread Handling events when using ode integrator in the same time, I’d like to point out that both item 6 & 7 are possible. Item 6 is doable by manually setting integration boundaries. Item 7 can only be done through our backend ODE integrators that support rootfinding as we haven’t exposed it yet.
By that logic, one would have a parameter say t_interupt where at that right before it state variables are exposed to the desired events and then put back into the integrator ?
Except, how does one guarantee that the time an event is desired to happen is available, since the underlying integrator might not access the specified time
There’s no such guarantee in general, much as most numerical algorithms that requires user discretion and prior knowledge. A simple example is when the event function has multiple roots and the integrator may only find one or even none of them. It eventually comes down to the user to apply right hammer to right nail.