Further ODE speedup opportunity with reduced sensitivity system?

[Sebastian Weber]

Hi!

I had an idea which seems to turn out very attractive for complex ODE modeling. In PKPD modeling and I believe generally in complex structural ODE modeling where Bayes is very attractive to use due to available physical knowledge of the system, the situation is very often that quite a number of ODEs N are needed to describe the system, but only very few No << N are observed and enter the likelihood.

In my example I am considering, I have a 2-cmt oral dosing PK system which has N=3, but only concentrations from the second compartement are observed and enter the likelihood. Currently Stan calculates the sensitivities with respect to parameters and initials for all N equations - but only Ns (=1 in my case) enter the likelihood. So the full coupled system in this case (M=4 parameters) is 24. Now, I thought why not reduce the coupled system to only what is needed, which is in this case 7 sensitivities for the second equation and 3 for the base system, so in total only 10 equations to solve.

I applied this idea to the cvode branch I worked out with the user-specified sensitivities and the results I get as parameters estimates are exactly the same. The theoretical speedup should be 2.4, but I end up getting 1.7:

Full coupled system, 24 ODEs (250 warmup/250 iterations):
# Elapsed Time: 207.61 seconds (Warm-up)
# 76.71 seconds (Sampling)
# 284.32 seconds (Total)

Reduced coupled system, 10 ODEs (same iterations):
# Elapsed Time: 122.37 seconds (Warm-up)
# 46.28 seconds (Sampling)
# 168.65 seconds (Total)

The advantage of this approach grows quadratic in N! I know that many modelers love to come up with tons of euqations to model their system, but only very few of those quantities are observed such that this (little) optimization will have a large impact, I think. It would enable to apply Stan to very complex systems with few observed outcomes - a typical situation from what I am aware of.

Best,
Sebastian

[Sebastian Weber]

Small addition: The idea to reduce the sensitivity system can also be applied to the current implementation (without user specified sensitivities, but it will likely need a user-interface addition or change). In fact, I would guess that the speedup under the AD based system now is closer to the theoretical speedup, since currently getting the sensitivities relies on auto-diffing the user-supplied base system - which happens inside the system functor (resulting in many repetitions of this process). Hence, the cost for calculating the sensitivities with AD is probably large such that removing many of those unnecessary computations will speedup considerably. Currently we always have the computation cost proportional to N + N*(N+M) which will be reduced to N + No * (N+M); sorry for my typo above where I wrote Ns=1, I meant No=1.

Sebastian

[Bob Carpenter]

I need very simple examples to understand this stuff. Is there a simple system
illustrating the problem?

Is this just a data vs. parameters distinction
if you don't need sensitivities? Otherwise, I think you need
sensitivities to all of the inputs that can affect the outputs.

- 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]

Ok, I can sent you what I have compiled. This is working code, but on a prototype level; I hope that's fine.

My empirical experiments seem to support what I have done, but is there some math written up and documented which shows what terms are exactly needed for NuTS?

Sebastian

On Monday, August 31, 2015 at 10:26:18 PM UTC+2, Bob Carpenter wrote:
> I need very simple examples to understand this stuff. Is there a simple system
> illustrating the problem?
>
> Is this just a data vs. parameters distinction
> if you don't need sensitivities? Otherwise, I think you need
> sensitivities to all of the inputs that can affect the outputs.
>
> - Bob
>
> > On Aug 29, 2015, at 5:26 AM, Sebastian Weber
> >

[Bob Carpenter]

Yes, the autodiff paper explains how things are done.
We're about to release it when we release 2.7.1, but
I can send it to you (Sebastian) before that.

- Bob